library(anovir)This example shows nll_proportional_virulence extended to comparing multiple treatments against a reference treatment for virulence.
Data from [1,2]
data01 <- data_lorenz
# create column 'ref_vir' in dataframe coded 0, 1 and 2
# for control, reference dose, and other dose treatments, respectively
data01$ref_vir <- ifelse(data01$g == 0, 0,
ifelse(data01$g == 1, ifelse(data01$Infectious.dose == 5000, 1, 2), 0)
)
# copy and modify 'nll_proportional_virulence' function
# where the virulence in the reference treatment is 'theta'
# which is multiplied by 'th10', 'th20', etc... in the other dose treatements
nll_proportional_virulence2 <- nll_proportional_virulence
body(nll_proportional_virulence2)[[6]] <- substitute(pfth <- theta *
ifelse(data01$Infectious.dose == 10000, th10,
ifelse(data01$Infectious.dose == 20000, th20,
ifelse(data01$Infectious.dose == 40000, th40,
ifelse(data01$Infectious.dose == 80000, th80,
ifelse(data01$Infectious.dose == 160000, th160,
exp(0)
)))))
)
# update formals
formals(nll_proportional_virulence2) <- alist(
a1 = a1, b1 = b1, a2 = a2, b2 = b2,
theta = theta,
th10 = th10, th20 = th20, th40 = th40, th80 = th80, th160 = th160,
data = data01,
time = t,
censor = censored,
infected_treatment = g,
reference_virulence = ref_vir,
d1 = 'Gumbel', d2 = 'Weibull')
# write 'prep function'
m01_prep_function <- function(
a1, b1, a2, b2, theta, th10, th20, th40, th80, th160){
nll_proportional_virulence2(
a1, b1, a2, b2, theta, th10, th20, th40, th80, th160)
}
# send 'prep function' to mle2
# NB fixing 'theta = 1' scales virulence in other treatments relative to
# that in the reference treatment
m01 <- mle2(m01_prep_function,
start = list(
a1 = 24, b1 = 5, a2 = 4, b2 = 0.2,
theta = 1,
th10 = 1, th20 = 1, th40 = 1, th80 = 1, th160 = 1),
fixed = list(theta = 1)
)
summary(m01)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m01_prep_function, start = list(a1 = 24, b1 = 5,
#> a2 = 4, b2 = 0.2, theta = 1, th10 = 1, th20 = 1, th40 = 1,
#> th80 = 1, th160 = 1), fixed = list(theta = 1))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 23.165445 0.652943 35.4785 < 2e-16 ***
#> b1 4.660055 0.384772 12.1112 < 2e-16 ***
#> a2 3.205421 0.084396 37.9807 < 2e-16 ***
#> b2 0.188362 0.020123 9.3607 < 2e-16 ***
#> th10 2.280418 1.142473 1.9960 0.04593 *
#> th20 2.303080 1.144282 2.0127 0.04415 *
#> th40 3.649006 1.798433 2.0290 0.04246 *
#> th80 4.590818 2.277040 2.0161 0.04379 *
#> th160 5.513346 2.706341 2.0372 0.04163 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1504.575The values of a2 and b2 define virulence at time t in the reference dose treatment of 5 (x103) spores larva-1 as,
\[\begin{equation} h_{5000}(t) = \frac{1}{0.188 t} \exp \left[ \frac{\log t - 3.205}{0.188} \right] \end{equation}\]which was estimated as being multiplied by 2.280, 2.303, 3.649, 4.591 and 5.513 times in the dose treatments of 10, 20, 40, 80 and 160 (x103) spores larva-1, respectively, assuming virulence is proportional among dose treatments.
This example shows how to manipulate the function nll_basic_logscale, which assumes location and scale parameters are on a logscale.
This function can avoid the generation of warning messages associated with parameters with negative values being given to logarithmic functions.
Data from [1,2]
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
nll_basic_logscale2 <- nll_basic_logscale
body(nll_basic_logscale2)[[4]]
#> pfa2 <- exp(a2)Here the location parameter a2 is to be made a linear function of log(dose),
a2 = a2i - a2ii*log(data01$Infectious.dose)
where a2i and a2ii are parameters to be estimated.
However instead of directly replacing a2 with the right hand side of the expression above, i.e.,
pfa2 <- exp(a2i - a2ii*log(data01$Infectious.dose))
both parameters need exponentiating, i.e.,
pfa2 <- exp(a2i) - exp(a2ii)*log(data01$Infectious.dose)
body(nll_basic_logscale2)[[4]] <- substitute(
pfa2 <- exp(a2i) - exp(a2ii)*log(data01$Infectious.dose)
)
body(nll_basic_logscale2)[[4]]
#> pfa2 <- exp(a2i) - exp(a2ii) * log(data01$Infectious.dose)
formals(nll_basic_logscale2) <- alist(
a1 = a1, b1 = b1, a2i = a2i, a2ii = a2ii, b2 = b2,
data = data01,
time = t,
censor = censored,
infected_treatment = g,
d1 = 'Gumbel', d2 = 'Weibull'
)
m01_prep_function <- function(a1, b1, a2i, a2ii, b2){
nll_basic_logscale2(a1, b1, a2i, a2ii, b2)
}
m01 <- mle2(m01_prep_function,
start = list(a1 = 3, b1 = 1.5, a2i = 1.4, a2ii = -3, b2 = -1.5)
)
summary(m01)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m01_prep_function, start = list(a1 = 3, b1 = 1.5,
#> a2i = 1.4, a2ii = -3, b2 = -1.5))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 3.143306 0.028334 110.937 < 2.2e-16 ***
#> b1 1.550561 0.081129 19.112 < 2.2e-16 ***
#> a2i 1.346133 0.049796 27.033 < 2.2e-16 ***
#> a2ii -2.511533 0.217948 -11.524 < 2.2e-16 ***
#> b2 -1.684226 0.104964 -16.046 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1506.184
exp(coef(m01))
#> a1 b1 a2i a2ii b2
#> 23.18036585 4.71411429 3.84253844 0.08114371 0.18558808Estimates made using nll_basic are similar, but generate Warning messages
nll_basic2 <- nll_basic
body(nll_basic2)[[4]]
#> pfa2 <- a2
body(nll_basic2)[[4]] <- substitute(
pfa2 <- a2i + a2ii*log(data01$Infectious.dose)
)
body(nll_basic2)[[4]]
#> pfa2 <- a2i + a2ii * log(data01$Infectious.dose)
formals(nll_basic2) <- alist(
a1 = a1, b1 = b1, a2i = a2i, a2ii = a2ii, b2 = b2,
data = data01,
time = t,
censor = censored,
infected_treatment = g,
d1 = 'Gumbel', d2 = 'Weibull'
)
m02_prep_function <- function(a1, b1, a2i, a2ii, b2){
nll_basic2(a1, b1, a2i, a2ii, b2)
}
m02 <- mle2(m02_prep_function,
start = list(a1 = 23, b1 = 4.5, a2i = 3.8, a2ii = -0.1, b2 = 0.4)
)
#> Warning in log(h1 + g * h2): NaNs produced
#> Warning in log(h1 + g * h2): NaNs produced
#> Warning in log(h1 + g * h2): NaNs produced
#> Warning in log(h1 + g * h2): NaNs produced
#> Warning in log(h1 + g * h2): NaNs produced
#> Warning in log(h1 + g * h2): NaNs produced
#> Warning in log(h1 + g * h2): NaNs produced
#> Warning in log(h1 + g * h2): NaNs produced
#> Warning in log(h1 + g * h2): NaNs produced
#> Warning in log(h1 + g * h2): NaNs produced
#> Warning in log(h1 + g * h2): NaNs produced
summary(m02)
#> Maximum likelihood estimation
#>
#> Call:
#> mle2(minuslogl = m02_prep_function, start = list(a1 = 23, b1 = 4.5,
#> a2i = 3.8, a2ii = -0.1, b2 = 0.4))
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(z)
#> a1 23.189717 0.657339 35.2782 < 2.2e-16 ***
#> b1 4.720309 0.383171 12.3191 < 2.2e-16 ***
#> a2i 3.833275 0.190004 20.1747 < 2.2e-16 ***
#> a2ii -0.080280 0.017580 -4.5666 4.956e-06 ***
#> b2 0.185412 0.019389 9.5627 < 2.2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> -2 log L: 1506.1861. 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)
2. 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)