Simulate "hhh4"
Count Time Series
hhh4_simulate.Rd
Simulates a multivariate time series of counts based on the Poisson/Negative Binomial model as described in Paul and Held (2011).
Arguments
- object
an object of class
"hhh4"
.- nsim
number of time series to simulate. Defaults to
1
.- seed
an object specifying how the random number generator should be initialized for simulation (via
set.seed
). The initial state will also be stored as an attribute"seed"
of the result. The original state of the.Random.seed
will be restored at the end of the simulation. By default (NULL
), neither initialization nor recovery will be done. This behaviour is copied from thesimulate.lm
method.- y.start
vector or matrix (with
ncol(object$stsObj)
columns) with starting counts for the epidemic components. IfNULL
, the observed means in the respective units of the data inobject
duringsubset
are used.- subset
time period in which to simulate data. Defaults to (and cannot exceed) the whole period defined by the underlying
"sts"
object.- coefs
coefficients used for simulation from the model in
object
. Default is to use the fitted parameters. Note that thecoefs
-vector must be in the same order and scaling ascoef(object)
, which especially meansreparamPsi = TRUE
(as per default when using thecoef
-method to extract the parameters). The overdispersion parameter incoefs
is the inverse of the dispersion parametersize
inrnbinom
.- components
character vector indicating which components of the fitted model
object
should be active during simulation. For instance, a simulation withcomponents="end"
is solely based on the fitted endemic mean.- simplify
logical indicating if only the simulated counts (
TRUE
) or the full"sts"
object (FALSE
) should be returned for every replicate. By default a full"sts"
object is returned iffnsim=1
.- ...
unused (argument of the generic).
Details
Simulates data from a Poisson or a Negative Binomial model
with mean
$$\mu_{it} = \lambda_{it} y_{i,t-1} +
\phi_{it} \sum_{j \neq i} w_{ji} y_{j,t-1} +
\nu_{it}$$
where
\(\lambda_{it}>0\), \(\phi_{it}>0\), and \(\nu_{it}>0\) are
parameters which are modelled parametrically.
The function uses the model and parameter estimates of the fitted
object
to simulate the time series.
With the argument coefs
it is possible to simulate from
the model as specified in object
, but with different
parameter values.
Value
If simplify=FALSE
: an object of class
"sts"
(nsim = 1
) or a list of those
(nsim > 1
).
If simplify=TRUE
: an object of class
"hhh4sims"
, which is an array of dimension
c(length(subset), ncol(object$stsObj), nsim)
.
The originally observed counts during the simulation period,
object$stsObj[subset,]
, are attached for reference
(used by the plot
-methods) as an attribute "stsObserved"
,
and the initial condition y.start
as attribute "initial"
.
The [
-method for "hhh4sims"
takes care of subsetting
these attributes appropriately.
References
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
See also
plot.hhh4sims
and scores.hhh4sims
and the examples therein for nsim > 1
.
Examples
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
# 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")))
# 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))