surveillance: twinstim_simulation – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

twinstim_simulation

Simulation of a Self-Exciting Spatio-Temporal Point Process

Description

The function simEpidataCS simulates events of a self-exciting spatio-temporal point process of the "twinstim" class. Simulation works via Ogata's modified thinning of the conditional intensity as described in Meyer et al. (2012). Note that simulation is limited to the spatial and temporal range of stgrid.

The simulate method for objects of class "twinstim" simulates new epidemic data using the model and the parameter estimates of the fitted object.

Usage

simEpidataCS(endemic, epidemic, siaf, tiaf, qmatrix, rmarks,
    events, stgrid, tiles, beta0, beta, gamma, siafpars, tiafpars,
    epilink = "log", t0 = stgrid$start[1], T = tail(stgrid$stop,1),
    nEvents = 1e5, control.siaf = list(F=list(), Deriv=list()),
    W = NULL, trace = 5, nCircle2Poly = 32, gmax = NULL, .allocate = 500,
    .skipChecks = FALSE, .onlyEvents = FALSE)

## S3 method for class 'twinstim'
simulate(object, nsim = 1, seed = NULL, data, tiles,
    newcoef = NULL, rmarks = NULL, t0 = NULL, T = NULL, nEvents = 1e5,
    control.siaf = object$control.siaf,
    W = data$W, trace = FALSE, nCircle2Poly = NULL, gmax = NULL,
    .allocate = 500, simplify = TRUE, ...)

Arguments

