# Function that describes negative binomial overdisperion from the mean
# and Taylor's Power Law parameters, a and b:
estimate_k <- function(mean) {
a = 1.830012
b = 1.218041 # a and b are Taylor's Power Law parameters
(mean^2) / ((a * mean^(b)) - mean)
}
# Generate a SPRT object with negative binomial and varying overdispersion
counts <- c(2, 5, 6, 2, 7)
test11 <- sprt(data = counts,
mu0 = 8,
mu1 = 10,
density_func = "negative binomial",
overdispersion = estimate_k(9),
alpha = 0.1,
beta = 0.1)
# Stochastic version of 'estimate k', where sd here is the the root of the
# mean square error for the regression used to estimate a and b.
estimate_k_stoch <- function(mean) {
a <- 1.830012
b <- 1.218041
(mean^2) /
((a * mean^(b) *
exp(truncdist::rtrunc(1,
"norm",
a = log(1 / (a * mean^(b - 1))),
b = Inf,
mean = 0,
sd = 0.3222354)))
- mean)
}
# Run model evaluation for test11 with varying overdispersion and
# added stochasticity.
eval4 <- SPRT.eval(test11, eval.range = seq(4, 13),
overdispersion.sim = "estimate_k_stoch", N = 30)
plot(seq(4, 13), eval4$AvgSamples, type = "o", xlab = "True population size",
ylab = "Average number of bouts")
plot(seq(4, 13), eval4$AcceptRate, type = "o", xlab = "True population size",
ylab = "Acceptance rate of H1")
## End (Not run)
Run the code above in your browser using DataLab