Nothing
## ----include = FALSE---------------------------------------------------------------
## load the "cool" package
library("surveillance")
## Compute everything or fetch cached results?
message("Doing computations: ",
COMPUTE <- !file.exists("hhh4_spacetime-cache.RData"))
if (!COMPUTE) load("hhh4_spacetime-cache.RData", verbose = TRUE)
## ----measlesWeserEms_components, echo=FALSE----------------------------------------
## extract components from measlesWeserEms to reconstruct
data("measlesWeserEms")
counts <- observed(measlesWeserEms)
map <- measlesWeserEms@map
populationFrac <- measlesWeserEms@populationFrac
## ----measlesWeserEms_neighbourhood, eval=FALSE-------------------------------------
# weserems_adjmat <- poly2adjmat(map)
# weserems_nbOrder <- nbOrder(weserems_adjmat)
## ----echo=FALSE--------------------------------------------------------------------
weserems_nbOrder <- measlesWeserEms@neighbourhood
## ----measlesWeserEms_construct-----------------------------------------------------
measlesWeserEms <- sts(counts, start = c(2001, 1), frequency = 52,
population = populationFrac, neighbourhood = weserems_nbOrder, map = map)
## ----measlesWeserEms, fig.cap="Measles infections in the Weser-Ems region, 2001--2002.", fig.subcap=c("Time series of weekly counts.","Disease incidence (per 100\\,000 inhabitants)."), fig.width=5, fig.height=5, out.width="0.5\\linewidth", fig.pos="htb", echo=-1----
par(mar = c(5,5,1,1))
plot(measlesWeserEms, type = observed ~ time)
plot(measlesWeserEms, type = observed ~ unit,
population = measlesWeserEms@map$POPULATION / 100000,
labels = list(font = 2), colorkey = list(space = "right"),
sp.layout = layout.scalebar(measlesWeserEms@map, corner = c(0.05, 0.05),
scale = 50, labels = c("0", "50 km"), height = 0.03))
## ----measlesWeserEms15, fig.cap=paste("Count time series of the", sum(colSums(observed(measlesWeserEms))>0), "affected districts."), out.width="\\linewidth", fig.width=10, fig.height=6, fig.pos="htb"----
if (require("ggplot2")) {
autoplot(measlesWeserEms, units = which(colSums(observed(measlesWeserEms)) > 0))
} else plot(measlesWeserEms, units = which(colSums(observed(measlesWeserEms)) > 0))
## ----measlesWeserEms_animation, eval=FALSE-----------------------------------------
# animation::saveHTML(
# animate(measlesWeserEms, tps = 1:52, total.args = list()),
# title = "Evolution of the measles epidemic in the Weser-Ems region, 2001",
# ani.width = 500, ani.height = 600)
## ----echo=FALSE, eval=FALSE--------------------------------------------------------
# ## to perform the following analysis using biweekly aggregated measles counts:
# measlesWeserEms <- aggregate(measlesWeserEms, by = "time", nfreq = 26)
## ----measlesModel_basic------------------------------------------------------------
measlesModel_basic <- list(
end = list(f = addSeason2formula(~1 + t, period = measlesWeserEms@freq),
offset = population(measlesWeserEms)),
ar = list(f = ~1),
ne = list(f = ~1, weights = neighbourhood(measlesWeserEms) == 1),
family = "NegBin1")
## ----measlesFit_basic--------------------------------------------------------------
measlesFit_basic <- hhh4(stsObj = measlesWeserEms, control = measlesModel_basic)
## ----measlesFit_basic_summary------------------------------------------------------
summary(measlesFit_basic, idx2Exp = TRUE, amplitudeShift = TRUE, maxEV = TRUE)
## ----measlesFit_basic_endseason, fig.width=6, fig.height=2.5, out.width=".5\\linewidth", fig.cap="Estimated multiplicative effect of seasonality on the endemic mean.", fig.pos="ht"----
plot(measlesFit_basic, type = "season", components = "end", main = "")
## ----measlesFitted_basic, fig.cap="Fitted components in the initial model \\code{measlesFit\\_basic} for the five districts with more than 50 cases as well as summed over all districts (bottom right). Dots are only drawn for positive weekly counts.", out.width="\\linewidth", fig.pos="htb"----
districts2plot <- which(colSums(observed(measlesWeserEms)) > 50)
par(mfrow = c(2,3), mar = c(3, 5, 2, 1), las = 1)
plot(measlesFit_basic, type = "fitted", units = districts2plot,
hide0s = TRUE, par.settings = NULL, legend = 1)
plot(measlesFit_basic, type = "fitted", total = TRUE,
hide0s = TRUE, par.settings = NULL, legend = FALSE) -> fitted_components
## ----------------------------------------------------------------------------------
fitted_components$Overall[20:22,]
## ----------------------------------------------------------------------------------
colSums(fitted_components$Overall)[3:5] / sum(fitted_components$Overall[,1])
## ----------------------------------------------------------------------------------
confint(measlesFit_basic, parm = "overdisp")
## ----measlesFit_basic_Poisson------------------------------------------------------
AIC(measlesFit_basic, update(measlesFit_basic, family = "Poisson"))
## ----Sprop-------------------------------------------------------------------------
Sprop <- matrix(1 - measlesWeserEms@map@data$vacc1.2004,
nrow = nrow(measlesWeserEms), ncol = ncol(measlesWeserEms), byrow = TRUE)
summary(Sprop[1, ])
## ----SmodelGrid--------------------------------------------------------------------
Soptions <- c("unchanged", "Soffset", "Scovar")
SmodelGrid <- expand.grid(end = Soptions, ar = Soptions)
row.names(SmodelGrid) <- do.call("paste", c(SmodelGrid, list(sep = "|")))
## ----measlesFits_vacc, eval=COMPUTE------------------------------------------------
# measlesFits_vacc <- apply(X = SmodelGrid, MARGIN = 1, FUN = function (options) {
# updatecomp <- function (comp, option) switch(option, "unchanged" = list(),
# "Soffset" = list(offset = comp$offset * Sprop),
# "Scovar" = list(f = update(comp$f, ~. + log(Sprop))))
# update(measlesFit_basic,
# end = updatecomp(measlesFit_basic$control$end, options[1]),
# ar = updatecomp(measlesFit_basic$control$ar, options[2]),
# data = list(Sprop = Sprop))
# })
## ----aics_vacc, eval=COMPUTE-------------------------------------------------------
# aics_vacc <- do.call(AIC, lapply(names(measlesFits_vacc), as.name),
# envir = as.environment(measlesFits_vacc))
## ----------------------------------------------------------------------------------
aics_vacc[order(aics_vacc[, "AIC"]), ]
## ----measlesFit_vacc---------------------------------------------------------------
measlesFit_vacc <- update(measlesFit_basic,
end = list(f = update(formula(measlesFit_basic)$end, ~. + log(Sprop))),
data = list(Sprop = Sprop))
coef(measlesFit_vacc, se = TRUE)["end.log(Sprop)", ]
## ----------------------------------------------------------------------------------
2^cbind("Estimate" = coef(measlesFit_vacc),
confint(measlesFit_vacc))["end.log(Sprop)",]
## ----measlesFit_nepop--------------------------------------------------------------
measlesFit_nepop <- update(measlesFit_vacc,
ne = list(f = ~log(pop)), data = list(pop = population(measlesWeserEms)))
## ----------------------------------------------------------------------------------
cbind("Estimate" = coef(measlesFit_nepop),
confint(measlesFit_nepop))["ne.log(pop)",]
## ----measlesFit_powerlaw-----------------------------------------------------------
measlesFit_powerlaw <- update(measlesFit_nepop,
ne = list(weights = W_powerlaw(maxlag = 5)))
## ----------------------------------------------------------------------------------
cbind("Estimate" = coef(measlesFit_powerlaw),
confint(measlesFit_powerlaw))["neweights.d",]
## ----measlesFit_np-----------------------------------------------------------------
measlesFit_np2 <- update(measlesFit_nepop,
ne = list(weights = W_np(maxlag = 2)))
## ----measlesFit_neweights, fig.width=5, fig.height=3.5, fig.cap="Estimated weights as a function of adjacency order.", out.width="0.47\\linewidth", fig.subcap=c("Normalized power-law weights.", "Non-normalized weights with 95\\% CIs."), echo=c(1,5)----
library("lattice")
trellis.par.set("reference.line", list(lwd=3, col="gray"))
trellis.par.set("fontsize", list(text=14))
set.seed(20200303)
plot(measlesFit_powerlaw, type = "neweights", plotter = stripplot,
panel = function (...) {panel.stripplot(...); panel.average(...)},
jitter.data = TRUE, xlab = expression(o[ji]), ylab = expression(w[ji]))
## non-normalized weights (power law and unconstrained second-order weight)
local({
colPL <- "#0080ff"
ogrid <- 1:5
par(mar=c(3.6,4,2.2,2), mgp=c(2.1,0.8,0))
plot(ogrid, ogrid^-coef(measlesFit_powerlaw)["neweights.d"], col=colPL,
xlab="Adjacency order", ylab="Non-normalized weight", type="b", lwd=2)
matlines(t(sapply(ogrid, function (x)
x^-confint(measlesFit_powerlaw, parm="neweights.d"))),
type="l", lty=2, col=colPL)
w2 <- exp(c(coef(measlesFit_np2)["neweights.d"],
confint(measlesFit_np2, parm="neweights.d")))
lines(ogrid, c(1,w2[1],0,0,0), type="b", pch=19, lwd=2)
arrows(x0=2, y0=w2[2], y1=w2[3], length=0.1, angle=90, code=3, lty=2)
legend("topright", col=c(colPL, 1), pch=c(1,19), lwd=2, bty="n", inset=0.1, y.intersp=1.5,
legend=c("Power-law model", "Second-order model"))
})
## ----------------------------------------------------------------------------------
AIC(measlesFit_nepop, measlesFit_powerlaw, measlesFit_np2)
## ----measlesFit_ri, results="hide"-------------------------------------------------
measlesFit_ri <- update(measlesFit_powerlaw,
end = list(f = update(formula(measlesFit_powerlaw)$end, ~. + ri() - 1)),
ar = list(f = update(formula(measlesFit_powerlaw)$ar, ~. + ri() - 1)),
ne = list(f = update(formula(measlesFit_powerlaw)$ne, ~. + ri() - 1)))
## ----measlesFit_ri_summary_echo, eval=FALSE----------------------------------------
# summary(measlesFit_ri, amplitudeShift = TRUE, maxEV = TRUE)
## ----------------------------------------------------------------------------------
head(ranef(measlesFit_ri, tomatrix = TRUE), n = 3)
## ----measlesFit_ri_map, out.width="0.31\\linewidth", fig.width=3.5, fig.height=3.7, fig.pos="htb", fig.cap="Estimated multiplicative effects on the three components.", fig.subcap=c("Autoregressive", "Spatio-temporal", "Endemic")----
for (comp in c("ar", "ne", "end")) {
print(plot(measlesFit_ri, type = "ri", component = comp, exp = TRUE,
labels = list(cex = 0.6)))
}
## ----------------------------------------------------------------------------------
exp(ranef(measlesFit_ri, intercept = TRUE)["03403", "ar.ri(iid)"])
## ----measlesFitted_ri, out.width="\\linewidth", fig.pos="htb", fig.cap="Fitted components in the random effects model \\code{measlesFit\\_ri} for the five districts with more than 50 cases as well as summed over all districts. Compare to Figure~\\ref{fig:measlesFitted_basic}."----
par(mfrow = c(2,3), mar = c(3, 5, 2, 1), las = 1)
plot(measlesFit_ri, type = "fitted", units = districts2plot,
hide0s = TRUE, par.settings = NULL, legend = 1)
plot(measlesFit_ri, type = "fitted", total = TRUE,
hide0s = TRUE, par.settings = NULL, legend = FALSE)
## ----measlesFitted_maps, fig.cap="Maps of the fitted component proportions averaged over all weeks.", fig.pos="hbt", fig.width=10, fig.height=3.7, out.width="0.93\\linewidth"----
plot(measlesFit_ri, type = "maps",
which = c("epi.own", "epi.neighbours", "endemic"),
prop = TRUE, labels = list(cex = 0.6))
## ----measlesPreds1, results="hide"-------------------------------------------------
tp <- c(65, 77)
models2compare <- paste0("measlesFit_", c("basic", "powerlaw", "ri"))
measlesPreds1 <- lapply(mget(models2compare), oneStepAhead,
tp = tp, type = "final")
## ----echo=FALSE--------------------------------------------------------------------
stopifnot(all.equal(measlesPreds1$measlesFit_powerlaw$pred,
fitted(measlesFit_powerlaw)[tp[1]:tp[2],],
check.attributes = FALSE))
## ----include=FALSE-----------------------------------------------------------------
stopifnot(all.equal(
measlesFit_powerlaw$loglikelihood,
-sum(scores(oneStepAhead(measlesFit_powerlaw, tp = 1, type = "final"),
which = "logs", individual = TRUE))))
## ----measlesScores1----------------------------------------------------------------
SCORES <- c("logs", "rps", "dss", "ses")
measlesScores1 <- lapply(measlesPreds1, scores, which = SCORES, individual = TRUE)
t(sapply(measlesScores1, colMeans, dims = 2))
## ----measlesPreds2, eval=COMPUTE, results="hide"-----------------------------------
# measlesPreds2 <- lapply(mget(models2compare), oneStepAhead,
# tp = tp, type = "rolling", which.start = "final")
## ----measlesPreds2_plot, fig.cap = "Fan charts of rolling one-week-ahead forecasts during the second quarter of 2002, as produced by the random effects model \\code{measlesFit\\_ri}, for the five most affected districts.", out.width="\\linewidth", echo=-1----
par(mfrow = sort(n2mfrow(length(districts2plot))), mar = c(4.5,4.5,2,1))
for (unit in names(districts2plot))
plot(measlesPreds2[["measlesFit_ri"]], unit = unit, main = unit,
key.args = if (unit == tail(names(districts2plot),1)) list())
## ----measlesScores2----------------------------------------------------------------
measlesScores2 <- lapply(measlesPreds2, scores, which = SCORES, individual = TRUE)
t(sapply(measlesScores2, colMeans, dims = 2))
## ----measlesScores_test------------------------------------------------------------
set.seed(321)
sapply(SCORES, function (score) permutationTest(
measlesScores2$measlesFit_ri[, , score],
measlesScores2$measlesFit_basic[, , score],
nPermutation = 999))
## ----measlesPreds2_calibrationTest_echo, eval=FALSE--------------------------------
# calibrationTest(measlesPreds2[["measlesFit_ri"]], which = "rps")
## ----measlesPreds2_pit, fig.width=8, fig.height=2.5, out.width="0.93\\linewidth", fig.cap="PIT histograms of competing models to check calibration of the one-week-ahead predictions during the second quarter of 2002.", echo=-1, fig.pos="hbt"----
par(mfrow = sort(n2mfrow(length(measlesPreds2))), mar = c(4.5,4.5,2,1), las = 1)
for (m in models2compare)
pit(measlesPreds2[[m]], plot = list(ylim = c(0, 1.25), main = m))
## ----measlesFit_powerlaw2, include = FALSE-----------------------------------------
## a simplified model which includes the autoregression in the power law
measlesFit_powerlaw2 <- update(measlesFit_powerlaw,
ar = list(f = ~ -1),
ne = list(weights = W_powerlaw(maxlag = 5, from0 = TRUE)))
AIC(measlesFit_powerlaw, measlesFit_powerlaw2)
## simpler is really worse; probably needs random effects
## ----measlesFit_ri_simulate--------------------------------------------------------
(y.start <- observed(measlesWeserEms)[52, ])
measlesSim <- simulate(measlesFit_ri,
nsim = 100, seed = 1, subset = 53:104, y.start = y.start)
## ----------------------------------------------------------------------------------
summary(colSums(measlesSim, dims = 2))
## ----measlesSim_plot_time, fig.cap="Simulation-based long-term forecast starting from the last week in 2001 (left-hand dot). The plot shows the weekly counts aggregated over all districts. The fan chart represents the 1\\% to 99\\% quantiles of the simulations in each week; their mean is displayed as a white line. The circles correspond to the observed counts.", fig.pos="htb"----
plot(measlesSim, "fan", means.args = list(), key.args = list())
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.