Nothing
### R code from vignette source 'hhh4.Rnw'
###################################################
### code chunk number 1: setup
###################################################
library("surveillance")
options(width=75)
## create directory for plots
dir.create("plots", showWarnings=FALSE)
######################################################
## Do we need to compute or can we just fetch results?
######################################################
compute <- !file.exists("hhh4-cache.RData")
message("Doing computations: ", compute)
if(!compute) load("hhh4-cache.RData")
###################################################
### code chunk number 2: loadInfluMen
###################################################
# load data
data("influMen")
# convert to sts class and print basic information about the time series
print(fluMen <- disProg2sts(influMen))
###################################################
### code chunk number 3: getMen
###################################################
meningo <- fluMen[, "meningococcus"]
dim(meningo)
###################################################
### code chunk number 4: plotfluMen
###################################################
getOption("SweaveHooks")[["fig"]]()
plot(fluMen, type = observed ~ time | unit, # type of plot (default)
same.scale = FALSE, # unit-specific ylim?
col = "grey") # color of bars
###################################################
### code chunk number 5: readInFlu
###################################################
# read in observed number of cases
flu.counts <- as.matrix(read.table(system.file("extdata/counts_flu_BYBW.txt",
package = "surveillance"),
check.names = FALSE))
###################################################
### code chunk number 6: nhoodByBw
###################################################
getOption("SweaveHooks")[["fig"]]()
# read in 0/1 adjacency matrix (1 if regions share a common border)
nhood <- as.matrix(read.table(system.file("extdata/neighbourhood_BYBW.txt",
package = "surveillance"),
check.names = FALSE))
library("Matrix")
print(image(Matrix(nhood)))
###################################################
### code chunk number 7: fluAsSTS
###################################################
# read in population fractions
popfracs <- read.table(system.file("extdata/population_2001-12-31_BYBW.txt",
package = "surveillance"),
header = TRUE)$popFrac
# create sts object
flu <- sts(flu.counts, start = c(2001, 1), frequency = 52,
population = popfracs, neighbourhood = nhood)
###################################################
### code chunk number 8: plot-flu-ByBw
###################################################
getOption("SweaveHooks")[["fig"]]()
data("fluBYBW")
plot(fluBYBW[year(fluBYBW) == 2001, ], # select year 2001
type = observed ~ unit, # total counts by region
population = fluBYBW@map$X31_12_01 / 100000, # per 100000 inhabitants
colorkey = list(title = "Incidence [per 100'000 inhabitants]"))
###################################################
### code chunk number 9: hhh4.Rnw:266-271
###################################################
# consistency check
local({
fluBYBW@map <- flu@map
stopifnot(all.equal(fluBYBW, flu))
})
###################################################
### code chunk number 10: measles2w
###################################################
data("measlesDE")
measles2w <- aggregate(measlesDE, nfreq = 26)
###################################################
### code chunk number 11: plot-measles
###################################################
getOption("SweaveHooks")[["fig"]]()
plot(measles2w, type = observed ~ time, # aggregate counts over all units
main = "Bi-weekly number of measles cases in Germany")
###################################################
### code chunk number 12: hhh4 (eval = FALSE)
###################################################
## hhh4(sts, control)
###################################################
### code chunk number 13: controlObj (eval = FALSE)
###################################################
## control = list(
## ar = list(f = ~ -1, # formula for log(lambda_it)
## offset = 1), # optional multiplicative offset
## ne = list(f = ~ -1, # formula for log(phi_it)
## offset = 1, # optional multiplicative offset
## weights = neighbourhood(stsObj) == 1), # (w_ji) matrix
## end = list(f = ~ 1, # formula for log(nu_it)
## offset = 1), # optional multiplicative offset e_it
## family = "Poisson", # Poisson or NegBin model
## subset = 2:nrow(stsObj), # subset of observations to be used
## optimizer = list(stop = list(tol = 1e-5, niter = 100), # stop rules
## regression = list(method = "nlminb"), # for penLogLik
## variance = list(method = "nlminb")), # for marLogLik
## verbose = FALSE, # level of progress reporting
## start = list(fixed = NULL, # list with initial values for fixed,
## random = NULL, # random, and
## sd.corr = NULL), # variance parameters
## data = list(t = epoch(stsObj)-1),# named list of covariates
## keep.terms = FALSE # whether to keep the model terms
## )
###################################################
### code chunk number 14: fitMeningo0
###################################################
# specify a formula object for the endemic component
( f_S1 <- addSeason2formula(f = ~ 1, S = 1, period = 52) )
# fit the Poisson model
result0 <- hhh4(meningo, control = list(end = list(f = f_S1),
family = "Poisson"))
summary(result0)
###################################################
### code chunk number 15: fitMeningo1
###################################################
result1 <- update(result0, family = "NegBin1")
###################################################
### code chunk number 16: hhh4.Rnw:496-497
###################################################
AIC(result0, result1)
###################################################
### code chunk number 17: fitMeningo2
###################################################
# fit an autoregressive model
result2 <- update(result1, ar = list(f = ~ 1))
###################################################
### code chunk number 18: hhh4.Rnw:510-514
###################################################
coef(result2, se = TRUE, # also return standard errors
amplitudeShift = TRUE, # transform sine/cosine coefficients
# to amplitude/shift parameters
idx2Exp = TRUE) # exponentiate remaining parameters
###################################################
### code chunk number 19: plot_result2
###################################################
getOption("SweaveHooks")[["fig"]]()
plot(result2)
###################################################
### code chunk number 20: neighbourhood_fluMen
###################################################
# no "transmission" from meningococcus to influenza
neighbourhood(fluMen)["meningococcus","influenza"] <- 0
neighbourhood(fluMen)
###################################################
### code chunk number 21: fitFluMen
###################################################
# create formula for endemic component
f.end <- addSeason2formula(f = ~ -1 + fe(1, unitSpecific = TRUE),
# disease-specific intercepts
S = c(3, 1), # S = 3 for flu, S = 1 for men
period = 52)
# specify model
m <- list(ar = list(f = ~ -1 + fe(1, unitSpecific = TRUE)),
ne = list(f = ~ 1, # phi, only relevant for meningococcus due to
weights = neighbourhood(fluMen)), # the weight matrix
end = list(f = f.end),
family = "NegBinM") # disease-specific overdispersion
# fit model
result <- hhh4(fluMen, control = m)
summary(result, idx2Exp=1:3)
###################################################
### code chunk number 22: plot-fit_men
###################################################
getOption("SweaveHooks")[["fig"]]()
plot(result, units = NULL, pch = 20, legend = 2, legend.args = list(
legend = c("influenza-driven", "autoregressive", "endemic")))
###################################################
### code chunk number 23: plot-fit_men_decomposed
###################################################
getOption("SweaveHooks")[["fig"]]()
plot(result, units = NULL, pch = 20, legend = 2,
decompose = TRUE, col = c(7, 4))
###################################################
### code chunk number 24: ri (eval = FALSE)
###################################################
## f.end <- ~ -1 + ri(type = "iid", corr = "all")
###################################################
### code chunk number 25: modelFluBYBW
###################################################
# endemic component: iid random effects, linear trend, S=3 seasonal terms
f.end <- addSeason2formula(f = ~ -1 + ri(type="iid", corr="all") +
I((t-208)/100),
S = 3, period = 52)
# model specification
model.B2 <- list(ar = list(f = ~ 1),
ne = list(f = ~ -1 + ri(type="iid", corr="all"),
weights = neighbourhood(fluBYBW),
normalize = TRUE), # all(rowSums(weights) == 1)
end = list(f = f.end, offset = population(fluBYBW)),
family = "NegBin1", verbose = TRUE,
optimizer = list(variance = list(method = "Nelder-Mead")))
# default start values for random effects are sampled from a normal
set.seed(42)
###################################################
### code chunk number 26: computeFluBYBW
###################################################
if(compute){
result.B2 <- hhh4(fluBYBW, model.B2)
s.B2 <- summary(result.B2, maxEV = TRUE, idx2Exp = 1:3)
#pred.B2 <- oneStepAhead(result.B2, tp = nrow(fluBYBW) - 2*52)
predfinal.B2 <- oneStepAhead(result.B2, tp = nrow(fluBYBW) - 2*52,
type = "final")
meanSc.B2 <- colMeans(scores(predfinal.B2))
save(s.B2, meanSc.B2, file="hhh4-cache.RData")
}
###################################################
### code chunk number 27: fitFluBYBW (eval = FALSE)
###################################################
## # fit the model (takes about 35 seconds)
## result.B2 <- hhh4(fluBYBW, model.B2)
## summary(result.B2, maxEV = TRUE, idx2Exp = 1:3)
###################################################
### code chunk number 28: hhh4.Rnw:670-671
###################################################
s.B2
###################################################
### code chunk number 29: oneStepAhead_rolling (eval = FALSE)
###################################################
## pred.B2 <- oneStepAhead(result.B2, tp = nrow(fluBYBW) - 2*52)
###################################################
### code chunk number 30: oneStepAhead_fake (eval = FALSE)
###################################################
## predfinal.B2 <- oneStepAhead(result.B2, tp = nrow(fluBYBW) - 2*52,
## type = "final")
###################################################
### code chunk number 31: scores (eval = FALSE)
###################################################
## colMeans(scores(predfinal.B2, which = c("logs", "rps")))
###################################################
### code chunk number 32: hhh4.Rnw:703-704
###################################################
meanSc.B2[c("logs", "rps")]
###################################################
### code chunk number 33: createVacc
###################################################
data(MMRcoverageDE)
cardVac1 <- MMRcoverageDE[1:16,3:4]
adjustVac <- function(cardVac, p=0.5,nrow=1){
card <- cardVac[,1]
vac <- cardVac[,2]
vacAdj <- vac*card + p*vac*(1-card)
return(matrix(vacAdj,nrow=nrow, ncol=length(vacAdj), byrow=TRUE))
}
vac0 <- 1-adjustVac(cardVac1,p=0.5,nrow=measles2w@freq*3)
colnames(vac0) <- colnames(measles2w)
###################################################
### code chunk number 34: hhh4.Rnw:750-751
###################################################
vac0[1:2, 1:6]
###################################################
### code chunk number 35: fitMeasles
###################################################
# endemic component: Intercept + sine/cosine terms
f.end <- addSeason2formula(f = ~ 1, S = 1, period = 26)
# autoregressive component: Intercept + vaccination coverage information
model.A0 <- list(ar = list(f = ~ 1 + logVac0),
end = list(f = f.end, offset = population(measles2w)),
data = list(t = epoch(measles2w), logVac0 = log(vac0)))
# fit the model
result.A0 <- hhh4(measles2w, model.A0)
summary(result.A0, amplitudeShift = TRUE)
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.