inst/doc/Seq-vignette.R

In sequential.pops: Sequential Analysis of Biological Population Sizes

## ----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)