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)
}
## simulate 2 realizations (over a short period, for speed)
## considering events from data(imdepi) before t=31 as prehistory
## IGNORE_RDIFF_BEGIN
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)
## IGNORE_RDIFF_END
# \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")
### compare the observed _cumulative_ number of cases in the first 90 days to
nsim <- 20
### simulations from the fitted model
if (!interactive()) nsim <- 2
sims <- simulate(imdepifit, nsim=nsim, 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.f=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)
## extend the 'stgrid' by replicating the last block 'nm2add' times
## (i.e., holding "popdensity" constant)
stgridext <- local({
gLast <- subset(imdepi$stgrid, BLOCK == max(BLOCK))
gAdd <- gLast[rep(1:nrow(gLast), nm2add),]; rownames(gAdd) <- NULL
newstart <- seq(origT, by=30, length.out=nm2add)
newstop <- c(newstart[-1], max(newstart) + 30)
gAdd$start <- rep(newstart, each=nlevels(gAdd$tile))
gAdd$stop <- rep(newstop, each=nlevels(gAdd$tile))
rbind(imdepi$stgrid, gAdd, make.row.names = FALSE)[,-1]
})
## create an updated "epidataCS" with the time-extended 'stgrid'
imdepiext <- update(imdepi, stgrid = stgridext)
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