# Load invasive meningococcal disease data
data("imdepi")
### first, fit a simple endemic-only model
m_noepi <- twinstim(
endemic = addSeason2formula(~ offset(log(popdensity)) + I(start/365-3.5),
S=1, period=365, timevar="start"),
data = imdepi, subset = !is.na(agegrp)
)
## look at the model summary
summary(m_noepi)
## there is no evidence for a type-dependent endemic intercept (LR test)
m_noepi_type <- update(m_noepi, endemic = ~(1|type) + .)
pchisq(2*c(logLik(m_noepi_type)-logLik(m_noepi)), df=1, lower.tail=FALSE)
### add an epidemic component with just the intercept, i.e.
### assuming uniform dispersal in time and space up to a distance of
### eps.s = 200 km and eps.t = 30 days (see summary(imdepi))
m0 <- update(m_noepi, epidemic=~1, start=c("e.(Intercept)"=-18), model=TRUE)
## summarize the model fit
s <- summary(m0, correlation = TRUE, symbolic.cor = TRUE)
s
# output the table of coefficients as LaTeX code
toLatex(s, digits=2)
# or, to report rate ratios
xtable(s)
## the default confint-method can be used for Wald-CI's
confint(m0, level=0.95)
## same "untrimmed" R0 for every event (simple epidemic intercept model)
summary(R0(m0, trimmed=FALSE))
## plot the path of the fitted total intensity
plot(m0, "total intensity", tgrid=500)
## the remaining examples are more time-consuming
if (surveillance.options("allExamples"))
{
## NB: in contrast to using nlminb() optim's BFGS would miss the
## likelihood maximum wrt the epidemic intercept if starting at -10
m0_BFGS <- update(m_noepi, epidemic=~1, start=c("e.(Intercept)"=-10),
optim.args = list(method="BFGS"))
format(cbind(nlminb=coef(m0), BFGS=coef(m0_BFGS)), digits=3, scientific=FALSE)
m0_BFGS$fisherinfo # singular Fisher information matrix here
m0$fisherinfo
logLik(m0_BFGS)
logLik(m0)
## nlminb is more powerful since we make use of the analytical fisherinfo
## as estimated by the model during optimization, which optim cannot
## extract "residual process" integrating over space (takes some seconds)
res <- residuals(m0)
# if the model describes the true CIF well _in the temporal dimension_,
# then this residual process should behave like a stationary Poisson
# process with intensity 1
plot(res, type="l"); abline(h=c(0, length(res)), lty=2)
# -> use the function checkResidualProcess
checkResidualProcess(m0)
## estimate an exponential temporal decay of infectivity
m1_tiaf <- update(m0, tiaf=tiaf.exponential())
plot(m1_tiaf, "tiaf", scaled=FALSE)
## estimate a step function for spatial interaction
m1_fstep <- update(m0, siaf=c(5,10,25,50))
plot(m1_fstep, "siaf", scaled=FALSE)
## estimate a continuously decreasing spatial interaction function, where
## likelihood evaluation takes much longer than for (piecewise) constant
## spatial spread; here we use the kernel of an isotropic bivariate
## Gaussian with a standard deviation of exp(2.8) = 16.4 km for both
## finetypes (fixed to speed up the example)
m1 <- update(m0,
siaf = siaf.gaussian(1, F.adaptive = TRUE),
start = c("e.siaf.1" = 2.8), # sigma = exp(2.8) = 16.4 km
optim.args = list(fixed="e.siaf.1"),
control.siaf = list(F=list(adapt=0.25)) # use adaptive eps=sigma/4
)
AIC(m_noepi, m0, m1) # further improvement
summary(m1, test.iaf=FALSE) # nonsense to test H0: log(sigma) = 0
plot(m1, "siaf", scaled=FALSE)
## add epidemic covariates
m2 <- update(m1, epidemic = ~ 1 + type + agegrp)
AIC(m1, m2) # further improvement
summary(m2)
## look at estimated R0 values by event type
tapply(R0(m2), imdepi$events@data[names(R0(m2)), "type"], summary)
}Run the code above in your browser using DataLab