### 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=frequency(measles2w)*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)