Learn R Programming
surveillance (version 1.4-1)

twinSIR_simulation: Simulation of Epidemic Data

Description

This function simulates the infection (and removal) times of an epidemic. Besides the classical SIR type of epidemic, also SI, SIRS and SIS epidemics are supported. Simulation works via the conditional intensity of infection of an individual, given some (time varying) endemic covariates and/or some distance functions (epidemic components) as well as the fixed positions of the individuals. The lengths of the infectious and removed periods are generated following a pre-specified function (can be deterministic). The simulate method for objects of class "twinSIR" simulates new epidemic data using the model and the parameter estimates of the fitted object.

Usage

simEpidata(formula, data, id.col, I0.col, coords.cols, subset,
           beta, h0, f = list(), alpha, infPeriod,
           remPeriod = function(ids) rep(Inf, length(ids)),
           end = Inf, trace = FALSE, .allocate = 500L)
              
## S3 method for class 'twinSIR':
simulate(object, nsim = 1, seed = 1,
         infPeriod = NULL, remPeriod = NULL,
         end = diff(range(object$intervals)), trace = FALSE, .allocate = 500L,
         data = object$data, ...)

Arguments

formula
an object of class "formula" (or one that can be coerced to that class): a symbolic description of the intensity model to be estimated. The details of model specification are given under Details.
data
a data.frame containing the variables in formula and the variables specified by id.col, I0.col and coords.col (see below). It represents the history of the endemic covariates to
id.col
only if data does not inherit from epidata: single index of the id column in data. Can be numeric (by column number) or character (by column name). The id column identifies the indi
I0.col
only if data does not inherit from epidata: single index of the I0 column in data. Can be numeric (by column number), character (by column name) or NULL. The I0 column
coords.cols
only if data does not inherit from epidata: indexes of the coords columns in data. Can be a numeric (by column number), a character (by column name) vector or NULL
subset
an optional vector specifying a subset of the covariate history to be used in the simulation.
beta
numeric vector of length equal the number of endemic (cox) terms on the rhs of formula. It contains the effects of the endemic predictor (excluding the log-baseline h0, see below) in the same order as in
h0
either a single number to specify a constant baseline hazard (equal to exp(h0)) or a list of functions named exact and upper. In the latter case, h0$exact is the true log-b
f
a named list of distance functions or list() (the default), if no epidemic component is desired (if additionally infPeriod (see below) always returns Inf, then one simulates from the Cox model).
alpha
numeric vector of length equal the number of epidemic terms, i.e. the number of distance functions in f. It contains the effects of the epidemic predictor in the same order as the distance functions in f or, if names
infPeriod
a function generating lengths of infectious periods. It should take one parameter (e.g. ids), which is a character vector of id's of individuals, and return appropriate infection periods for those individuals. Therefore, the va
remPeriod
a function generating lengths of removal periods. Per default, once an individual was removed it will stay in this state forever (Inf). Therefore, it will not become at-risk (S) again and re-infections are not possible. Alternat
end
a single positive numeric value specifying the time point at which the simulation should be forced to end. By default, this is Inf, i.e. the simulation continues until there is no susceptible individual left. For the simulate
trace
logical (or integer) indicating if (or how often) the sets of susceptible and infected individuals as well as the rejection indicator (of the rejection sampling step) should be cated. Defaults to FALSE.
.allocate
number of blocks to initially allocate for the event history; defaults to 500, i.e. 500*N rows. Each time the simulated epidemic exceeds the allocated space, the event history will be enlarged by .allocate blocks.
object
an object of class "twinSIR". This must contain the original data used for the fit (see data).
nsim
number of epidemics to simulate. Defaults to 1.
seed
an integer that will be used in the call to set.seed before simulating the epidemics.
...
unused (argument of the generic).

