data(influMen)
# convert to sts class and extract meningococcal disease time series
meningo <- disProg2sts(influMen)[,2]
# fit model
fit <- hhh4(meningo, control = list(
ar = list(f = ~ 1),
end = list(f = addSeason2formula(~1, period = 52)),
family = "NegBin1"))
plot(fit)
# simulate from model (generates an "sts" object)
simData <- simulate(fit, seed=1234)
# plot simulated data
plot(simData, main = "simulated data", xaxis.labelFormat=NULL)
# use simplify=TRUE to return an array of simulated counts
simCounts <- simulate(fit, seed=1234, simplify=TRUE)
dim(simCounts) # nTime x nUnit x nsim
stopifnot(observed(simData) == c(simCounts))
# plot the first year of simulated counts (+ initial + observed)
plot(simCounts[1:52,,], type = "time", xaxis.labelFormat = NULL)
# see help(plot.hhh4sims) for other plots, mainly useful for nsim > 1
# simulate from a Poisson instead of a NegBin model
# keeping all other parameters fixed at their original estimates
coefs <- replace(coef(fit), "overdisp", 0)
simData2 <- simulate(fit, seed=123, coefs = coefs)
plot(simData2, main = "simulated data: Poisson model", xaxis.labelFormat = NULL)
# simulate from a model with higher autoregressive parameter
coefs <- replace(coef(fit), "ar.1", log(0.9))
simData3 <- simulate(fit, seed=321, coefs = coefs)
plot(simData3, main = "simulated data: lambda = 0.5", xaxis.labelFormat = NULL)
## more sophisticated: simulate beyond initially observed time range
# extend data range by one year (non-observed domain), filling with NA values
nextend <- 52
timeslots <- c("observed", "state", "alarm", "upperbound", "populationFrac")
addrows <- function (mat, n) mat[c(seq_len(nrow(mat)), rep(NA, n)),,drop=FALSE]
extended <- Map(function (x) addrows(slot(meningo, x), n = nextend), x = timeslots)
# create new sts object with extended matrices
meningo2 <- do.call("sts", c(list(start = meningo@start, frequency = meningo@freq,
map = meningo@map), extended))
# fit to the observed time range only, via the 'subset' argument
fit2 <- hhh4(meningo2, control = list(
ar = list(f = ~ 1),
end = list(f = addSeason2formula(~1, period = 52)),
family = "NegBin1",
subset = 2:(nrow(meningo2) - nextend)))
# the result is the same as before
stopifnot(all.equal(fit, fit2, ignore = c("stsObj", "control")))
# \dontshow{
# one-week-ahead prediction only "works" for the first non-observed time point
# because the autoregressive component relies on non-missing past counts
oneStepAhead(fit2, tp = rep(nrow(meningo2)-nextend, 2), type = "final", verbose = FALSE)
# however, methods won't work as observed is NA
# }
# long-term probabilistic forecast via simulation for non-observed time points
meningoSim <- simulate(fit2, nsim = 100, seed = 1,
subset = seq(nrow(meningo)+1, nrow(meningo2)),
y.start = tail(observed(meningo), 1))
apply(meningoSim, 1:2, function (ysim) quantile(ysim, c(0.1, 0.5, 0.9)))
# three plot types are available for "hhh4sims", see also ?plot.hhh4sims
plot(meningoSim, type = "time", average = median)
plot(meningoSim, type = "size", observed = FALSE)
if (requireNamespace("fanplot"))
plot(meningoSim, type = "fan", means.args = list(),
fan.args = list(ln = c(.1,.9), ln.col = 8))
Run the code above in your browser using DataLab