# 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 the`simulate.lm`

method.- y.start
vector or matrix (with

`ncol(object$stsObj)`

columns) with starting counts for the epidemic components. If`NULL`

, the observed means in the respective units of the data in`object`

during`subset`

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 the`coefs`

-vector must be in the same order and scaling as`coef(object)`

, which especially means`reparamPsi = TRUE`

(as per default when using the`coef`

-method to extract the parameters). The overdispersion parameter in`coefs`

is the inverse of the dispersion parameter`size`

in`rnbinom`

.- components
character vector indicating which components of the fitted model

`object`

should be active during simulation. For instance, a simulation with`components="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 iff`nsim=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
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))
```