Value

  • An object of class "simEpidata", which is a data.frame with the columns "id", "start", "stop", "atRiskY", "event", "Revent" and the coordinate columns (with the original names from data), which are all obligatory. These columns are followed by all the variables appearing on the rhs of the formula. Last but not least, the generated columns with epidemic variables corresponding to the functions in the list f are appended, if length(f) > 0.

    Note that objects of class "simEpidata" also inherit from class "epidata", thus all "epidata" methods can be applied. The data.frame is given the additional attributes

  • "eventTimes"numeric vector of infection time points (sorted chronologically).
  • "timeRange"numeric vector of length 2: c(min(start), max(stop)).
  • "coords.cols"numeric vector containing the column indices of the coordinate columns in the resulting data-frame.
  • "f"this equals the argument f.
  • "config"a list with elements h0 = h0$exact, beta and alpha.
  • callthe matched call.
  • termsthe terms object used.

encoding

latin1

Details

A model is specified through the formula, which has the form

cbind(start, stop) ~ cox(endemicVar1) * cox(endemicVar2), i.e. the right hand side has the usual form as in lm, but all variables are marked as being endemic by the special function cox. The effects of those predictor terms are specified by beta. The left hand side of the formula denotes the start and stop columns in data. This can be omitted, if data inherits from class "epidata" in which case cbind(start, stop) will be used. The epidemic model component is specified by the parameter f (and by the epidemic effects alpha).

The simulation algorithm used is Ogata's modified thinning. For details, see H�hle{Hoehle} (2009), Section 4.

References

H�hle{Hoehle}, M. (2009), Additive-Multiplicative Regression Models for Spatio-Temporal Epidemics, Biometrical Journal, 51(6):961-978.

See Also

The plot.epidata and animate.epidata methods for plotting and animating (simulated) epidemic data, respectively. The intensityplot.simEpidata method for plotting paths of infection intensities.

Function twinSIR for fitting spatio-temporal epidemic intensity models to epidemic data.

Examples

Run this code
set.seed(1234)
# Generate a data frame containing a hypothetic population with 100 individuals
n <- 100
pos <- matrix(rnorm(n*2), ncol=2, dimnames=list(NULL, c("x", "y")))
pop <- data.frame(id=1:n,x=pos[,1], y=pos[,2], 
                  gender=sample(0:1, n, replace=TRUE),
                  I0col=rep(0,n),start=rep(0,n),stop=rep(Inf,n))

# Simulate an epidemic in this population
epi <- simEpidata(cbind(start,stop) ~ cox(gender), 
                  data = pop, 
                  id = "id", I0.col = "I0col", coords.cols = c("x","y"),
                  beta = c(-3), h0 = -2, alpha = c(B1 = 0.1),
                  f = list(B1 = function(u) u <= 1),
                  infPeriod = function(ids) rexp(length(ids), rate=1))

# Plot the numbers of susceptible, infectious and removed individuals
plot(epi)

# load data of an observed epidemic
data("fooepidata")
summary(fooepidata)

# simulate a new evolution of the epidemic
simepi <- simEpidata(cbind(start, stop) ~ cox(z1) + cox(z1):cox(z2),
    data = fooepidata, 
    beta = c(1,0.5), h0 = -7, alpha = c(B2 = 0.1, B1 = 0.05),
    f = list(B1 = function(u) u<=1, B2 = function(u) u>1 & is.finite(u)),
    infPeriod = function(ids) rexp(length(ids), rate=1), trace = FALSE)
summary(simepi)
plot(simepi)

animate(simepi)

intensityplot(simepi)

# load a fitted model object, which must contain the original data 
# (use 'keep.data = TRUE' in the call of 'twinSIR')
data("foofit")
foofit

# plot original epidemic
plot(foofit$data)

## simulate a new epidemic using the model and parameter estimates of 'foofit'
## and set simulation period = observation period
# with observed infPeriods (i.e. fixed length 3 days):
simfitepi1 <- simulate(foofit, nsim = 1)[[1]]
plot(simfitepi1)
# with new infPeriods (simuluated from the Exp(0.3) distribution):
simfitepi2 <- simulate(foofit, nsim = 1,
                       infPeriod=function(ids) rexp(length(ids), rate=0.3))[[1]]
plot(simfitepi2)
intensityplot(simfitepi2, which="total", aggregate=FALSE,
              col=rgb(0,0,0,alpha=if(getRversion() < "2.7.0") 1 else 0.1))

Run the code above in your browser using DataLab