The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
library(anovir)The code below repeats analyses found in the supplementary files of Agnew 2019 [1].
Data from [2]
data01 <- subset(data_blanford,
(data_blanford$block == 3) &
((data_blanford$treatment == 'cont') | (data_blanford$treatment == 'Bb06')) &
(data_blanford$day > 0)
)
head(data01, 3)
#> block treatment replicate_cage day censor d inf t fq
#> 450 3 cont 1 1 0 1 0 1 0
#> 451 3 cont 1 2 0 1 0 2 0
#> 452 3 cont 1 3 0 1 0 3 0
m01_prep_function <- function(a1, b1, a2, b2){
nll_basic(a1, b1, a2, b2,
data = data01,
time = t,
censor = censor,
infected_treatment = inf,
d1 = 'Weibull', d2 = 'Weibull')
}
# starting values taken from linear regression of complementary
# log-log transformed cumulative survival data
m01 <- mle2(m01_prep_function,
start = list(a1 = 3.343, b1 = 0.792, a2 = 2.508, b2 = 0.493)
)
summary(m01)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m01_prep_function, start = list(a1 = 3.343,
#> b1 = 0.792, a2 = 2.508, b2 = 0.493))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 2.845028 0.069054 41.2000 < 2.2e-16 ***
#> b1 0.482718 0.043025 11.2195 < 2.2e-16 ***
#> a2 2.580839 0.028630 90.1457 < 2.2e-16 ***
#> b2 0.183062 0.031634 5.7868 7.174e-09 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1293.144
confint(m01)
#> 2.5 % 97.5 %
#> a1 2.723105 2.9989058
#> b1 0.407109 0.5787759
#> a2 2.524759 2.6396903
#> b2 0.121188 0.2507871
# log-likelihood based on estimates of linear regression
m02<- mle2(m01_prep_function,
start = list(a1 = 3.343, b1 = 0.792, a2 = 2.508, b2 = 0.493),
eval.only = TRUE
)
summary(m02)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m01_prep_function, start = list(a1 = 3.343,
#> b1 = 0.792, a2 = 2.508, b2 = 0.493), eval.only = TRUE)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 3.343 NA NA NA
#> b1 0.792 NA NA NA
#> a2 2.508 NA NA NA
#> b2 0.493 NA NA NA
#>
#> -2 log L: 1381.426
AICc(m01, m02, nobs = sum(data01$fq))
#> AICc df
#> 1 1301.286 4
#> 2 1389.568 4Data from [3,4]
data01 <- data_lorenz
head(data01, 3)
#> Infectious.dose Food Sex Spore.Count t censored d g
#> 1 10000 100 F 425000 13.0 0 1 1
#> 2 10000 50 F 22000 3.5 0 1 1
#> 3 10000 100 F 465000 20.5 0 1 1
### Model 1
m01_prep_function <- function(a1, b1, a2, b2){
nll_basic(a1, b1, a2, b2,
data = data01,
time = t,
censor = censored,
infected_treatment = g,
d1 = 'Gumbel', d2 = 'Weibull')
}
m01 <- mle2(m01_prep_function,
start = list(a1 = 20, b1 = 5, a2 = 3, b2 = 1)
)
summary(m01)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m01_prep_function, start = list(a1 = 20, b1 = 5,
#> a2 = 3, b2 = 1))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 23.215721 0.665292 34.8955 < 2.2e-16 ***
#> b1 4.674867 0.406836 11.4908 < 2.2e-16 ***
#> a2 3.019835 0.028913 104.4438 < 2.2e-16 ***
#> b2 0.210731 0.023552 8.9475 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1528.421
### Model 2
# copy 'nll_basic' & make each parameter function a function of food treatment
nll_basic2 <- nll_basic
body(nll_basic2)[[2]] <- substitute(
pfa1 <- a1 + ifelse(data01$Food == 50, a1i, -a1i)
)
body(nll_basic2)[[3]] <- substitute(
pfb1 <- b1 + ifelse(data01$Food == 50, b1i, -b1i)
)
body(nll_basic2)[[4]] <- substitute(
pfa2 <- a2 + ifelse(data01$Food == 50, a2i, -a2i)
)
body(nll_basic2)[[5]] <- substitute(
pfb2 <- b2 + ifelse(data01$Food == 50, b2i, -b2i)
)
# update formals
formals(nll_basic2) <- alist(
a1 = a1, a1i = a1i,
b1 = b1, b1i = b1i,
a2 = a2, a2i = a2i,
b2 = b2, b2i = b2i,
data = data01,
time = t,
censor = censored,
infected_treatment = g,
d1 = 'Gumbel', d2 = 'Weibull')
m02_prep_function <- function(a1, a1i, b1, b1i, a2, a2i, b2, b2i){
nll_basic2(a1, a1i, b1, b1i, a2, a2i, b2, b2i)
}
m02 <- mle2(m02_prep_function,
start = list(
a1 = 23.2, a1i = 0,
b1 = 4.6, b1i = 0,
a2 = 3.0, a2i = 0,
b2 = 0.2, b2i = 0)
)
summary(m02)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m02_prep_function, start = list(a1 = 23.2, a1i = 0,
#> b1 = 4.6, b1i = 0, a2 = 3, a2i = 0, b2 = 0.2, b2i = 0))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 23.218273 0.673537 34.4721 <2e-16 ***
#> a1i -0.336986 0.673538 -0.5003 0.6168
#> b1 4.686262 0.409784 11.4359 <2e-16 ***
#> b1i -0.312469 0.409784 -0.7625 0.4457
#> a2 3.014350 0.028045 107.4842 <2e-16 ***
#> a2i -0.042484 0.028045 -1.5149 0.1298
#> b2 0.202349 0.023291 8.6880 <2e-16 ***
#> b2i -0.016894 0.023291 -0.7253 0.4682
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1521.028
### Model 3
# create columns in data frame
# with dummary variables for dose treatments
data01$d5 <- ifelse(data01$Infectious.dose == 5000, 1, 0)
data01$d10 <- ifelse(data01$Infectious.dose == 10000, 1, 0)
data01$d20 <- ifelse(data01$Infectious.dose == 20000, 1, 0)
data01$d40 <- ifelse(data01$Infectious.dose == 40000, 1, 0)
data01$d80 <- ifelse(data01$Infectious.dose == 80000, 1, 0)
data01$d160 <- ifelse(data01$Infectious.dose == 160000, 1, 0)
head(data01)
#> Infectious.dose Food Sex Spore.Count t censored d g d5 d10 d20 d40 d80
#> 1 10000 100 F 425000 13.0 0 1 1 0 1 0 0 0
#> 2 10000 50 F 22000 3.5 0 1 1 0 1 0 0 0
#> 3 10000 100 F 465000 20.5 0 1 1 0 1 0 0 0
#> 8 0 50 F NA 17.0 0 1 0 0 0 0 0 0
#> 10 160000 50 F 295000 17.5 0 1 1 0 0 0 0 0
#> 11 160000 100 F 0 4.0 0 1 1 0 0 0 0 0
#> d160
#> 1 0
#> 2 0
#> 3 0
#> 8 0
#> 10 1
#> 11 1
# make parameter functions 'pfa2' and 'pfb2' functions of dose
# using columns in data01
nll_basic3 <- nll_basic
body(nll_basic3)[[4]] <- substitute(
pfa2 <- a2 + a5 * data01$d5
+ a10 * data01$d10
+ a20 * data01$d20
+ a40 * data01$d40
+ a80 * data01$d80
- (a5 + a10 + a20 + a40 + a80) * data01$d160
)
body(nll_basic3)[[5]] <- substitute(
pfb2 <- b2 + b5 * data01$d5
+ b10 * data01$d10
+ b20 * data01$d20
+ b40 * data01$d40
+ b80 * data01$d80
- (b5 + b10 + b20 + b40 + b80) * data01$d160
)
# update formals
formals(nll_basic3) <- alist(
a1 = a1, b1 = b1,
a2 = a2, a5 = a5, a10 = a10, a20 = a20, a40 = a40, a80 = a80,
b2 = b2, b5 = b5, b10 = b10, b20 = b20, b40 = b40, b80 = b80,
data = data01,
time = t, censor = censored, infected_treatment = g,
d1 = 'Gumbel', d2 = 'Weibull'
)
m03_prep_function <- function(
a1, b1,
a2, a5, a10, a20, a40, a80,
b2, b5, b10, b20, b40, b80){
nll_basic3(
a1, b1,
a2, a5, a10, a20, a40, a80,
b2, b5, b10, b20, b40, b80)
}
m03 <- mle2(m03_prep_function,
start = list(
a1 = 23, b1 = 4.6,
a2 = 3, a5 = 0, a10 = 0, a20 = 0, a40 = 0, a80 = 0,
b2 = 0.2, b5 = 0, b10 = 0, b20 = 0, b40 = 0, b80 = 0)
)
summary(m03)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m03_prep_function, start = list(a1 = 23, b1 = 4.6,
#> a2 = 3, a5 = 0, a10 = 0, a20 = 0, a40 = 0, a80 = 0, b2 = 0.2,
#> b5 = 0, b10 = 0, b20 = 0, b40 = 0, b80 = 0))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 23.2051619 0.6525030 35.5633 < 2.2e-16 ***
#> b1 4.5880595 0.4189079 10.9524 < 2.2e-16 ***
#> a2 3.0078037 0.0283159 106.2232 < 2.2e-16 ***
#> a5 0.1840405 0.0590506 3.1167 0.001829 **
#> a10 0.0292572 0.0674387 0.4338 0.664409
#> a20 0.0397519 0.0458397 0.8672 0.385836
#> a40 -0.0503364 0.0465339 -1.0817 0.279379
#> a80 -0.0857037 0.0423325 -2.0245 0.042915 *
#> b2 0.1941075 0.0255850 7.5868 3.28e-14 ***
#> b5 -0.0332491 0.0411073 -0.8088 0.418609
#> b10 0.0928078 0.0797457 1.1638 0.244506
#> b20 -0.0136344 0.0416699 -0.3272 0.743516
#> b40 0.0068783 0.0374779 0.1835 0.854381
#> b80 -0.0221821 0.0349606 -0.6345 0.525762
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1501.587
### Model 4
# this model estimates the full dose * food interaction
# and replaces the approximate approach used in the original text
data01 <- data_lorenz
head(data01)
#> Infectious.dose Food Sex Spore.Count t censored d g
#> 1 10000 100 F 425000 13.0 0 1 1
#> 2 10000 50 F 22000 3.5 0 1 1
#> 3 10000 100 F 465000 20.5 0 1 1
#> 8 0 50 F NA 17.0 0 1 0
#> 10 160000 50 F 295000 17.5 0 1 1
#> 11 160000 100 F 0 4.0 0 1 1
# create columns with dummy variables
data01$d5 <- ifelse(data01$Infectious.dose == 5000, 1, 0)
data01$d10 <- ifelse(data01$Infectious.dose == 10000, 1, 0)
data01$d20 <- ifelse(data01$Infectious.dose == 20000, 1, 0)
data01$d40 <- ifelse(data01$Infectious.dose == 40000, 1, 0)
data01$d80 <- ifelse(data01$Infectious.dose == 80000, 1, 0)
data01$d160 <- ifelse(data01$Infectious.dose == 160000, 1, 0)
data01$af50 <- ifelse(data01$Food == 50, 1, 0)
data01$af100 <- ifelse(data01$Food == 100, 1, 0)
data01$f50d5 <- ifelse(data01$Infectious.dose == 5000 & data01$Food == 50, 1, 0)
data01$f50d10 <- ifelse(data01$Infectious.dose == 10000 & data01$Food == 50, 1, 0)
data01$f50d20 <- ifelse(data01$Infectious.dose == 20000 & data01$Food == 50, 1, 0)
data01$f50d40 <- ifelse(data01$Infectious.dose == 40000 & data01$Food == 50, 1, 0)
data01$f50d80 <- ifelse(data01$Infectious.dose == 80000 & data01$Food == 50, 1, 0)
data01$f50d160 <- ifelse(data01$Infectious.dose == 160000 & data01$Food == 50, 1, 0)
data01$f100d5 <- ifelse(data01$Infectious.dose == 5000 & data01$Food == 100, 1, 0)
data01$f100d10 <- ifelse(data01$Infectious.dose == 10000 & data01$Food == 100, 1, 0)
data01$f100d20 <- ifelse(data01$Infectious.dose == 20000 & data01$Food == 100, 1, 0)
data01$f100d40 <- ifelse(data01$Infectious.dose == 40000 & data01$Food == 100, 1, 0)
data01$f100d80 <- ifelse(data01$Infectious.dose == 80000 & data01$Food == 100, 1, 0)
data01$f100d160 <- ifelse(data01$Infectious.dose == 160000 & data01$Food == 100, 1, 0)
head(data01)
#> Infectious.dose Food Sex Spore.Count t censored d g d5 d10 d20 d40 d80
#> 1 10000 100 F 425000 13.0 0 1 1 0 1 0 0 0
#> 2 10000 50 F 22000 3.5 0 1 1 0 1 0 0 0
#> 3 10000 100 F 465000 20.5 0 1 1 0 1 0 0 0
#> 8 0 50 F NA 17.0 0 1 0 0 0 0 0 0
#> 10 160000 50 F 295000 17.5 0 1 1 0 0 0 0 0
#> 11 160000 100 F 0 4.0 0 1 1 0 0 0 0 0
#> d160 af50 af100 f50d5 f50d10 f50d20 f50d40 f50d80 f50d160 f100d5 f100d10
#> 1 0 0 1 0 0 0 0 0 0 0 1
#> 2 0 1 0 0 1 0 0 0 0 0 0
#> 3 0 0 1 0 0 0 0 0 0 0 1
#> 8 0 1 0 0 0 0 0 0 0 0 0
#> 10 1 1 0 0 0 0 0 0 1 0 0
#> 11 1 0 1 0 0 0 0 0 0 0 0
#> f100d20 f100d40 f100d80 f100d160
#> 1 0 0 0 0
#> 2 0 0 0 0
#> 3 0 0 0 0
#> 8 0 0 0 0
#> 10 0 0 0 0
#> 11 0 0 0 1
nll_basic4 <- nll_basic
# make 'pfa2' and 'pfb2' functions of food-by-dose interaction
body(nll_basic4)[[4]] <- substitute(
pfa2 <- a2 + a5 * data01$d5
+ a10 * data01$d10
+ a20 * data01$d20
+ a40 * data01$d40
+ a80 * data01$d80
- (a5 + a10 + a20 + a40 + a80) * data01$d160
+ af * data01$af50
- af * data01$af100
+ afd5 * data01$f50d5
+ afd10 * data01$f50d10
+ afd20 * data01$f50d20
+ afd40 * data01$f50d40
+ afd80 * data01$f50d80
- (afd5 + afd10 + afd20 + afd40 + afd80) * data01$f50d160
- afd5 * data01$f100d5
- afd10 * data01$f100d10
- afd20 * data01$f100d20
- afd40 * data01$f100d40
- afd80 * data01$f100d80
+ (afd5 + afd10 + afd20 + afd40 + afd80) * data01$f100d160
)
body(nll_basic4)[[5]] <- substitute(
pfb2 <- b2 + b5 * data01$d5
+ b10 * data01$d10
+ b20 * data01$d20
+ b40 * data01$d40
+ b80 * data01$d80
- (b5 + b10 + b20 + b40 + b80) * data01$d160
+ bf * data01$af50
- bf * data01$af100
+ bfd5 * data01$f50d5
+ bfd10 * data01$f50d10
+ bfd20 * data01$f50d20
+ bfd40 * data01$f50d40
+ bfd80 * data01$f50d80
- (bfd5 + bfd10 + bfd20 + bfd40 + bfd80) * data01$f50d160
- bfd5 * data01$f100d5
- bfd10 * data01$f100d10
- bfd20 * data01$f100d20
- bfd40 * data01$f100d40
- bfd80 * data01$f100d80
+ (bfd5 + bfd10 + bfd20 + bfd40 + bfd80) * data01$f100d160
)
formals(nll_basic4) <- alist(
a1 = a1, b1 = b1,
a2 = a2, a5 = a5, a10 = a10, a20 = a20, a40 = a40, a80 = a80,
af = af, afd5 = afd5, afd10 = afd10, afd20 = afd20, afd40 = afd40, afd80 = afd80,
b2 = b2, b5 = b5, b10 = b10, b20 = b20, b40 = b40, b80 = b80,
bf = bf, bfd5 = bfd5, bfd10 = bfd10, bfd20 = bfd20, bfd40 = bfd40, bfd80 = bfd80,
data = data01,
time = t,
censor = censored,
infected_treatment = g,
d1 = 'Gumbel', d2 = 'Weibull')
m04_prep_function <- function(
a1 = a1, b1 = b1,
a2 = a2, a5 = a5, a10 = a10, a20 = a20, a40 = a40, a80 = a80,
af = af, afd5 = afd5, afd10 = afd10, afd20 = afd20, afd40 = afd40, afd80 = afd80,
b2 = b2, b5 = b5, b10 = b10, b20 = b20, b40 = b40, b80 = b80,
bf = bf, bfd5 = bfd5, bfd10 = bfd10, bfd20 = bfd20, bfd40 = bfd40, bfd80 = bfd80
){nll_basic4(
a1 = a1, b1 = b1,
a2 = a2, a5 = a5, a10 = a10, a20 = a20, a40 = a40, a80 = a80,
af = af, afd5 = afd5, afd10 = afd10, afd20 = afd20, afd40 = afd40, afd80 = afd80,
b2 = b2, b5 = b5, b10 = b10, b20 = b20, b40 = b40, b80 = b80,
bf = bf, bfd5 = bfd5, bfd10 = bfd10, bfd20 = bfd20, bfd40 = bfd40, bfd80 = bfd80
)}
m04 <- mle2(m04_prep_function,
start = list(
a1 = 23, b1 = 4.6,
a2 = 3, a5 = 0.18, a10 = 0.03, a20 = 0.04, a40 = -0.05, a80 = -0.08,
af = 0, afd5 = 0, afd10 = 0, afd20 = 0, afd40 = 0, afd80 = 0,
b2 = 0.2, b5 = -0.03, b10 = 0.09, b20 = -0.01, b40 = 0.01, b80 = -0.02,
bf = 0, bfd5 = 0, bfd10 = 0, bfd20 = 0, bfd40 = 0, bfd80 = 0)
)
coef(m04)
#> a1 b1 a2 a5 a10 a20
#> 23.159088319 4.589283055 2.993905503 0.143349586 0.034624488 0.054680576
#> a40 a80 af afd5 afd10 afd20
#> -0.033525078 -0.092740294 -0.039210763 -0.078448071 -0.043172902 0.019762051
#> afd40 afd80 b2 b5 b10 b20
#> 0.062786054 0.032678441 0.169070068 -0.107516377 0.065835289 0.004724730
#> b40 b80 bf bfd5 bfd10 bfd20
#> 0.031260873 0.010136435 0.003320434 0.050601433 0.046724363 -0.035313347
#> bfd40 bfd80
#> -0.002633637 -0.095967757
### Model 5
# recall 'nll_basic3' from Model 3 above
# return 'pfb2' to original form
body(nll_basic3)
#> {
#> pfa1 <- a1
#> pfb1 <- b1
#> pfa2 <- a2 + a5 * data01$d5 + a10 * data01$d10 + a20 * data01$d20 +
#> a40 * data01$d40 + a80 * data01$d80 - (a5 + a10 + a20 +
#> a40 + a80) * data01$d160
#> pfb2 <- b2 + b5 * data01$d5 + b10 * data01$d10 + b20 * data01$d20 +
#> b40 * data01$d40 + b80 * data01$d80 - (b5 + b10 + b20 +
#> b40 + b80) * data01$d160
#> t <- data[[deparse(substitute(time))]]
#> g <- data[[deparse(substitute(infected_treatment))]]
#> d <- data[[deparse(substitute(censor))]] * -1 + 1
#> z1 <- P_get_zx(t, pfa1, pfb1, d1)
#> z2 <- P_get_zx(t, pfa2, pfb2, d2)
#> h1 <- P_get_hx(t, z1, pfb1, d1)
#> h2 <- P_get_hx(t, z2, pfb2, d2)
#> S1 <- P_get_Sx(t, z1, d1)
#> S2 <- P_get_Sx(t, z2, d2)
#> logl <- 0
#> model <- (d * log(h1 + g * h2) + log(S1) + g * log(S2))
#> if ("fq" %in% colnames(data)) {
#> logl <- sum(data$fq * model)
#> }
#> else {
#> logl <- sum(model)
#> }
#> return(-logl)
#> }
body(nll_basic3)[[5]] <- substitute(pfb2 <- b2)
formals(nll_basic3) <- alist(
a1 = a1, b1 = b1,
a2 = a2, a5 = a5, a10 = a10, a20 = a20, a40 = a40, a80 = a80,
b2 = b2,
data = data01,
time = t, censor = censored, infected_treatment = g,
d1 = 'Gumbel', d2 = 'Weibull'
)
m05_prep_function <- function(a1, b1, a2, a5, a10, a20, a40, a80, b2){
nll_basic3(a1, b1, a2, a5, a10, a20, a40, a80, b2)
}
m05 <- mle2(m05_prep_function,
start = list(
a1 = 23, b1 = 4.6,
a2 = 3, a5 = 0.18, a10 = 0.03, a20 = 0.04, a40 = -0.05, a80 = -0.08,
b2 = 0.2)
)
summary(m05)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m05_prep_function, start = list(a1 = 23, b1 = 4.6,
#> a2 = 3, a5 = 0.18, a10 = 0.03, a20 = 0.04, a40 = -0.05, a80 = -0.08,
#> b2 = 0.2))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 23.165483 0.652940 35.4788 < 2.2e-16 ***
#> b1 4.660075 0.384771 12.1113 < 2.2e-16 ***
#> a2 3.011271 0.027298 110.3099 < 2.2e-16 ***
#> a5 0.194128 0.070401 2.7575 0.005825 **
#> a10 0.038874 0.048343 0.8041 0.421324
#> a20 0.037011 0.046981 0.7878 0.430830
#> a40 -0.049675 0.044162 -1.1248 0.260656
#> a80 -0.092924 0.043514 -2.1355 0.032719 *
#> b2 0.188361 0.020122 9.3609 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1504.575
### Model 6
# make 'pfa2' a linear function of log(Infectious dose)
nll_basic6 <- nll_basic
body(nll_basic6)[[4]] <- substitute(pfa2 <- a2i + a2ii * log(data01$Infectious.dose))
formals(nll_basic6) <- alist(
a1 = a1, b1 = b1, a2i = a2i, a2ii = a2ii, b2 = b2,
data = data01, time = t, censor = censored, infected_treatment = g,
d1 = 'Gumbel', d2 = 'Weibull')
m06_prep_function <- function(a1, b1, a2i, a2ii, b2){
nll_basic6(a1, b1, a2i, a2ii, b2)
}
m06 <- mle2(m06_prep_function,
start = list(a1 = 23, b1 = 4.6, a2i = 4, a2ii = -0.1, b2 = 0.2)
)
summary(m06)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m06_prep_function, start = list(a1 = 23, b1 = 4.6,
#> a2i = 4, a2ii = -0.1, b2 = 0.2))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 23.193699 0.656664 35.3205 < 2.2e-16 ***
#> b1 4.718324 0.383030 12.3184 < 2.2e-16 ***
#> a2i 3.826938 0.189127 20.2348 < 2.2e-16 ***
#> a2ii -0.079691 0.017512 -4.5507 5.346e-06 ***
#> b2 0.185438 0.019378 9.5696 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1506.19
AICc(m01, m02, m03, m04, m05, m06, nobs = 256)
#> AICc df
#> 1 1536.580 4
#> 2 1537.611 8
#> 3 1531.329 14
#> 4 1531.097 26
#> 5 1523.306 9
#> 6 1516.430 5
### Model 6 with different distributions
m06b_prep_function <- function(a1, b1, a2i, a2ii, b2){
nll_basic6(a1, b1, a2i, a2ii, b2, d1 = 'Weibull', d2 = 'Gumbel')
}
m06b <- mle2(m06b_prep_function,
start = list(a1 = 2, b1 = 0.5, a2i = 23, a2ii = -0.1, b2 = 4)
)
m06c_prep_function <- function(a1, b1, a2i, a2ii, b2){
nll_basic6(a1, b1, a2i, a2ii, b2, d1 = 'Weibull', d2 = 'Weibull')
}
m06c <- mle2(m06c_prep_function,
start = list(a1 = 2, b1 = 0.5, a2i = 3.8, a2ii = -0.1, b2 = 0.2)
)
m06d_prep_function <- function(a1, b1, a2i, a2ii, b2){
nll_basic6(a1, b1, a2i, a2ii, b2, d1 = 'Gumbel', d2 = 'Gumbel')
}
m06d <- mle2(m06d_prep_function,
start = list(a1 = 23, b1 = 5, a2i = 23, a2ii = -0.1, b2 = 5)
)
AICc(m06, m06b, m06c, m06d, nobs = 256)
#> AICc df
#> 1 1516.430 5
#> 2 1534.636 5
#> 3 1541.630 5
#> 4 1517.887 5Data from [2]
# data01 <- data_blanford_bl5
# subset data for 'block = 5', treatments 'cont', 'Ma06', 'Ma07', 'Ma08'
data01 <- data_blanford
data01 <- subset(data01,
(data01$block == 5) & (
(data01$treatment == 'cont') |
(data01$treatment == 'Ma06') |
(data01$treatment == 'Ma07') |
(data01$treatment == 'Ma08') ) &
(data01$day > 0)
)
# create column 'g' as index of infected treatment
data01$g <- data01$inf
head(data01)
#> block treatment replicate_cage day censor d inf t fq g
#> 1026 5 cont 1 1 0 1 0 1 2 0
#> 1027 5 cont 1 2 0 1 0 2 0 0
#> 1028 5 cont 1 3 0 1 0 3 0 0
#> 1029 5 cont 1 4 0 1 0 4 2 0
#> 1030 5 cont 1 5 0 1 0 5 3 0
#> 1031 5 cont 1 6 0 1 0 6 3 0
nll_basic2 <- nll_basic
# make 'pfa2' and 'pfb2' functions of fungal treatment
# NB to avoid problems with log(0), set final ifelse 'false' values to 'exp(0)'
body(nll_basic2)[[4]] <- substitute(pfa2 <-
ifelse(((data01$g == 1) & (data01$treatment == 'Ma06')), a2,
ifelse(((data01$g == 1) & (data01$treatment == 'Ma07')), a3,
ifelse(((data01$g == 1) & (data01$treatment == 'Ma08')), a4,
exp(0)
))))
body(nll_basic2)[[5]] <- substitute(pfb2 <-
ifelse(((data01$g == 1) & (data01$treatment == 'Ma06')), b2,
ifelse(((data01$g == 1) & (data01$treatment == 'Ma07')), b3,
ifelse(((data01$g == 1) & (data01$treatment == 'Ma08')), b4,
exp(0)
))))
formals(nll_basic2) <- alist(
a1 = a1, b1 = b1,
a2 = a2, b2 = b2,
a3 = a3, b3 = b3,
a4 = a4, b4 = b4,
data = data01,
time = t, censor = censor, infected_treatment = g,
d1 = 'Weibull', d2 = 'Fréchet')
m01_prep_function <- function(a1, b1, a2, b2, a3, b3, a4, b4){
nll_basic2(a1, b1, a2, b2, a3, b3, a4, b4)
}
m01 <- mle2(m01_prep_function,
start = list(a1 = 2, b1 = 1, a2 = 2, b2 = 1, a3 = 2, b3 = 1, a4 = 2, b4 = 1)
)
summary(m01)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m01_prep_function, start = list(a1 = 2, b1 = 1,
#> a2 = 2, b2 = 1, a3 = 2, b3 = 1, a4 = 2, b4 = 1))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 3.503561 0.102580 34.154 < 2.2e-16 ***
#> b1 0.686278 0.046712 14.692 < 2.2e-16 ***
#> a2 1.871541 0.036112 51.826 < 2.2e-16 ***
#> b2 0.514390 0.031659 16.248 < 2.2e-16 ***
#> a3 1.637038 0.023854 68.627 < 2.2e-16 ***
#> b3 0.335179 0.018259 18.357 < 2.2e-16 ***
#> a4 2.423527 0.091296 26.546 < 2.2e-16 ***
#> b4 0.916577 0.089208 10.275 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 4611.834
confint(m01)
#> 2.5 % 97.5 %
#> a1 3.3227930 3.7290348
#> b1 0.6023677 0.7869933
#> a2 1.8024006 1.9450014
#> b2 0.4580684 0.5836821
#> a3 1.5908580 1.6846702
#> b3 0.3014103 0.3737186
#> a4 2.2615941 2.6305102
#> b4 0.7642829 1.1232615Data from [5,6]
data01 <- data_parker
head(data01, 3)
#> Genotype SD Fecundity Sporulation Status Time dose censored d t g
#> 1 721 16 46 1 1 6 2 0 1 6 1
#> 2 721 0 56 0 1 7 0 0 1 7 0
#> 3 721 8 68 0 0 18 1 1 0 18 1
### Model 1
nll_basic2 <- nll_basic
body(nll_basic2)[[4]] <- substitute(pfa2 <-
a2 + ifelse(data01$dose == 1, a2i,
ifelse(data01$dose == 2, a2ii,
ifelse(data01$dose == 3, -(a2i + a2ii),
exp(0)
)))
)
body(nll_basic2)[[5]] <- substitute(pfb2 <-
b2 + ifelse(data01$dose == 1, b2i,
ifelse(data01$dose == 2, b2ii,
ifelse(data01$dose == 3, -(b2i + b2ii),
exp(0)
)))
)
formals(nll_basic2) <- alist(
a1 = a1, b1 = b1,
a2 = a2, a2i = a2i, a2ii = a2ii,
b2 = b2, b2i = b2i, b2ii = b2ii,
data = data01,
time = t, censor = censored, infected_treatment = g,
d1 = 'Frechet', d2 = 'Frechet')
m01_prep_function <- function(
a1, b1, a2, a2i, a2ii, b2, b2i, b2ii){
nll_basic2(a1, b1, a2, a2i, a2ii, b2, b2i, b2ii)
}
m01 <- mle2(m01_prep_function,
start = list(
a1 = 2, b1 = 1,
a2 = 2, a2i = 0, a2ii = 0,
b2 = 1, b2i = 0, b2ii = 0)
)
summary(m01)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m01_prep_function, start = list(a1 = 2, b1 = 1,
#> a2 = 2, a2i = 0, a2ii = 0, b2 = 1, b2i = 0, b2ii = 0))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 2.3518841 0.0829811 28.3424 < 2.2e-16 ***
#> b1 0.6585920 0.0566255 11.6307 < 2.2e-16 ***
#> a2 2.0466969 0.0692055 29.5742 < 2.2e-16 ***
#> a2i 0.2784921 0.1043580 2.6686 0.007616 **
#> a2ii -0.0291581 0.0684531 -0.4260 0.670139
#> b2 0.4802869 0.0504207 9.5256 < 2.2e-16 ***
#> b2i 0.1919610 0.0811746 2.3648 0.018040 *
#> b2ii -0.0018653 0.0565409 -0.0330 0.973682
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1557.2
### Model 2
# make 'pfa2' and 'pfb2' linear functions of log(dose)
body(nll_basic2)[[4]] <- substitute(pfa2 <-
a2 + ifelse(data01$dose == 1, a2i,
ifelse(data01$dose == 2, 0,
ifelse(data01$dose == 3, -a2i,
exp(0)
)))
)
body(nll_basic2)[[5]] <- substitute(pfb2 <-
b2 + ifelse(data01$dose == 1, b2i,
ifelse(data01$dose == 2, 0,
ifelse(data01$dose == 3, -b2i,
exp(0)
)))
)
formals(nll_basic2) <- alist(
a1 = a1, b1 = b1,
a2 = a2, a2i = a2i,
b2 = b2, b2i = b2i,
data = data01,
time = t, censor = censored, infected_treatment = g,
d1 = 'Frechet', d2 = 'Frechet')
m02_prep_function <- function(
a1, b1, a2, a2i, b2, b2i){
nll_basic2(a1, b1, a2, a2i, b2, b2i)
}
m02 <- mle2(m02_prep_function,
start = list(
a1 = 2, b1 = 1, a2 = 2, a2i = 0, b2 = 1, b2i = 0)
)
summary(m02)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m02_prep_function, start = list(a1 = 2, b1 = 1,
#> a2 = 2, a2i = 0, b2 = 1, b2i = 0))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 2.355514 0.082213 28.6515 < 2.2e-16 ***
#> b1 0.660935 0.056451 11.7082 < 2.2e-16 ***
#> a2 2.040629 0.064046 31.8621 < 2.2e-16 ***
#> a2i 0.247347 0.062753 3.9416 8.095e-05 ***
#> b2 0.477763 0.047596 10.0379 < 2.2e-16 ***
#> b2i 0.187200 0.049455 3.7853 0.0001535 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1557.45
### Model 3
nll_basic3 <- nll_basic
body(nll_basic3)[[4]] <- substitute(pfa2 <-
a2 + ifelse(data01$g == 1,
ifelse(data01$Sporulation == 1, a2sp, -a2sp),
exp(0))
)
body(nll_basic3)[[5]] <- substitute(pfb2 <-
b2 + ifelse(data01$g == 1,
ifelse(data01$Sporulation == 1, b2sp, -b2sp),
exp(0))
)
formals(nll_basic3) <- alist(
a1 = a1, b1 = b1,
a2 = a2, a2sp = a2sp,
b2 = b2, b2sp = b2sp,
data = data01,
time = t, censor = censored, infected_treatment = g,
d1 = 'Fréchet', d2 = 'Weibull')
m03_prep_function <- function(a1, b1, a2, a2sp, b2, b2sp){
nll_basic3(a1, b1, a2, a2sp, b2, b2sp)
}
m03 <- mle2(m03_prep_function,
start = list(
a1 = 2.35, b1 = 0.66,
a2 = 2, a2sp = 0,
b2 = 0.5, b2sp = 0)
)
summary(m03)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m03_prep_function, start = list(a1 = 2.35, b1 = 0.66,
#> a2 = 2, a2sp = 0, b2 = 0.5, b2sp = 0))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 2.170488 0.066675 32.5531 < 2.2e-16 ***
#> b1 0.609945 0.047171 12.9304 < 2.2e-16 ***
#> a2 2.735482 0.181460 15.0748 < 2.2e-16 ***
#> a2sp -0.652283 0.176444 -3.6968 0.0002183 ***
#> b2 0.328575 0.079697 4.1228 3.743e-05 ***
#> b2sp -0.126966 0.078484 -1.6177 0.1057213
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1513.319
### Model 4
nll_basic4 <- nll_basic
data01 <- data_parker
body(nll_basic4)[[4]] <- substitute(pfa2 <-
a2 + ifelse(data01$g == 1,
ifelse(data01$Sporulation == 1,
ifelse(data01$dose == 1, a2sp + a2d1,
ifelse(data01$dose == 3, a2sp - a2d1, a2sp)),
exp(0)),
exp(0))
)
body(nll_basic4)[[5]] <- substitute(pfb2 <-
b2 + ifelse(data01$g == 1,
ifelse(data01$Sporulation == 1,
ifelse(data01$dose == 1, b2sp + b2d1,
ifelse(data01$dose == 3, b2sp - b2d1, b2sp)),
exp(0)),
exp(0))
)
formals(nll_basic4) <- alist(
a1 = a1, b1 = b1,
a2 = a2, a2sp = a2sp, a2d1 = a2d1,
b2 = b2, b2sp = b2sp, b2d1 = b2d1,
data = data01,
time = t, censor = censored, infected_treatment = g,
d1 = 'Fréchet', d2 = 'Weibull')
m04_prep_function <- function(a1, b1, a2, a2sp, a2d1, b2, b2sp, b2d1){
nll_basic4(a1, b1, a2, a2sp, a2d1, b2, b2sp, b2d1)
}
m04 <- mle2(m04_prep_function,
start = list(
a1 = 2.35, b1 = 0.66,
a2 = 2, a2sp = 0, a2d1 = 0,
b2 = 0.5, b2sp = 0, b2d1 = 0)
)
summary(m04)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m04_prep_function, start = list(a1 = 2.35, b1 = 0.66,
#> a2 = 2, a2sp = 0, a2d1 = 0, b2 = 0.5, b2sp = 0, b2d1 = 0))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 2.178448 0.065615 33.2003 < 2.2e-16 ***
#> b1 0.613280 0.047463 12.9212 < 2.2e-16 ***
#> a2 2.369249 0.332842 7.1182 1.093e-12 ***
#> a2sp -0.275297 0.328486 -0.8381 0.401986
#> a2d1 0.083912 0.031086 2.6994 0.006947 **
#> b2 -0.537586 0.105418 -5.0996 3.404e-07 ***
#> b2sp 0.727421 0.105199 6.9147 4.687e-12 ***
#> b2d1 -0.024774 0.018365 -1.3490 0.177339
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1506.256
### Model 5
# NB here analyse sporulating vs unsporulating hosts
# i.e. pool uninfected + sporulation = 0 hosts together
# so for index of 'infected_treatment' use 'Sporulation'
nll_basic5 <- nll_basic
body(nll_basic5)[[4]] <- substitute(pfa2 <-
ifelse(data01$dose == 1, a2 + a2d1,
ifelse(data01$dose == 3, a2 - a2d1,
a2))
)
formals(nll_basic5) <- alist(
a1 = a1, b1 = b1,
a2 = a2, a2d1 = a2d1, b2 = b2,
data = data01,
time = t,
censor = censored,
infected_treatment = Sporulation,
d1 = 'Fréchet', d2 = 'Weibull')
m05_prep_function <- function(a1, b1, a2, a2d1, b2){
nll_basic5(a1, b1, a2, a2d1, b2)
}
m05 <- mle2(m05_prep_function,
start = list(a1 = 2, b1 = 1, a2 = 2, a2d1 = 0.1, b2 = 0.5)
)
summary(m05)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m05_prep_function, start = list(a1 = 2, b1 = 1,
#> a2 = 2, a2d1 = 0.1, b2 = 0.5))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 2.121152 0.041035 51.6916 <2e-16 ***
#> b1 0.575389 0.034000 16.9234 <2e-16 ***
#> a2 2.099621 0.027668 75.8858 <2e-16 ***
#> a2d1 0.069198 0.031473 2.1986 0.0279 *
#> b2 0.196991 0.016106 12.2312 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1510.239
AICc(m01, m02, m03, m04, m05, nobs = 328)
#> AICc df
#> 1 1573.651 8
#> 2 1569.711 6
#> 3 1525.581 6
#> 4 1522.707 8
#> 5 1520.426 5Data from [5,6]
data01 <- data_parker
head(data01)
#> Genotype SD Fecundity Sporulation Status Time dose censored d t g
#> 1 721 16 46 1 1 6 2 0 1 6 1
#> 2 721 0 56 0 1 7 0 0 1 7 0
#> 3 721 8 68 0 0 18 1 1 0 18 1
#> 4 721 8 73 0 0 18 1 1 0 18 1
#> 5 721 0 57 0 1 8 0 0 1 8 0
#> 6 721 8 60 0 1 10 1 0 1 10 1
# Infection status known
# Here a host's infection status is defined by whether
# it had visible signs of sporulation at the time of its death
# or right-censoring
# non-sporulating hosts are assumed to only experience same background
# mortality as uninfected hosts
# this is equivalent to taking infection_treatment = Sporulation
m01_prep_function <- function(a1, b1, a2, b2){
nll_basic(a1, b1, a2, b2,
data = data01,
time = t,
censor = censored,
infected_treatment = Sporulation,
d1 = 'Fréchet', d2 = 'Weibull')
}
m01 <- mle2(m01_prep_function,
start = list(a1 = 2.56, b1 = 0.72, a2 = 2, b2 = 0.5)
)
summary(m01)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m01_prep_function, start = list(a1 = 2.56, b1 = 0.72,
#> a2 = 2, b2 = 0.5))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 2.115934 0.041024 51.578 < 2.2e-16 ***
#> b1 0.573568 0.033886 16.927 < 2.2e-16 ***
#> a2 2.092844 0.027533 76.011 < 2.2e-16 ***
#> b2 0.199005 0.016864 11.801 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1515.303
# Infection status unknown
# the infection status of individual hosts is not always known
# the observed survival data may suggest an infected treatment
# harboured exposed-infected and exposed-uninfected hosts
# nll_exposed_infected estimates the proportion of hosts in
# an infected treatment experiencing increased rates of mortality
# due to infection (p1), and
# the proportion experiencing only background mortality (1 - p1)
m02_prep_function <- function(a1, b1, a2, b2, p1){
nll_exposed_infected(a1, b1, a2, b2, p1,
data = data01,
time = t,
censor = censored,
infected_treatment = g,
d1 = 'Frechet', d2 = 'Weibull')
}
m02 <- mle2(m02_prep_function,
start = list(a1 = 2, b1 = 1, a2 = 2, b2 = 0.5, p1 = 0.5)
)
summary(m02)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m02_prep_function, start = list(a1 = 2, b1 = 1,
#> a2 = 2, b2 = 0.5, p1 = 0.5))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 2.204858 0.058196 37.8868 < 2.2e-16 ***
#> b1 0.531848 0.037554 14.1622 < 2.2e-16 ***
#> a2 1.882150 0.050340 37.3891 < 2.2e-16 ***
#> b2 0.166504 0.026792 6.2146 5.145e-10 ***
#> p1 0.479973 0.072172 6.6504 2.922e-11 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1581.703
# in the Parker et al data, the proportion of sporulating hosts in the
# infected treatments were, 119 / (119 + 127) = 0.484
aggregate(data01, by = list(data01$g, data01$Sporulation), length)
#> Group.1 Group.2 Genotype SD Fecundity Sporulation Status Time dose censored
#> 1 0 0 82 82 82 82 82 82 82 82
#> 2 1 0 127 127 127 127 127 127 127 127
#> 3 1 1 119 119 119 119 119 119 119 119
#> d t g
#> 1 82 82 82
#> 2 127 127 127
#> 3 119 119 119
# extend the above to the proportion sporulating within dose treatments
nll_exposed_infected2 <- nll_exposed_infected
body(nll_exposed_infected2)[[6]] <- substitute(pfp1 <-
p1 + ifelse(data01$g == 1 & data01$dose == 1, -p1d,
ifelse(data01$g == 1 & data01$dose == 3, + p1d,
ifelse(data01$g == 1 & data01$dose == 2, 0,
exp(0)
)))
)
formals(nll_exposed_infected2) <- alist(
a1 = a1, b1 = b1, a2 = a2, b2 = b2,
p1 = p1, p1d = p1d,
data = data01,
time = t,
censor = censored,
infected_treatment = g,
d1 = 'Frechet', d2 = 'Weibull'
)
m03_prep_function <- function(a1, b1, a2, b2, p1, p1d){
nll_exposed_infected2(a1, b1, a2, b2, p1, p1d)
}
m03 <- mle2(m03_prep_function,
start = list(a1 = 2.2, b1 = 0.53, a2 = 1.88, b2 = 0.17,
p1 = 0.48, p1d = 0)
)
summary(m03)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m03_prep_function, start = list(a1 = 2.2, b1 = 0.53,
#> a2 = 1.88, b2 = 0.17, p1 = 0.48, p1d = 0))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 2.220077 0.058715 37.8113 < 2.2e-16 ***
#> b1 0.543872 0.038349 14.1823 < 2.2e-16 ***
#> a2 1.914102 0.050109 38.1989 < 2.2e-16 ***
#> b2 0.183066 0.026799 6.8311 8.429e-12 ***
#> p1 0.513907 0.066363 7.7438 9.646e-15 ***
#> p1d 0.243259 0.052237 4.6569 3.211e-06 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1561.351
# estimated proportions infected;
# dose 1 = 0.51 - 0.24 = 0.27
# dose 2 = 0.51 + 0 = 0.51
# dose 3 = 0.51 + 0.24 = 0.75
# observed proportions sporulating;
aggregate(data01, by = list(data01$g, data01$Sporulation, data01$dose), length)
#> Group.1 Group.2 Group.3 Genotype SD Fecundity Sporulation Status Time dose
#> 1 0 0 0 82 82 82 82 82 82 82
#> 2 1 0 1 56 56 56 56 56 56 56
#> 3 1 1 1 25 25 25 25 25 25 25
#> 4 1 0 2 42 42 42 42 42 42 42
#> 5 1 1 2 40 40 40 40 40 40 40
#> 6 1 0 3 29 29 29 29 29 29 29
#> 7 1 1 3 54 54 54 54 54 54 54
#> censored d t g
#> 1 82 82 82 82
#> 2 56 56 56 56
#> 3 25 25 25 25
#> 4 42 42 42 42
#> 5 40 40 40 40
#> 6 29 29 29 29
#> 7 54 54 54 54
# dose 1 = 25 / (25 + 56) = 0.31
# dose 2 = 40 / (40 + 42) = 0.49
# dose 3 = 54 / (54 + 29) = 0.65
AICc (m01, m02, m03, nobs = 328)
#> AICc df
#> 1 1523.427 4
#> 2 1591.890 5
#> 3 1573.613 6Data from [3,4]
data01 <- data_lorenz
m01_prep_function <- function(
a1 = a1, b1 = b1, a2 = a2, b2 = b2){
nll_basic(
a1 = a1, b1 = b1, a2 = a2, b2 = b2,
data = data01,
time = t,
censor = censored,
infected_treatment = g,
d1 = "Gumbel",
d2 = "Frechet"
)}
m01 <- mle2(m01_prep_function,
start = list(a1 = 20, b1 = 5, a2 = 3, b2 = 0.5)
)
summary(m01)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m01_prep_function, start = list(a1 = 20, b1 = 5,
#> a2 = 3, b2 = 0.5))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 22.693577 0.553906 40.9701 < 2.2e-16 ***
#> b1 4.760078 0.329663 14.4392 < 2.2e-16 ***
#> a2 2.843594 0.026058 109.1272 < 2.2e-16 ***
#> b2 0.217803 0.021846 9.9698 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1515.153
m02_prep_function <- function(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, theta = theta){
nll_frailty(
a1 = a1, b1 = b1, a2 = a2, b2 = b2, theta = theta,
data = data_lorenz,
time = t,
censor = censored,
infected_treatment = g,
d1 = "Gumbel",
d2 = "Weibull",
d3 = "Gamma"
)}
m02 <- mle2(m02_prep_function,
start = list(a1 = 20, b1 = 5, a2 = 3, b2 = 0.1, theta = 2)
)
summary(m02)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m02_prep_function, start = list(a1 = 20, b1 = 5,
#> a2 = 3, b2 = 0.1, theta = 2))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 22.789421 0.599282 38.0279 < 2.2e-16 ***
#> b1 4.763390 0.339709 14.0220 < 2.2e-16 ***
#> a2 2.857107 0.036650 77.9566 < 2.2e-16 ***
#> b2 0.090079 0.019450 4.6314 3.632e-06 ***
#> theta 2.618774 1.093566 2.3947 0.01663 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1514.97
AICc(m01, m02, nobs = 256)
#> AICc df
#> 1 1523.312 4
#> 2 1525.210 5Data from [3,4]
data01 <- data_lorenz
m01_prep_function <- function(a1, b1, a2, b2, theta){
nll_frailty_shared(a1, b1, a2, b2, theta,
data = data01,
time = t,
censor = censored,
infected_treatment = g,
d1 = "Gumbel", d2 = "Gumbel")
}
m01 <- mle2(m01_prep_function,
start = list(a1 = 23, b1 = 5, a2 = 10, b2 = 1, theta = 1),
method = "Nelder-Mead",
control = list(maxit = 5000)
)
summary(m01)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m01_prep_function, start = list(a1 = 23, b1 = 5,
#> a2 = 10, b2 = 1, theta = 1), method = "Nelder-Mead", control = list(maxit = 5000))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 22.89581 0.77752 29.4472 < 2.2e-16 ***
#> b1 3.60055 0.45224 7.9616 1.698e-15 ***
#> a2 18.72795 0.70588 26.5313 < 2.2e-16 ***
#> b2 3.34756 0.42925 7.7986 6.259e-15 ***
#> theta 0.34593 0.15975 2.1655 0.03035 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1526.578
m02_prep_function <- function(a1, b1, a2, b2, theta01, theta02, rho){
nll_frailty_correlated(a1, b1, a2, b2, theta01, theta02, rho,
data = data01,
time = t,
censor = censored,
infected_treatment = g,
d1 = "Gumbel",
d2 = "Gumbel")
}
m02 <- mle2(m02_prep_function,
start = list(
a1 = 20, b1 = 5, a2 = 20, b2 = 4, theta01 = 1, theta02 = 1, rho = 1),
method = "L-BFGS-B",
lower = list(
a1 = 1e-6, b1 = 1e-6, a2 = 1e-6, b2 = 1e-6,
theta01 = 1e-6, theta02 = 1e-6, rho = 1e-6)
)
summary(m02)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m02_prep_function, start = list(a1 = 20, b1 = 5,
#> a2 = 20, b2 = 4, theta01 = 1, theta02 = 1, rho = 1), method = "L-BFGS-B",
#> lower = list(a1 = 1e-06, b1 = 1e-06, a2 = 1e-06, b2 = 1e-06,
#> theta01 = 1e-06, theta02 = 1e-06, rho = 1e-06))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 22.880954 NA NA NA
#> b1 4.773423 NA NA NA
#> a2 16.914601 NA NA NA
#> b2 1.267076 NA NA NA
#> theta01 0.000001 NA NA NA
#> theta02 3.857864 NA NA NA
#> rho 1.043373 NA NA NA
#>
#> -2 log L: 1515.177
# NB no standard errors estimated and estimate of 'theta01' at lower boundary
# rerun model with theta01 set at lower boundary
m02b <- mle2(m02_prep_function,
start = list(
a1 = 20, b1 = 5, a2 = 20, b2 = 4,
theta01 = 1, theta02 = 1, rho = 1),
fixed = list(theta01 = 1e-6),
method = "L-BFGS-B",
lower = list(
a1 = 1e-6, b1 = 1e-6, a2 = 1e-6, b2 = 1e-6,
theta02 = 1e-6, rho = 1e-6)
)
summary(m02b)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m02_prep_function, start = list(a1 = 20, b1 = 5,
#> a2 = 20, b2 = 4, theta01 = 1, theta02 = 1, rho = 1), method = "L-BFGS-B",
#> fixed = list(theta01 = 1e-06), lower = list(a1 = 1e-06, b1 = 1e-06,
#> a2 = 1e-06, b2 = 1e-06, theta02 = 1e-06, rho = 1e-06))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 22.87859 0.62613 36.5399 < 2.2e-16 ***
#> b1 4.77634 0.35150 13.5886 < 2.2e-16 ***
#> a2 16.91544 0.70807 23.8896 < 2.2e-16 ***
#> b2 1.26642 0.32547 3.8911 9.98e-05 ***
#> theta02 3.86031 1.80024 2.1443 0.03201 *
#> rho 0.99938 834.99077 0.0012 0.99905
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1515.177
# NB standard error of 'rho' crosses zero (0)
# rerun model with rho set to lower limit
m02c <- mle2(m02_prep_function,
start = list(
a1 = 20, b1 = 5, a2 = 20, b2 = 4,
theta01 = 1, theta02 = 1, rho = 1),
fixed = list(theta01 = 1e-6, rho = 1e-6),
method = "L-BFGS-B",
lower = list(
a1 = 1e-6, b1 = 1e-6, a2 = 1e-6, b2 = 1e-6,
theta02 = 1e-6)
)
summary(m02c)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m02_prep_function, start = list(a1 = 20, b1 = 5,
#> a2 = 20, b2 = 4, theta01 = 1, theta02 = 1, rho = 1), method = "L-BFGS-B",
#> fixed = list(theta01 = 1e-06, rho = 1e-06), lower = list(a1 = 1e-06,
#> b1 = 1e-06, a2 = 1e-06, b2 = 1e-06, theta02 = 1e-06))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 22.87892 0.61671 37.0983 < 2.2e-16 ***
#> b1 4.77634 0.34855 13.7035 < 2.2e-16 ***
#> a2 16.91596 0.64164 26.3635 < 2.2e-16 ***
#> b2 1.26657 0.32409 3.9081 9.303e-05 ***
#> theta02 3.86069 1.57823 2.4462 0.01444 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1515.176
# result of m02c corresponds with estimates from the 'nll_frailty' model,
# where it is assumed there is no unobserved variation in the rate of background mortality
# and where the gamma distribution describes the unobserved variation in virulence
m03_prep_function <- function(a1, b1, a2, b2, theta){
nll_frailty(a1, b1, a2, b2, theta,
data = data01, time = t,
censor = censored, infected_treatment = g,
d1 = "Gumbel", d2 = "Gumbel", d3 = "Gamma")
}
m03 <- mle2(m03_prep_function,
start = list(a1 = 20, b1 = 4, a2 = 20, b2 = 4, theta = 1)
)
summary(m03)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m03_prep_function, start = list(a1 = 20, b1 = 4,
#> a2 = 20, b2 = 4, theta = 1))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 22.87860 0.61668 37.0999 < 2.2e-16 ***
#> b1 4.77647 0.34855 13.7037 < 2.2e-16 ***
#> a2 16.91591 0.64161 26.3649 < 2.2e-16 ***
#> b2 1.26648 0.32402 3.9087 9.28e-05 ***
#> theta 3.86109 1.57816 2.4466 0.01442 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1515.176
AICc(m02c, m03, nobs = 256)
#> AICc df
#> 1 1525.416 5
#> 2 1525.416 5
coef(m02c)
#> a1 b1 a2 b2 theta01 theta02 rho
#> 22.878916 4.776338 16.915958 1.266575 0.000001 3.860685 0.000001
coef(m03)
#> a1 b1 a2 b2 theta
#> 22.878598 4.776473 16.915907 1.266480 3.8610871. Agnew P. 2019 Estimating virulence from relative survival. bioRxiv (doi:10.1101/530709)
2. Blanford S, Jenkins NE, Read AF, Thomas MB. 2012 Evaluating the lethal and pre-lethal effects of a range of fungi against adult anopheles stephensi mosquitoes. Malaria Journal 11, 10. (doi:10.1186/1475-2875-11-365)
3. Lorenz LM, Koella JC. 2011 The microsporidian parasite vavraia culicis as a potential late life-acting control agent of malaria. Evolutionary Applications 4, 783–790. (doi:10.1111/j.1752-4571.2011.00199.x)
4. Lorenz LM, Koella JC. 2011 Data from: The microsporidian parasite vavraia culicis as a potential late life-acting control agent of malaria. Dryad Digital Repository (doi:10.5061/dryad.2s231)
5. Parker BJ, Garcia JR, Gerardo NM. 2014 Genetic variation in resistance and fecundity tolerance in a natural host-pathogen interaction. Evolution 68, 2421–2429. (doi:10.1111/evo.12418)
6. Parker BJ, Garcia JR, Gerardo NM. 2014 Data from: Genetic variation in resistance and fecundity tolerance in a natural host-pathogen interaction. Dryad Digital Repository (doi:10.5061/dryad.24gq7)
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.