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

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

back to top

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

back to top

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

back to top

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)