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 = NULL)
## S3 method for class 'twinSIR':
simulate(object, nsim = 1, seed = 1,
         infPeriod = NULL, remPeriod = NULL,
         end = diff(range(object$intervals)), trace = FALSE, .allocate = NULL,
         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

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

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 (i.e.
    .allocate*N rows). By default (NULL), this number is set to
    max(500, ceiling(nBlocks/100)*100), i.e. 500 but at least the
    number of block

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.
If nsim > 1 epidemics are simulated by the
  simulate-method for fitted "twinSIR" models, these are
  returned in a list.

`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## Generate a data frame containing a hypothetic population with 100 individuals
set.seed(1234)
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
set.seed(1)
epi <- simEpidata(cbind(start,stop) ~ cox(gender), 
                  data = pop, 
                  id = "id", I0.col = "I0col", coords.cols = c("x","y"),
                  beta = c(-2), h0 = -1, 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)


if (surveillance.options("allExamples"))
{
  ## load data of an artificial epidemic
  data("fooepidata")
  summary(fooepidata)

  ## simulate a new evolution of the epidemic
  set.seed(1)
  simepi <- simEpidata(cbind(start, stop) ~ cox(z1) + cox(z1):cox(z2),
      data = fooepidata, 
      beta = c(1,0.5), h0 = -7, alpha = c(B2 = 0.01, B1 = 0.005),
      f = list(B1 = function(u) u<=1, B2 = function(u) u>1 & is.finite(u)),
      infPeriod = function(ids) rexp(length(ids), rate=0.2), trace = FALSE)
  summary(simepi)
  plot(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
# a) with observed infPeriods (i.e. fixed length 3 days):
simfitepi1 <- simulate(foofit)
plot(simfitepi1)
# b) with new infPeriods (simuluated from the Exp(0.3) distribution):
simfitepi2 <- simulate(foofit,
                       infPeriod=function(ids) rexp(length(ids), rate=0.3))
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