|
#--------------------------------# |
|
# Algorithm for Poisson model |
|
#--------------------------------# |
|
alpha <- 0.10 |
|
beta <- 0.10 |
|
eps <- 0.01 # epsilon |
|
v <- 0.01 # v_0 |
|
lambdaA <- 12 |
|
c <- qpois(alpha, lambda = v*10, lower.tail = FALSE) |
|
power<- 1 - ppois(c, lambda = v*lambdaA) |
|
while (power < 1 - beta) { |
|
v <- v + eps |
|
c <- qpois(alpha, lambda = v*10, lower.tail = FALSE) |
|
power <- 1 - ppois(c, lambda = v*lambdaA) |
|
} |
|
round(v, 2) # sample volume |
|
c # compliance threshold |
|
round(power, 2) # detection power for lambdaA=12 |
|
|
|
# optional commands |
|
# detection power for lambdaA= 11.5, 12.5 and 13 |
|
lambdaA <- 11.5 |
|
round(1 - ppois(c, lambda = v*lambdaA), 2) |
|
lambdaA <- 12.5 |
|
round(1 - ppois(c, lambda = v*lambdaA), 2) |
|
lambdaA <- 13 |
|
round(1 - ppois(c, lambda = v*lambdaA), 2) |
|
|
|
#------------------------------------------# |
|
# Algorithm for stratified Poisson model |
|
#------------------------------------------# |
|
H <- 4 # number of strata |
|
e <- 1 # maximum estimation error |
|
Vh <- c(135, 75, 40, 20) # volume by stratum |
|
alpha.h <- c(0.02, 0.01, 0.01, 0.01)# alpha by stratum |
|
sum(alpha.h) # alpha |
|
a.h <- c(1, 1, 1, 1) # minimum value for the concentration by stratum |
|
b.h <- c(25, 40, 30, 60) # maximum value for the concentration by stratum |
|
V <- sum(Vh) # total volume |
|
Wh <- Vh/V # weights |
|
Wh |
|
eh <- e/(Wh*H) # maximum estimation error by stratum |
|
eh |
|
|
|
# computation of the sample sizes n_1, n_2, n_3 and n_4 |
|
source("https://raw.githubusercontent.com/eliardocosta/samplesize/master/sample.sizeP.R") |
|
sample.sizeP(d = alpha.h[1], a = a.h[1], b = b.h[1], ea = eh[1]) |
|
sample.sizeP(d = alpha.h[2], a = a.h[2], b = b.h[2], ea = eh[2]) |
|
sample.sizeP(d = alpha.h[3], a = a.h[3], b = b.h[3], ea = eh[3]) |
|
sample.sizeP(d = alpha.h[4], a = a.h[4], b = b.h[4], ea = eh[4]) |
|
|
|
#-------------------------------------------# |
|
# Algorithm for the Poisson process model |
|
#-------------------------------------------# |
|
# Code for fitting model |
|
# $\lambda(t|\delta,\gamma)=\delta\gamma t^{\gamma-1}$. |
|
V <- 270 # total volume |
|
t <- seq(5, V-5, length = 30) # deballasting volumes |
|
X <- c(0, 1, 6, 4, 8, 9, 9, 13, 12, 17, 18, 19, |
|
14, 21, 30, 30, 25, 28, 24, 31, 31, 24, 45, |
|
38, 31, 36, 40, 46, 53, 45) # organism counts |
|
lambda <- function(d, g, t) { # concentration function |
|
d*g*t^(g - 1) |
|
} |
|
d <- 0.05 # delta parameter |
|
g <- 2.1 # gamma parameter |
|
round(d*V^(g - 1), 3) # true concentration |
|
|
|
# (-1)*log-likelihood function |
|
logver <- function(d, g) { |
|
d*g*sum(t^(g - 1)) - log(d*g)*sum(X) - (g - 1)*sum(X*log(t)) + sum(log(factorial(X))) |
|
} |
|
if (is.element("stats4", installed.packages()[,1]) == FALSE) { |
|
install.packages("stats4") |
|
library(stats4) |
|
} else { |
|
library(stats4) |
|
} |
|
fit <- mle(logver, start = list(d = 0.1, g = 0.1)) # fitting the model |
|
summary(fit, digits = 2) |
|
d.est <- coef(fit)[1] # delta estimate |
|
d.est |
|
g.est <- coef(fit)[2] # gamma estimate |
|
g.est |
|
round(sqrt(vcov(fit)), 2) # standard errors |
|
round(confint(fit), 2) # confidence intervals (95%) |
|
|
|
# concentration estimate |
|
lambda.est <- d.est*V^(g.est - 1) |
|
round(lambda.est, 2) |
|
|
|
# variance for lambda estimator via the delta method |
|
var.lambda <- function(d, g, fit) { |
|
t(c(V^(g - 1),d*log(V)*V^(g - 1)))%*%vcov(fit)%*%c(V^(g - 1),d*log(V)*V^(g - 1)) |
|
} |
|
|
|
# std error for lambda estimate |
|
se.lambda <- sqrt(var.lambda(d.est, g.est, fit)) |
|
round(se.lambda, 2) |
|
|
|
# 95% confidence interval for lambda |
|
round(lambda.est + c(-1, 1)*qnorm(1 - 0.025)*se.lambda, 2) |
|
|
|
# Z statistics |
|
round((lambda.est - 10)/se.lambda, 2) |
|
|
|
# Simulation of the data depicted in Table 1. |
|
X <- numeric() |
|
con <- numeric() |
|
j <- 1 |
|
for (i in t) { |
|
con[j] <- lambda(d, g, i) |
|
set.seed(2014 + j) |
|
X[j] <- rpois(1, con[j]) |
|
j <- j + 1 |
|
} |
|
X # simulated data |
|
|
|
# Code for fitting model |
|
# $\lambda(t|\delta,\gamma,\eta)=\eta+[(t-\gamma)/\delta]^2$ |
|
V <- 270 # total volume |
|
t <- seq(5, V - 5, length = 30) # deballasting volumes |
|
X <- c(14, 16, 26, 16, 20, 18, 15, 18, 13, 17, 15, 14, |
|
9, 12, 18, 16, 12, 13, 9, 14, 13, 8, 22, 17, 11, 15, |
|
17, 22, 27, 21) # organism counts |
|
lambda <- function(d, e, g, t) { # concentration function |
|
e + ((t - g)/d)^2 |
|
} |
|
d <- 40 # delta parameter |
|
e <- 12.5 # eta parameter |
|
g <- 130 # gamma parameter |
|
round((e*V + (V - g)^3/(3*d^2) + g^3/(3*d^2))/V, 2) # true concentration |
|
|
|
# (-1)*log-likelihood function |
|
logver <- function(d, e, g) { |
|
length(X)*e + sum(((t - g)/d)^2) - sum(X*log(e + ((t - g)/d)^2)) + sum(log(factorial(X))) |
|
} |
|
fit <- mle(logver, start = list(d = 76, e = min(X), g=t[which.min(X)])) # fitting the model |
|
summary(fit, digits = 3) |
|
d.est <- coef(fit)[1] # delta estimate |
|
d.est |
|
e.est <- coef(fit)[2] # eta estimate |
|
e.est |
|
g.est <- coef(fit)[3] # gamma estimate |
|
g.est |
|
round(sqrt(vcov(fit)), 2) # standard errors |
|
round(confint(fit), 2) # confidence intervals (95%) |
|
|
|
# concentration estimate |
|
lambda.est <- (e.est*V + (V - g.est)^3/(3*d.est^2) + g.est^3/(3*d.est^2))/V |
|
round(lambda.est, 2) |
|
|
|
# variance for lambda estimator via delta method |
|
var.lambda <- function(d, e, g, fit) { |
|
t(c(-2*((V - g)^3 + g^3)/(3*V*d^3), ((V - g)^3 + g^3)/(3*V*d^2), |
|
(g^2 - (V - g)^2)/(V*d^2)))%*%vcov(fit)%*%c(-2*((V - g)^3 + g^3)/(3*V*d^3), |
|
((V - g)^3 + g^3)/(3*V*d^2), (g^2 - (V - g)^2)/(V*d^2)) |
|
} |
|
|
|
# std error for lambda estimate |
|
se.lambda <- sqrt(var.lambda(d.est, e.est, g.est, fit)) |
|
round(se.lambda, 2) |
|
|
|
# 95% confidence interval for lambda |
|
round(lambda.est + c(-1, 1)*qnorm(1 - 0.025)*se.lambda, 2) |
|
|
|
# Z statistics |
|
round((lambda.est - 10)/se.lambda, 2) |
|
|
|
# Simulation of the data depicted in Table 2. |
|
X <- numeric() |
|
con <- numeric() |
|
j <- 1 |
|
for (i in t) { |
|
con[j] <- lambda(d, e, g, i) |
|
set.seed(2014 + j) |
|
X[j] <- rpois(1, con[j]) |
|
j <- j + 1 |
|
} |
|
X # simulated data |
|
|
|
#---------------------------------------------# |
|
# Algorithm for the negative binomial model |
|
#---------------------------------------------# |
|
alpha <- 0.10 |
|
beta <- 0.10 |
|
phi <- 10 |
|
w <- 0.001 |
|
lambdaA <- 12 |
|
n <- 1 |
|
c <- qnbinom(alpha, mu = n*w*10, size = n*phi, lower.tail = FALSE) |
|
power <- 1 - pnbinom(c, mu = n*w*lambdaA, size = n*phi) |
|
while (power < 1 - beta) { |
|
n <- n + 1 |
|
c <- qnbinom(alpha, mu = n*w*10, size = n*phi, lower.tail = FALSE) |
|
power <- 1 - pnbinom(c, mu = n*w*lambdaA, size = n*phi) |
|
} |
|
round(n*w, 2) # sample volume |
|
c # compliance threshold |
|
round(power, 2) # detection power for lambdaA=12 |
|
|
|
# optional commands |
|
# detection powers for lambdaA= 11.5, 12.5 and 13 |
|
lambdaA <- 11.5 |
|
round(1 - pnbinom(c, mu = n*w*lambdaA, size = n*phi), 2) |
|
lambdaA <- 12.5 |
|
round(1 - pnbinom(c, mu = n*w*lambdaA, size = n*phi), 2) |
|
lambdaA <- 13 |
|
round(1 - pnbinom(c, mu = n*w*lambdaA, size = n*phi), 2) |
|
|
|
#---------------------------------------------------------# |
|
# R codes for the analysis of the Gollasch & David data |
|
#---------------------------------------------------------# |
|
# Poisson model |
|
alpha <- 0.05 |
|
beta <- 0.10 |
|
w <- 0.27 # volume sampled by Gollasch&David |
|
n <- 9 # number of replicates |
|
lambdaA <- 12 |
|
c <- qpois(alpha, lambda = n*w*10, lower.tail = FALSE) |
|
c |
|
power <- 1 - ppois(c, lambda = n*w*lambdaA) |
|
round(power, 2) |
|
while (power < 1 - beta) { |
|
n <- n + 1 |
|
c <- qpois(alpha, lambda = n*w*10, lower.tail = FALSE) |
|
power <- 1 - ppois(c, lambda = n*w*lambdaA) |
|
} |
|
n # number of replicates |
|
c # compliance threshold |
|
round(n*w, 2) # sample volume |
|
round(power, 2) # detection power for lambdaA=12 |
|
|
|
# optional commands |
|
# detection powers for lambdaA= 11.5, 12.5 and 13 |
|
lambdaA <- 11.5 |
|
round(1 - ppois(c, lambda = n*w*lambdaA), 2) |
|
lambdaA <- 12.5 |
|
round(1 - ppois(c, lambda = n*w*lambdaA), 2) |
|
lambdaA <- 13 |
|
round(1 - ppois(c, lambda = n*w*lambdaA), 2) |
|
|
|
# estimating the parameter phi of the negative binomial model |
|
if (is.element("MASS", installed.packages()[,1])==FALSE) { |
|
install.packages("MASS") |
|
library(MASS) |
|
} else { |
|
library(MASS) |
|
} |
|
case1 <- c(4, 2, 1, 2, 2, 2, 6, 5, 3) # uptake |
|
mod1 <- fitdistr(case1, "negative binomial") |
|
mod1 |
|
confint(mod1) |
|
case2 <- c(8, 6, 2, 6, 3, 4, 29, 17, 25) # discharge |
|
mod2 <- fitdistr(case2, "negative binomial") |
|
mod2 |
|
confint(mod2) |
|
|
|
# Negative binomial model |
|
alpha <- 0.05 |
|
beta <- 0.10 |
|
phi <- 1.66 # estimate of phi |
|
w <- 0.27 # volume sampled by Gollasch & David |
|
n <- 9 # number of replicates |
|
lambdaA <- 12 |
|
c <- qnbinom(alpha, mu = n*w*10, size = n*phi, lower.tail = FALSE) |
|
c |
|
power <- 1 - pnbinom(c, mu = n*w*lambdaA, size = n*phi) |
|
round(power, 2) |
|
while (power < 1 - beta) { |
|
n <- n + 1 |
|
c <- qnbinom(alpha, mu = n*w*10, size = n*phi, lower.tail = FALSE) |
|
power <- 1 - pnbinom(c, mu = n*w*lambdaA, size = n*phi) |
|
} |
|
n # number of replicates |
|
c # compliance threshold |
|
round(n*w, 2) # sample volume |
|
round(power, 2) # detection power for lambdaA=12 |
|
|
|
# optional commands |
|
# detection powers for lambdaA= 11.5, 12.5 and 13 |
|
lambdaA <- 11.5 |
|
round(1 - pnbinom(c, mu = n*w*lambdaA, size = n*phi), 2) |
|
lambdaA <- 12.5 |
|
round(1 - pnbinom(c, mu = n*w*lambdaA, size = n*phi), 2) |
|
lambdaA <- 13 |
|
round(1 - pnbinom(c, mu = n*w*lambdaA, size = n*phi), 2) |