Nothing
################################################################################
### Demo of hhh4() modelling of influenza in Southern Germany - data("fluBYBW")
### based on
###
### Paul, M. and Held, L. (2011): Predictive assessment of a non-linear random
### effects model for multivariate time series of infectious disease counts.
### Statistics in Medicine, 30, 1118-1136.
###
### RUNNING THE WHOLE SCRIPT TAKES ~20 MINUTES!
###
### Copyright (C) 2009-2012 Michaela Paul, 2012-2013,2016-2019 Sebastian Meyer
###
### This file is part of the R package "surveillance",
### free software under the terms of the GNU General Public License, version 2,
### a copy of which is available at http://www.r-project.org/Licenses/.
################################################################################
library("surveillance")
## Reproducibility
set.seed(1) # affects initial values for ri() terms
options(digits = 5) # avoid machine-specific output
## Weekly counts of influenza in 140 districts of Bavaria and Baden-Wuerttemberg
data("fluBYBW") # data corrected in surveillance 1.6-0
# -> minor differences to original results in the paper
##################################################
# Fit the models from the Paul & Held (2011) paper
##################################################
## generate formula for temporal and seasonal trends
f.end <- addSeason2formula(f = ~ -1 + ri(type="iid", corr="all") +
I((t-208)/100), S=3, period=52)
## settings for the optimizer
opt <- list(stop = list(tol=1e-5, niter=200),
regression = list(method="nlminb"),
variance = list(method="nlminb"))
## models
# A0
cntrl_A0 <- list(ar = list(f = ~ -1),
end = list(f =f.end, offset = population(fluBYBW)),
family = "NegBin1", optimizer = opt, verbose = 1)
summary(res_A0 <- hhh4(fluBYBW,cntrl_A0))
# B0
cntrl_B0 <- list(ar = list(f = ~ 1),
end = list(f =f.end, offset = population(fluBYBW)),
family = "NegBin1", optimizer = opt, verbose=1)
res_B0 <- hhh4(fluBYBW,cntrl_B0)
# C0
cntrl_C0 <- list(ar = list(f = ~ -1 + ri(type="iid", corr="all")),
end = list(f =f.end, offset = population(fluBYBW)),
family = "NegBin1", optimizer = opt, verbose=1)
res_C0 <- hhh4(fluBYBW,cntrl_C0)
#A1
# weight matrix w_ji = 1/(No. neighbors of j) if j ~ i, and 0 otherwise
wji <- neighbourhood(fluBYBW)/rowSums(neighbourhood(fluBYBW))
cntrl_A1 <- list(ar = list(f = ~ -1),
ne = list(f = ~ 1, weights = wji),
end = list(f =f.end, offset = population(fluBYBW)),
family = "NegBin1", optimizer = opt, verbose=1)
res_A1 <- hhh4(fluBYBW,cntrl_A1)
# B1
cntrl_B1 <- list(ar = list(f = ~ 1),
ne = list(f = ~ 1, weights = wji),
end = list(f =f.end, offset = population(fluBYBW)),
family = "NegBin1", optimizer = opt, verbose=1)
res_B1 <- hhh4(fluBYBW,cntrl_B1)
# C1
cntrl_C1 <- list(ar = list(f = ~ -1 + ri(type="iid", corr="all")),
ne = list(f = ~ 1, weights = wji),
end = list(f =f.end, offset = population(fluBYBW)),
family = "NegBin1", optimizer = opt, verbose=1)
res_C1 <- hhh4(fluBYBW,cntrl_C1)
#A2
cntrl_A2 <- list(ar = list(f = ~ -1),
ne = list(f = ~ -1 + ri(type="iid",corr="all"), weights=wji),
end = list(f =f.end, offset = population(fluBYBW)),
family = "NegBin1", optimizer = opt, verbose=1)
res_A2 <- hhh4(fluBYBW,cntrl_A2)
# B2
cntrl_B2 <- list(ar = list(f = ~ 1),
ne = list(f = ~ -1 + ri(type="iid",corr="all"), weights =wji),
end = list(f =f.end, offset = population(fluBYBW)),
family = "NegBin1", optimizer = opt, verbose=1)
res_B2 <- hhh4(fluBYBW,cntrl_B2)
# C2
cntrl_C2 <- list(ar = list(f = ~ -1 + ri(type="iid", corr="all")),
ne = list(f = ~ -1 + ri(type="iid",corr="all"), weights =wji),
end = list(f =f.end, offset = population(fluBYBW)),
family = "NegBin1", optimizer = opt, verbose=1,
start=list(fixed=fixef(res_B0),random=c(rep(0,140),
ranef(res_B0)), sd.corr=c(-.5,res_B0$Sigma.orig,0)))
res_C2 <- hhh4(fluBYBW,cntrl_C2)
# D
cntrl_D <- list(ar = list(f = ~ 1),
ne = list(f = ~ -1 + ri(type="iid"), weights = wji),
end = list(f =addSeason2formula(f = ~ -1 + ri(type="car") +
I((t-208)/100), S=3, period=52),
offset = population(fluBYBW)),
family = "NegBin1", optimizer = opt, verbose=1)
res_D <- hhh4(fluBYBW,cntrl_D)
###########################################################
## Exemplary summary of model B2
## (compare with Paul & Held, 2011, Table III and Figure 5)
###########################################################
summary(res_B2, idx2Exp = 1:2, maxEV = TRUE)
## Note: as of surveillance 1.6-0, results differ slightly from the paper
## (see penalized log-likelihood), because a superfluous row of zeros
## has been removed from the fluBYBW data
.idx <- c(113, 111, 46, 77)
plot(res_B2, units = .idx, names = fluBYBW@map@data[.idx, "name"],
legend = 2, legend.args = list(x = "topleft"), legend.observed = TRUE)
######################################################################
# Compare the predictive performance of the models by computing
# one-step-ahead predictions to be assessed by proper scoring rules
######################################################################
## do 1-step ahead predictions for the last two years
tp <- nrow(fluBYBW)-2*52
verbose <- interactive()
## for this demo: only calculate pseudo-predictions based on the final fit
## to avoid the time-consuming sequential refitting at each step.
TYPE <- "final"
## use "rolling" for true one-step-ahead predictions => TAKES ~8 HOURS!
val_A0 <- oneStepAhead(res_A0, tp=tp, type=TYPE, verbose=verbose)
val_B0 <- oneStepAhead(res_B0, tp=tp, type=TYPE, verbose=verbose)
val_C0 <- oneStepAhead(res_C0, tp=tp, type=TYPE, verbose=verbose)
val_A1 <- oneStepAhead(res_A1, tp=tp, type=TYPE, verbose=verbose)
val_B1 <- oneStepAhead(res_B1, tp=tp, type=TYPE, verbose=verbose)
val_C1 <- oneStepAhead(res_C1, tp=tp, type=TYPE, verbose=verbose)
val_A2 <- oneStepAhead(res_A2, tp=tp, type=TYPE, verbose=verbose)
val_B2 <- oneStepAhead(res_B2, tp=tp, type=TYPE, verbose=verbose)
val_C2 <- oneStepAhead(res_C2, tp=tp, type=TYPE, verbose=verbose)
val_D <- oneStepAhead(res_D, tp=tp, type=TYPE, verbose=verbose)
## compute scores
vals <- ls(pattern="val_")
nam <- substring(vals,first=5,last=6)
whichScores <- c("logs", "rps", "ses")
scores_i <- vector(mode="list", length=length(vals))
meanScores <- NULL
for(i in seq_along(vals)){
sc <- scores(get(vals[i]), which=whichScores, individual=TRUE, reverse=TRUE)
## reverse=TRUE => same permutation test results as in surveillance < 1.16.0
scores_i[[i]] <- sc
meanScores <- rbind(meanScores,colMeans(sc, dims=2))
}
names(scores_i) <- nam
rownames(meanScores) <- nam
print(meanScores)
## Note that the above use of "final" fitted values instead of "rolling"
## one-step-ahead predictions leads to different mean scores than reported
## in Paul & Held (2011, Table IV).
## assess statistical significance of score differences
compareWithBest <- function(best, whichModels, nPermut=9999, seed=1234){
set.seed(seed)
pVals <- NULL
for(score in seq_along(whichScores)){
p <- c()
for(model in whichModels){
p <- c(p, if(model==best) NA else
permutationTest(scores_i[[model]][,,score],scores_i[[best]][,,score],
plot=interactive(),nPermutation=nPermut, verbose=TRUE)$pVal.permut)
}
pVals <- cbind(pVals,p)
}
return(pVals)
}
pVals_flu <- compareWithBest(best=9, whichModels=1:10,
nPermut=999, # reduced for this demo
seed=2059710987)
rownames(pVals_flu) <- nam
colnames(pVals_flu) <- whichScores
print(pVals_flu)
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.