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.

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)

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.