data("imdepi", "imdepifit")
## load borders of Germany's districts (originally obtained from
## the German Federal Agency for Cartography and Geodesy,
## https://gdz.bkg.bund.de/), simplified by the "modified Visvalingam"
## algorithm (level=6.6%) using MapShaper.org (v. 0.1.17):
load(system.file("shapes", "districtsD.RData", package="surveillance"))
if (surveillance.options("allExamples")) {
plot(districtsD)
plot(stateD, add=TRUE, border=2, lwd=2)
# 'stateD' was obtained as 'rgeos::gUnaryUnion(districtsD)'
}
## simulate 2 realizations (over a short period, for speed)
## considering events from data(imdepi) before t=31 as prehistory
mysims <- simulate(imdepifit, nsim=2, seed=1, data=imdepi,
tiles=districtsD, newcoef=c("e.typeC"=-1),
t0=31, T=if (interactive()) 180 else 45, # for CRAN
simplify=TRUE)
# \dontshow{
## check construction and selection from "simEpidataCSlist"
local({
mysim_from_list <- mysims[[1]]
capture.output(mysim_single <- eval("[[<-"(attr(mysims, "call"), "nsim", 1)))
mysim_from_list$runtime <- mysim_single$runtime <- NULL
stopifnot(all.equal(mysim_single, mysim_from_list,
check.attributes = FALSE))
})
## check equivalence of Lambdag from simulation and residuals via twinstim
stopifnot(all.equal(
residuals(mysims[[1]]),
suppressMessages(surveillance:::residuals.twinstim(surveillance:::as.twinstim.simEpidataCS(mysims[[1]])))
))
# }
## plot both simulations using the plot-method for simEpidataCSlist's
mysims
plot(mysims, aggregate="time")
## extract the second realization -> object of class simEpidataCS
mysim2 <- mysims[[2]]
summary(mysim2)
plot(mysim2, aggregate="space")
if (surveillance.options("allExamples")) {
### compare the observed _cumulative_ number of cases during the
### first 90 days to 20 simulations from the fitted model
sims <- simulate(imdepifit, nsim=20, seed=1, data=imdepi, t0=0, T=90,
tiles=districtsD, simplify=TRUE)
## extract cusums
getcsums <- function (events) {
tapply(events$time, events@data["type"],
function (t) cumsum(table(t)), simplify=FALSE)
}
csums_observed <- getcsums(imdepi$events)
csums_simulated <- lapply(sims$eventsList, getcsums)
## plot it
plotcsums <- function (csums, ...) {
mapply(function (csum, ...) lines(as.numeric(names(csum)), csum, ...),
csums, ...)
invisible()
}
plot(c(0,90), c(0,35), type="n", xlab="Time [days]",
ylab="Cumulative number of cases")
plotcsums(csums_observed, col=c(2,4), lwd=3)
legend("topleft", legend=levels(imdepi$events$type), col=c(2,4), lwd=1)
invisible(lapply(csums_simulated, plotcsums,
col=adjustcolor(c(2,4), alpha=0.5)))
}
if (FALSE) {
### Experimental code to generate 'nsim' simulations of 'nm2add' months
### beyond the observed time period:
nm2add <- 24
nsim <- 5
### The events still infective by the end of imdepi$stgrid will be used
### as the prehistory for the continued process.
origT <- tail(imdepi$stgrid$stop, 1)
## create a time-extended version of imdepi
imdepiext <- local({
## first we have to expand stgrid (assuming constant "popdensity")
g <- imdepi$stgrid
g$stop <- g$BLOCK <- NULL
gadd <- data.frame(start=rep(seq(origT, by=30, length.out=nm2add),
each=nlevels(g$tile)),
g[rep(seq_len(nlevels(g$tile)), nm2add), -1])
## now create an "epidataCS" using this time-extended stgrid
as.epidataCS(events=imdepi$events, # the replacement warnings are ok
W=imdepi$W, qmatrix=imdepi$qmatrix,
stgrid=rbind(g, gadd), T=max(gadd$start) + 30)
})
newT <- tail(imdepiext$stgrid$stop, 1)
## simulate beyond the original period
simsext <- simulate(imdepifit, nsim=nsim, seed=1, t0=origT, T=newT,
data=imdepiext, tiles=districtsD, simplify=TRUE)
## Aside to understand the note from checking events and tiles:
# marks(imdepi)["636",] # tile 09662 is attributed to this event, but:
# plot(districtsD[c("09678","09662"),], border=1:2, lwd=2, axes=TRUE)
# points(imdepi$events["636",])
## this mismatch is due to polygon simplification
## plot the observed and simulated event numbers over time
plot(imdepiext, breaks=c(unique(imdepi$stgrid$start),origT),
cumulative=list(maxat=330))
for (i in seq_along(simsext$eventsList))
plot(simsext[[i]], add=TRUE, legend.types=FALSE,
breaks=c(unique(simsext$stgrid$start),newT),
subset=!is.na(source), # have to exclude the events of the prehistory
cumulative=list(offset=c(table(imdepi$events$type)), maxat=330, axis=FALSE),
border=NA, density=0) # no histogram
abline(v=origT, lty=2, lwd=2)
}
Run the code above in your browser using DataLab