# Load invasive meningococcal disease data
data("imdepi")
### first, fit a simple endemic-only model
## Calculate initial value for the endemic intercept (the crude estimate
## which is used by default if no initial value is supplied).
## Assume the model only had a single-cell endemic component
## (rate of homogeneous Poisson process scaled for population density):
popdensity.overall <- with(subset(imdepi$stgrid, BLOCK == 1),
weighted.mean(popdensity, area))
popdensity.overall # pop. density in Germany is ~230 inhabitants/km^2
W.area <- with(subset(imdepi$stgrid, BLOCK == 1), sum(area))
W.area # area of Germany is about 357100 km^2
# this should actually be the same as
sum(sapply(imdepi$W@polygons, slot, "area"))
# which here is not exactly the case because of polygon simplification
## initial value for the endemic intercept
h.intercept <- with(summary(imdepi),
log(nEvents/popdensity.overall/diff(timeRange)/W.area))
## fit the endemic-only model
m_noepi <- twinstim(
endemic = ~ 1 + offset(log(popdensity)) + I(start/365) +
sin(start * 2 * pi/365) + cos(start * 2 * pi/365),
data = imdepi, subset = !is.na(agegrp),
optim.args = list(par=c(h.intercept,rep(0,3)), # this is the default
method="nlminb", control = list(trace=1)),
model = FALSE, cumCIF = FALSE # for reasons of speed
)
## look at the model summary
summary(m_noepi)
if (surveillance.options("allExamples"))
{
## type-dependent endemic intercept?
m_noepi_type <- update(m_noepi, endemic = ~(1|type) + .)
summary(m_noepi_type)
# LR-test for a type-dependent intercept in the endemic-only model
pchisq(2*(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)
## NOTE: 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(coef(m0), 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
## 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, withAIC=FALSE)
## the default confint-method can be used for Wald-CI's
confint(m0, level=0.95)
## same R0 estimate for every event (epidemic intercept only)
summary(R0(m0, trimmed=FALSE))
## plot the path of the fitted total intensity
plot(m0, "total intensity", tgrid=500)
if (surveillance.options("allExamples"))
{
## 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)
## NB: the model fit environment is kept in the workspace with model=TRUE
sort(sapply(environment(m0), object.size))
# saving m0 to RData will include this model environment, thus might
# consume quiet a lot of disk space (as it does use memory), mainly
# depending on the size of "stgrid", the number of "events" and the
# polygonal resolution of "W"
}
if (surveillance.options("allExamples"))
{
### try with spatially decreasing interaction kernel
## Likelihood evaluation takes much longer than for constant spatial spread
## Here we use the kernel of an isotropic bivariate Gaussian with same
## standard deviation for both finetypes
m1 <- update(m0,
siaf = siaf.gaussian(1, F.adaptive = TRUE),
start = c("e.siaf.1" = 2.8), # exp(2.8) = sigma = 16 km
optim.args = list(fixed="e.siaf.1"),
# for reasons of speed of the example, the siaf
# parameter log(sigma) is fixed to 2.8 here
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", xlim=c(0,200), types=1) # the same for types=2
### 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