The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
library(anovir)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)
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.