## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) set.seed(000) ## ----setup-------------------------------------------------------------------- a30 <- matrix(rep(0, 120), 30, 4) head(a30, 3) # only first three rows out of 30 are shown ## ----------------------------------------------------------------------------- library(sequential.pops) test30 <- stbp_simple(data = a30, density_func = "poisson", prior = 0.5, upper_bnd = Inf, lower_criterion = 0, upper_criterion = 0.9999) test30 ## ----------------------------------------------------------------------------- a3 <- matrix(rep(0, 30), 3, 10) a3 ## ----------------------------------------------------------------------------- test3 <- stbp_simple(data = a3, density_func = "poisson", prior = 0.5, upper_bnd = Inf, lower_criterion = 0, upper_criterion = 0.9999) test3 ## ----------------------------------------------------------------------------- plot(test30) ## ----------------------------------------------------------------------------- plot(test3) ## ----------------------------------------------------------------------------- counts1 <- rnbinom(20, mu = 5, size = 4.5) counts1 ## ----------------------------------------------------------------------------- test_counts1 <- stbp_composite(data = counts1, greater_than = TRUE, hypothesis = 9, density_func = "negative binomial", overdispersion = 4.5, prior = 0.5, lower_bnd = 0, upper_bnd = Inf, lower_criterion = 0.01, upper_criterion = 0.99) test_counts1 ## ----------------------------------------------------------------------------- estimate_k <- function(mean) { a = 1.83 b = 1.21 (mean^2) / ((a * mean^(b)) - mean) } ## ----------------------------------------------------------------------------- counts1a <- rnbinom(20, mu = 5, size = estimate_k(5)) counts1a ## ----------------------------------------------------------------------------- test_counts1a <- stbp_composite(data = counts1a, greater_than = TRUE, hypothesis = 9, density_func = "negative binomial", overdispersion = "estimate_k", prior = 0.5, lower_bnd = 0, upper_bnd = Inf, lower_criterion = 0.01, upper_criterion = 0.99) test_counts1a ## ----------------------------------------------------------------------------- plot(test_counts1) ## ----------------------------------------------------------------------------- plot(test_counts1a) ## ----------------------------------------------------------------------------- library(emdbook) # for rbetabinom() binom1 <- list() for(i in 1: 7) { binom1[[i]] <- matrix(c(emdbook::rbetabinom(3, prob = 0.2, size = 10, theta = 6.5), rep(10, 3)), 3, 2) } head(binom1, 3) # only first three sampling bouts out of seven are shown. ## ----------------------------------------------------------------------------- test_bin1 <- stbp_composite(data = binom1, greater_than = TRUE, hypothesis = 0.15, density_func = "beta-binomial", overdispersion = 6.5, prior = 0.5, lower_bnd = 0, upper_bnd = 1, lower_criterion = 0.001, upper_criterion = 0.999) test_bin1 ## ----------------------------------------------------------------------------- estimate_theta <- function(mean) { A = 780.72 b = 1.61 R = 10 (1 / (R - 1)) * (((A * R^(1 - b))/ ((mean * (1 - mean))^(1 - b))) - 1) } ## ----------------------------------------------------------------------------- binom2 <- list() for(i in 1: 7) { binom2[[i]] <- matrix(c(rbetabinom(3, prob = 0.2, size = 10, theta = estimate_theta(0.2)), rep(10, 3)), 3, 2) } head(binom2, 3) # only first three sampling bouts out of seven are shown. ## ----------------------------------------------------------------------------- test_bin2 <- stbp_composite(data = binom2, greater_than = TRUE, hypothesis = 0.15, density_func = "beta-binomial", overdispersion = "estimate_theta", prior = 0.5, lower_bnd = 0, upper_bnd = 1, lower_criterion = 0.001, upper_criterion = 0.999) test_bin2 ## ----------------------------------------------------------------------------- eval1 <- STBP.eval(test_counts1, eval.range = seq(3, 12), n = 1, prior = 0.5, N = 20) eval1 ## ----------------------------------------------------------------------------- plot(seq(3, 12), eval1$AvgSamples, type = "o", xlab = "True population size", ylab = "Average number of bouts") abline(v = 9, lty = 2) plot(seq(3, 12), eval1$AcceptRate, type = "o", xlab = "True population size", ylab = "Acceptance rate of H") abline(v = 9, lty = 2) ## ----------------------------------------------------------------------------- estimate_k_stoch <- function(mean) { a = 1.83 b = 1.21 RMSE = 0.32 (mean^2) / ((a * mean^(b) * exp(truncdist::rtrunc(1, "norm", a = log(1 / (a * mean^(b - 1))), b = Inf, mean = 0, sd = RMSE))) - mean) } ## ----------------------------------------------------------------------------- eval1a <- STBP.eval(test_counts1, eval.range = seq(3, 12), n = 1, prior = 0.5, overdispersion.sim = "estimate_k_stoch", N = 20) eval1a ## ----------------------------------------------------------------------------- eval3a <- STBP.eval(test_counts1, eval.range = seq(3, 12), n = 3, prior = 0.5, overdispersion.sim = "estimate_k_stoch", N = 20) eval6a <- STBP.eval(test_counts1, eval.range = seq(3, 12), n = 6, prior = 0.5, overdispersion.sim = "estimate_k_stoch", N = 20) ## ----------------------------------------------------------------------------- plot(seq(3, 12), eval1a$AvgSamples, type = "o", xlab = "True population size", ylab = "Average number of bouts") points(seq(3, 12), eval3a$AvgSamples, type = "o", pch = 19) points(seq(3, 12), eval6a$AvgSamples, type = "o", pch = 0) abline(v = 9, lty = 2) plot(seq(3, 12), eval1a$AcceptRate, type = "o", xlab = "True population size", ylab = "Acceptance rate of H") points(seq(3, 12), eval3a$AcceptRate, type = "o", pch = 19) points(seq(3, 12), eval6a$AcceptRate, type = "o", pch = 0) abline(v = 9, lty = 2) ## ----------------------------------------------------------------------------- estimate_theta_stoch <- function(mean) { A = 780.72 b = 1.61 R = 10 RMSE = 0.41 ((1 / (R - 1)) * (((A * R^(1 - b) * exp(rnorm(1, mean = 0, sd = RMSE))) / ((mean * (1 - mean))^(1 - b))) - 1)) } ## ----------------------------------------------------------------------------- evalB10 <- STBP.eval(test_bin2, eval.range = seq(0.01, 0.25, 0.02), n = 10, prior = 0.5, overdispersion.sim = "estimate_theta_stoch", N = 20) evalB20 <- STBP.eval(test_bin2, eval.range = seq(0.01, 0.25, 0.02), n = 20, prior = 0.5, overdispersion.sim = "estimate_theta_stoch", N = 20) evalB30 <- STBP.eval(test_bin2, eval.range = seq(0.01, 0.25, 0.02), n = 30, prior = 0.5, overdispersion.sim = "estimate_theta_stoch", N = 20) ## ----------------------------------------------------------------------------- plot(seq(0.01, 0.25, 0.02), evalB10$AvgSamples, type = "o", xlab = "True incidence", ylab = "Average number of bouts") points(seq(0.01, 0.25, 0.02), evalB20$AvgSamples, type = "o", pch = 19) points(seq(0.01, 0.25, 0.02), evalB30$AvgSamples, type = "o", pch = 0) abline(v = 0.15, lty = 2) plot(seq(0.01, 0.25, 0.02), evalB10$AcceptRate, type = "o", xlab = "True incidence", ylab = "Acceptance rate of H") points(seq(0.01, 0.25, 0.02), evalB20$AcceptRate, type = "o", pch = 19) points(seq(0.01, 0.25, 0.02), evalB30$AcceptRate, type = "o", pch = 0) abline(v = 0.15, lty = 2) ## ----------------------------------------------------------------------------- pop.outB <- function(t) { N0 = 15 K = 90 r = 0.25 K / (1 + (K / N0 - 1) * exp(-r * t)) } OB_pop <- pop.outB(t = seq(1, 20)) plot(seq(1, 20), OB_pop, xlab = "Time (weeks)", ylab = "Population density", lwd = 2, type = "o", pch = 19) ## ----------------------------------------------------------------------------- plot(seq(1, 20), OB_pop, ylim = c(0, 95), type = "o", lwd = 2, ylab = "Population size", xlab = "Time (weeks)", pch = 19) points(seq(1, 20), OB_pop * 0.9, type = "o") points(seq(1, 20), OB_pop * 0.7, type = "o", pch = 0) points(seq(1, 20), OB_pop * 0.5, type = "o", pch = 2) points(seq(1, 20), OB_pop * 1.1, type = "o", pch = 5) ## ----------------------------------------------------------------------------- pop70 <- OB_pop * 0.7 count70 <- matrix(NA, 10, 20) for(i in 1:20){ count70[, i] <- rnbinom(10, mu = pop70[i], size = estimate_k(pop70[i])) } head(count70, 3) # only first three samples out of ten are shown ## ----------------------------------------------------------------------------- plot(seq(1, 20), OB_pop, ylim = c(0, 95), type = "o", lwd = 2, ylab = "Population size", xlab = "Time (weeks)", pch = 19) points(seq(1, 20), colMeans(count70), type = "o") ## ----------------------------------------------------------------------------- test_dyn <- stbp_composite(data = count70, greater_than = TRUE, hypothesis = OB_pop, density_func = "negative binomial", overdispersion = "estimate_k", prior = 0.5, lower_bnd = 0, upper_bnd = Inf, lower_criterion = 0.001, upper_criterion = 0.999) test_dyn ## ----------------------------------------------------------------------------- eval_dyn <- STBP.eval(test_dyn, eval.range = seq(0.2, 1.5, 0.2), n = 10, prior = 0.5, N = 20) ## ----------------------------------------------------------------------------- plot(seq(0.2, 1.5, 0.2), eval_dyn$AvgSamples, type = "o", xlab = "True outbreak intensity", ylab = "Average number of bouts") abline(v = 1, lty = 2) plot(seq(0.2, 1.5, 0.2), eval_dyn$AcceptRate, type = "o", xlab = "True outbreak intensity", ylab = "Acceptance rate of H") abline(v = 1, lty = 2) ## ----------------------------------------------------------------------------- testPR10 <- sprt(mu0 = 8, mu1 = 10, density_func = "negative binomial", overdispersion = estimate_k(9), alpha = 0.1, beta = 0.1) ## ----------------------------------------------------------------------------- plot(testPR10) ## ----------------------------------------------------------------------------- testPR10 ## ----------------------------------------------------------------------------- counts1a <- rnbinom(20, mu = 5, size = estimate_k(5)) counts1a ## ----------------------------------------------------------------------------- testPR1 <- sprt(data = counts1a, mu0 = 8, mu1 = 10, density_func = "negative binomial", overdispersion = estimate_k(9), alpha = 0.1, beta = 0.1) testPR1 ## ----------------------------------------------------------------------------- plot(testPR1) ## ----------------------------------------------------------------------------- binPR <- rbinom(40, prob = 0.5, size = 1) binPR ## ----------------------------------------------------------------------------- testPR20 <- sprt(mu0 = 0.35, mu1 = 0.45, density_func = "binomial", alpha = 0.1, beta = 0.1) plot(testPR20) ## ----------------------------------------------------------------------------- testPR20 ## ----------------------------------------------------------------------------- testPR2 <- sprt(data = binPR, mu0 = 0.35, mu1 = 0.45, density_func = "binomial", alpha = 0.1, beta = 0.1) testPR2 ## ----------------------------------------------------------------------------- plot(testPR2) ## ----------------------------------------------------------------------------- evalPR1 <- SPRT.eval(testPR1, eval.range = seq(4, 12), N = 20) ## ----------------------------------------------------------------------------- evalPR1a <- SPRT.eval(testPR1, eval.range = seq(4, 12), overdispersion.sim = "estimate_k_stoch", N = 20) ## ----------------------------------------------------------------------------- plot(seq(4, 12), evalPR1a$AvgSamples, xlab = "True population size", ylab = "Average number of bouts", type = "o") points(seq(4, 12), evalPR1$AvgSamples, type = "o", pch = 2) abline(v = 9, lty = 2) plot(seq(4, 12), evalPR1a$AcceptRate, xlab = "True population size", ylab = "Acceptance rate of H1", type = "o") points(seq(4, 12), evalPR1$AcceptRate, type = "o", pch = 2) abline(v = 9, lty = 2) ## ----------------------------------------------------------------------------- evalPR2 <- SPRT.eval(testPR2, eval.range = seq(0.1, 0.8, 0.1), N = 20) plot(seq(0.1, 0.8, 0.1), evalPR2$AvgSamples, xlab = "True incidence", ylab = "Average number of bouts", type = "o") abline(v = 0.4, lty = 2) plot(seq(0.1, 0.8, 0.1), evalPR2$AcceptRate, xlab = "True incidence", ylab = "Acceptance rate of H1", type = "o") abline(v = 0.4, lty = 2)