surveillance (version 1.5-4)

simulate.ah4: Simulates data based on the model proposed by Paul and Held (2011)

Description

Simulates a multivariate time series of counts based on the Poisson/Negative Binomial model as described in Paul and Held (2011).

Usage

## S3 method for class 'ah4':
simulate(object, nsim = 1, seed = NULL, y.start = NULL, coefs = NULL, ...)

Arguments

object

an object of class "ah4".

nsim

number of time series to simulate. Defaults to 1.

seed

an integer that will be used in the call to set.seed before simulating the time series.

y.start

vector with starting counts for the epidemic components. If NULL, the observed means in the respective units of the data in object are used.

coefs

if not NULL, these values (in the same order as coef(object)) are used to simulate from the model specified in object.

...

unused (argument of the generic).

Value

• An object of class "sts" in the case of nsim = 1, and a list of "sts" objects otherwise.

encoding

latin1

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.

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

hhh4, simHHH

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(S = 1, period = 52)),
            family = "NegBin1"))
plot(fit)

# simulate from model
set.seed(1234)
simData <- simulate(fit)

# plot simulated data
plot(simData, main = "simulated data", legend.opts = NULL, xaxis.years = FALSE)

# consider a Poisson instead of a NegBin model
coefs <- coef(fit)
coefs["overdisp"] <- 0

simData2 <- simulate(fit, coefs = coefs)
plot(simData2, main = "simulated data: Poisson model", 
     legend.opts = NULL, xaxis.years = FALSE)

# consider a model with higher autoregressive parameter
coefs <- coef(fit)
coefs[1] <- log(0.5)

simData3 <- simulate(fit, coefs = coefs)
plot(simData3, main = "simulated data: lambda = 0.5", 
     legend.opts = NULL, xaxis.years = FALSE)