Worked examples I

library(anovir)

Introduction

The code below repeats analyses found in the supplementary files of Agnew 2019 [1].

S06 | S08 | S10 | S11 | S12 | S14 | S15

S06. Analysis of Blanford et al. data (i)

Data from [2]

data01 <- subset(data_blanford,
(data_blanford$block == 3) & ((data_blanford$treatment == 'cont') | (data_blanford$treatment == 'Bb06')) & (data_blanford$day > 0)
)

#>     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 4 back to top S08. Analysis of the Lorenz & Koella data Data 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 5 back to top S10. Analysis of Blanford et al. data (ii) Data 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.1232615 back to top S11. Analysis of Parker et al. data Data 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 5 back to top S12. Exposed-but-uninfected hosts model Data 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  6

S14. Lorenz & Koella pooled data

Data 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  5

S15. Shared and correlated frailty models

Data 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),
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.861087

References

1. 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)