`endemic`	see `twinstim`. Note that type-specific endemic intercepts are specified by `beta0` here, not by the term `(1\|type)`.
`epidemic`	see `twinstim`. Marks appearing in this formula must be returned by the generating function `rmarks`.
`siaf`	see `twinstim`. In addition to what is required for fitting with `twinstim`, the `siaf` specification must also contain the element `simulate`, a function which draws random locations following the spatial kernel `siaf$f`. The first argument of the function is the number of points to sample (say `n`), the second one is the vector of parameters `siafpars`, the third one is the type indicator (a character string matching a type name as specified by `dimnames(qmatrix)`). With the current implementation there will always be simulated only one location at a time, i.e. `n=1`. The predefined siaf's all provide simulation.
`tiaf`	e.g. what is returned by the generating function `tiaf.constant` or `tiaf.exponential`. See also `twinstim`.
`qmatrix`	see `epidataCS`. Note that this square matrix and its `dimnames` determine the number and names of the different event types. In the simplest case, there is only a single type of event, i.e. `qmatrix = diag(1)`.
`rmarks`	function of single time (1st argument) and location (2nd argument) returning a one-row `data.frame` of marks (named according to the variables in `epidemic`) for an event at this point. This must include the columns `eps.s` and `eps.t`, i.e. the values of the spatial and temporal interaction ranges at this point. Only `"numeric"` and `"factor"` columns are allowed. Assure that factor variables are coded equally (same levels and level order) for each new sample. For the `simulate.twinstim` method, the default (`NULL`) means sampling from the empirical distribution function of the (non-missing) marks in `data` restricted to events in the simulation period (`t0`;`T`]. If there are no events in this period, e.g., if simulating beyond the original observation period, `rmarks` will sample marks from all of `data$events`.
`events`	`NULL` or missing (default) in case of an empty prehistory, or a `SpatialPointsDataFrame` containing events of the prehistory (-Inf;`t0`] of the process (required for the epidemic to start in case of no endemic component in the model). The `SpatialPointsDataFrame` must have the same `proj4string` as `tiles` and `W`). The attached `data.frame` (data slot) must contain the typical columns as described in `as.epidataCS` (`time`, `tile`, `eps.t`, `eps.s`, and, for type-specific models, `type`) and all marks appearing in the `epidemic` specification. Note that some column names are reserved (see `as.epidataCS`). Only events up to time `t0` are selected and taken as the prehistory.
`stgrid`	see `as.epidataCS`. Simulation only works inside the spatial and temporal range of `stgrid`.
`tiles`	object inheriting from `"SpatialPolygons"` with `row.names` matching the `tile` names in `stgrid` and having the same `proj4string` as `events` and `W`. This is necessary to sample the spatial location of events generated by the endemic component.
`beta0,beta,gamma,siafpars,tiafpars`	these are the parameter subvectors of the `twinstim`. `beta` and `gamma` must be given in the same order as they appear in `endemic` and `epidemic`, respectively. `beta0` is either a single endemic intercept or a vector of type-specific endemic intercepts in the same order as in `qmatrix`.
`epilink`	a character string determining the link function to be used for the `epidemic` linear predictor of event marks. By default, the log-link is used. The experimental alternative is `epilink = "identity"`. Note that the identity link does not guarantee the force of infection to be positive. If this leads to a negative total intensity (endemic + epidemic), the point process is not well defined and simulation cannot proceed.
`t0`	`events` having occurred during (-Inf;`t0`] are regarded as part of the prehistory H_0 of the process. For `simEpidataCS`, by default and also if `t0=NULL`, the beginning of `stgrid` is used as `t0`. For the `simulate.twinstim` method, `NULL` means to use the fitted time range of the `"twinstim"` `object`.
`T, nEvents`	simulate a maximum of `nEvents` events up to time `T`, then stop. For `simEpidataCS`, by default, and also if `T=NULL`, `T` equals the last stop time in `stgrid` (it cannot be greater) and `nEvents` is bounded above by 10000. For the `simulate.twinstim` method, `T=NULL` means to use the same same time range as for the fitting of the `"twinstim"` `object`.
`W`	see `as.epidataCS`. When simulating from `twinstim`-fits, `W` is by default taken from the original `data$W`. If specified as `NULL`, `W` is generated automatically via `unionSpatialPolygons(tiles)`. However, since the result of such a polygon operation should always be verified, it is recommended to do that in advance. It is important that `W` and `tiles` cover the same region: on the one hand direct offspring is sampled in the spatial influence region of the parent event, i.e., in the intersection of `W` and a circle of radius the `eps.s` of the parent event, after which the corresponding tile is determined by overlay with `tiles`. On the other hand endemic events are sampled from `tiles`.
`trace`	logical (or integer) indicating if (or how often) the current simulation status should be `cat`ed. For the `simulate.twinstim` method, `trace` currently only applies to the first of the `nsim` simulations.
`.allocate`	number of rows (events) to initially allocate for the event history; defaults to 500. Each time the simulated epidemic exceeds the allocated space, the event `data.frame` will be enlarged by `.allocate` rows.
`.skipChecks,.onlyEvents`	these logical arguments are not meant to be set by the user. They are used by the `simulate`-method for `"twinstim"` objects.
`object`	an object of class `"twinstim"`.
`nsim`	number of epidemics (i.e. spatio-temporal point patterns inheriting from class `"epidataCS"`) to simulate. Defaults to 1 when the result is a simple object inheriting from class `"simEpidataCS"` (as if `simEpidataCS` would have been called directly). If `nsim > 1`, the result will be a list the structure of which depends on the argument `simplify`.
`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.
`data`	an object of class `"epidataCS"`, usually the one to which the `"twinstim"` `object` was fitted. It carries the `stgrid` of the endemic component, but also `events` for use as the prehistory, and defaults for `rmarks` and `nCircle2Poly`.
`newcoef`	an optional named numeric vector of (a subset of) parameters to replace the original point estimates in `coef(object)`. Elements which do not match any model parameter by name are silently ignored. The `newcoef`s may also be supplied in a list following the same conventions as for the `start` argument in `twinstim`.
`simplify`	logical. It is strongly recommended to set `simplify = TRUE` (default) if `nsim` is large. This saves space and computation time, because for each simulated epidemic only the `events` component is saved. All other components, which do not vary between simulations, are only stored from the first run. In this case, the runtime of each simulation is stored as an attribute `"runtime"` to each simulated `events`. See also the “Value” section below.
`control.siaf`	see `twinstim`.
`nCircle2Poly`	see `as.epidataCS`. For `simulate.twinstim`, `NULL` means to use the same value as for `data`.
`gmax`	maximum value the temporal interaction function `tiaf$g` can attain. If `NULL`, then it is assumed as the maximum value of the type-specific values at 0, i.e. `max(tiaf$g(rep.int(0,nTypes), tiafpars, 1:nTypes))`.
`...`	unused (arguments of the generic).

Value

The function simEpidataCS returns a simulated epidemic of class "simEpidataCS", which enhances the class "epidataCS" by the following additional components known from objects of class "twinstim": bbox, timeRange, formula, coefficients, npars, control.siaf, call, runtime. It has corresponding coeflist, residuals, R0, and intensityplot methods.

The simulate.twinstim method has some additional attributes set on its result: call, seed, and runtime. If nsim > 1, it returns an object of class "simEpidataCSlist", the form of which depends on the value of simplify (which is stored as an attribute simplified): if simplify = FALSE, then the return value is just a list of sequential simulations, each of class "simEpidataCS". However, if simplify = TRUE, then the sequential simulations share all components but the simulated events, i.e. the result is a list with the same components as a single object of class "simEpidataCS", but with events replaced by an eventsList containing the events returned by each of the simulations.

The stgrid component of the returned "simEpidataCS" will be truncated to the actual end of the simulation, which might be <T, if the upper bound nEvents is reached during simulation.

CAVE: Currently, simplify=TRUE in simulate.twinstim ignores that multiple simulated epidemics (nsim > 1) may have different stgrid time ranges. In a "simEpidataCSlist", the stgrid shared by all of the simulated epidemics is just the stgrid returned by the first simulation.

Note

The more detailed the polygons in tiles are the slower is the algorithm. You are advised to sacrifice some shape details for speed by reducing the polygon complexity, for example via the mapshaper JavaScript library wrapped by the R package rmapshaper. Alternative tools are provided by the packages maptools (thinnedSpatialPoly) and spatstat.geom (simplify.owin).

Author(s)

Sebastian Meyer, with contributions by Michael Höhle

References

Douglas, D. H. and Peucker, T. K. (1973): Algorithms for the reduction of the number of points required to represent a digitized line or its caricature. Cartographica: The International Journal for Geographic Information and Geovisualization, 10, 112-122

Meyer, S., Elias, J. and Höhle, M. (2012): A space-time conditional intensity model for invasive meningococcal disease occurrence. Biometrics, 68, 607-616. doi: 10.1111/j.1541-0420.2011.01684.x

Examples

data("imdepi", "imdepifit")

## load borders of Germany's districts (originally obtained from
## the German Federal Agency for Cartography and Geodesy,
## https://gdz.bkg.bund.de/), simplified by the "modified Visvalingam"
## algorithm (level=6.6%) using MapShaper.org (v. 0.1.17):
load(system.file("shapes", "districtsD.RData", package="surveillance"))
plot(districtsD)
plot(stateD, add=TRUE, border=2, lwd=2)
# 'stateD' was obtained as 'rgeos::gUnaryUnion(districtsD)'

## simulate 2 realizations (over a short period, for speed)
## considering events from data(imdepi) before t=31 as prehistory
mysims <- simulate(imdepifit, nsim=2, seed=1, data=imdepi,
                   tiles=districtsD, newcoef=c("e.typeC"=-1),
                   t0=31, T=if (interactive()) 180 else 45, # for CRAN
                   simplify=TRUE)


## plot both simulations using the plot-method for simEpidataCSlist's
mysims
plot(mysims, aggregate="time")

## extract the second realization -> object of class simEpidataCS
mysim2 <- mysims[[2]]
summary(mysim2)
plot(mysim2, aggregate="space")

if (surveillance.options("allExamples")) {
### compare the observed _cumulative_ number of cases during the
### first 90 days to 20 simulations from the fitted model
### (performing these simulations takes about 30 seconds)

sims <- simulate(imdepifit, nsim=20, seed=1, data=imdepi, t0=0, T=90,
                 tiles=districtsD, simplify=TRUE)

## extract cusums
getcsums <- function (events) {
    tapply(events$time, events@data["type"],
           function (t) cumsum(table(t)), simplify=FALSE)
}
csums_observed <- getcsums(imdepi$events)
csums_simulated <- lapply(sims$eventsList, getcsums)

## plot it
plotcsums <- function (csums, ...) {
    mapply(function (csum, ...) lines(as.numeric(names(csum)), csum, ...),
           csums, ...)
    invisible()
}
plot(c(0,90), c(0,35), type="n", xlab="Time [days]",
     ylab="Cumulative number of cases")
plotcsums(csums_observed, col=c(2,4), lwd=3)
legend("topleft", legend=levels(imdepi$events$type), col=c(2,4), lwd=1)
invisible(lapply(csums_simulated, plotcsums,
                 col=adjustcolor(c(2,4), alpha=0.5)))
}


## Not run: 

### Experimental code to generate 'nsim' simulations of 'nm2add' months
### beyond the observed time period:
nm2add <- 24
nsim <- 5
### With these settings, simulations will take about 30 seconds.
### The events still infective by the end of imdepi$stgrid will be used
### as the prehistory for the continued process.

origT <- tail(imdepi$stgrid$stop, 1)
## create a time-extended version of imdepi
imdepiext <- local({
    ## first we have to expand stgrid (assuming constant "popdensity")
    g <- imdepi$stgrid
    g$stop <- g$BLOCK <- NULL
    gadd <- data.frame(start=rep(seq(origT, by=30, length.out=nm2add),
                                 each=nlevels(g$tile)),
                       g[rep(seq_len(nlevels(g$tile)), nm2add), -1])
    ## now create an "epidataCS" using this time-extended stgrid
    as.epidataCS(events=imdepi$events,  # the replacement warnings are ok
                 W=imdepi$W, qmatrix=imdepi$qmatrix,
                 stgrid=rbind(g, gadd), T=max(gadd$start) + 30)
})
newT <- tail(imdepiext$stgrid$stop, 1)

## simulate beyond the original period
simsext <- simulate(imdepifit, nsim=nsim, seed=1, t0=origT, T=newT,
                    data=imdepiext, tiles=districtsD, simplify=TRUE)

## Aside to understand the note from checking events and tiles:
# marks(imdepi)["636",]  # tile 09662 is attributed to this event, but:
# plot(districtsD[c("09678","09662"),], border=1:2, lwd=2, axes=TRUE)
# points(imdepi$events["636",])
## this mismatch is due to polygon simplification

## plot the observed and simulated event numbers over time
plot(imdepiext, breaks=c(unique(imdepi$stgrid$start),origT),
     cumulative=list(maxat=330))
for (i in seq_along(simsext$eventsList))
    plot(simsext[[i]], add=TRUE, legend.types=FALSE,
         breaks=c(unique(simsext$stgrid$start),newT),
         subset=!is.na(source),  # have to exclude the events of the prehistory
         cumulative=list(offset=c(table(imdepi$events$type)), maxat=330, axis=FALSE),
         border=NA, density=0)  # no histogram
abline(v=origT, lty=2, lwd=2)


## End(Not run)

surveillance

Temporal and Spatio-Temporal Modeling and Monitoring of Epidemic Phenomena

v1.19.1

GPL-2

Authors

Michael H<f6>hle [aut, ths] (<https://orcid.org/0000-0002-0423-6702>), Sebastian Meyer [aut, cre] (<https://orcid.org/0000-0002-1791-9449>), Michaela Paul [aut], Leonhard Held [ctb, ths], Howard Burkom [ctb], Thais Correa [ctb], Mathias Hofmann [ctb], Christian Lang [ctb], Juliane Manitz [ctb], Andrea Riebler [ctb], Daniel Saban<e9>s Bov<e9> [ctb], Ma<eb>lle Salmon [ctb], Dirk Schumacher [ctb], Stefan Steiner [ctb], Mikko Virtanen [ctb], Wei Wei [ctb], Valentin Wimmer [ctb], R Core Team [ctb] (A few code segments are modified versions of code from base R)

Initial release

2021-03-30