R: Simulate Infectious Disease Time Series

seir {fastbeta}

R Documentation

Simulate Infectious Disease Time Series

Description

Simulates incidence time series based on an SEIR model with user-defined forcing and a simple model for observation error.

Note that simulation code depends on availability of suggested packages adaptivetau and deSolve. If the dependency cannot be loaded then an error is signaled.

Usage

seir(length.out = 1L,
     beta, nu, mu, sigma = gamma, gamma = 1, delta = 0,
     init, m = length(init) - n - 2L, n = 1L,
     stochastic = TRUE, prob = 1, delay = 1,
     useCompiled = TRUE, ...)

## A basic wrapper for the 'm = 0L' case:

 sir(length.out = 1L,
     beta, nu, mu, gamma = 1, delta = 0,
     init, n = 1L,
     stochastic = TRUE, prob = 1, delay = 1,
     useCompiled = TRUE, ...)

Arguments

`length.out`	a non-negative integer indicating the time series length.
`beta`, `nu`, `mu`	functions of one or more arguments returning transmission, birth, and natural death rates at the time point indicated by the first argument. Arguments after the first must be strictly optional. The functions need not be vectorized.
`sigma`, `gamma`, `delta`	non-negative numbers. `msigma`, `ngamma`, and `delta` are the rates of removal from each latent, infectious, and recovered compartment.
`init`	a numeric vector of length `1+m+n+1` giving an initial state with compartments ordered as `(S, E, I, R)`.
`m`	a non-negative integer indicating a number of latent stages.
`n`	a positive integer indicating a number of infectious stages.
`stochastic`	a logical indicating if the simulation should be stochastic; see ‘Details’.
`prob`	a numeric vector of length `n` such that `prob[i]` is the probability that an infection during interval `i` is eventually observed. `prob` of length 1 is recycled.
`delay`	a numeric vector of positive length such that `delay[i]` is the probability that an infection during interval `j` is observed during interval `j+i-1`, given that it is eventually observed. `delay` need not sum to 1 but must not sum to 0.
`useCompiled`	a logical indicating if derivatives should be computed by compiled C functions rather than by R functions (which may be byte-compiled). Set to `FALSE` only if `TRUE` seems to cause problems, and in that case please report the problems with `bug.report(package = "fastbeta")`.
`...`	optional arguments passed to `lsoda` (directly) or `ssa.adaptivetau` (via its list argument `tl.params`), depending on `stochastic`.

Details

Simulations are based on an SEIR model with

m latent stages (E^{i}, i = 1,\ldots,m);
n infectious stages (I^{j}, j = 1,\ldots,n);
time-varying rates \beta, \nu, and \mu of transmission, birth, and natural death; and
constant rates m \sigma, n \gamma, and \delta of removal from each latent, infectious, and recovered compartment, where removal from the recovered compartment implies return to the susceptible compartment (loss of immunity).

seir(stochastic = FALSE) works by numerically integrating the system of ordinary differential equations

\begin{alignedat}{10} \text{d} & S &{} / \text{d} t &{} = {}& \delta &R &{} - ( && \lambda(t) &{} + \mu(t)) S &{} + \nu(t) \\ \text{d} & E^{ 1} &{} / \text{d} t &{} = {}& \lambda(t) &S &{} - ( && m \sigma &{} + \mu(t)) E^{ 1} &{} \\ \text{d} & E^{i + 1} &{} / \text{d} t &{} = {}& m \sigma &E^{i} &{} - ( && m \sigma &{} + \mu(t)) E^{i + 1} &{} \\ \text{d} & I^{ 1} &{} / \text{d} t &{} = {}& m \sigma &E^{m} &{} - ( && n \gamma &{} + \mu(t)) I^{ 1} &{} \\ \text{d} & I^{j + 1} &{} / \text{d} t &{} = {}& n \gamma &I^{j} &{} - ( && n \gamma &{} + \mu(t)) I^{j + 1} &{} \\ \text{d} & R &{} / \text{d} t &{} = {}& n \gamma &I^{n} &{} - ( && \delta &{} + \mu(t)) R &{} \end{alignedat} \\ \lambda(t) = \beta(t) \sum_{j} I^{j}

where it is understood that the independent variable t is a unitless measure of time relative to an observation interval. To get time series of incidence and births, the system is augmented with two equations describing cumulative incidence and births

\begin{aligned} \text{d} Z / \text{dt} &{} = \lambda(t) S \text{d} B / \text{dt} &{} = \nu(t) \end{aligned}

and the augmented system is numerically integrated. Observed incidence is simulated from incidence by scaling the latter by prob and convolving the result with delay.

seir(stochastic = TRUE) works by simulating a Markov process corresponding to the augmented system, as described in the reference. Observed incidence is simulated from incidence by binning binomial samples taken with probabilities prob over future observation intervals according to multinomial samples taken with probabilities delay.

Value

A “multiple time series” object, inheriting from class mts. Beneath the class, it is a length.out-by-1+m+n+1+2 numeric matrix with columns S, E, I, R, Z, and B, where Z and B specify incidence and births as the number of infections and births since the previous time point.

If prob or delay is not missing, then there is an additional column Z.obs specifying observed incidence as the number of infections observed since the previous time point. The first length(delay) elements of this column contain partial counts.

References

Cao, Y., Gillespie, D. T., & Petzold, L. R. (2007). Adaptive explicit-implicit tau-leaping method with automatic tau selection. Journal of Chemical Physics, 126(22), Article 224101, 1-9. doi:10.1063/1.2745299

Examples


if (requireNamespace("adaptivetau")) withAutoprint({

beta <- function (t, a = 1e-01, b = 1e-05) b * (1 + a * sinpi(t / 26))
nu   <- function (t) 1e+03
mu   <- function (t) 1e-03

sigma <- 0.5
gamma <- 0.5
delta <- 0

init <- c(S = 50200, E = 1895, I = 1892, R = 946011)

length.out <- 250L
prob <- 0.1
delay <- diff(pgamma(0:8, 2.5))

set.seed(0L)
X <- seir(length.out, beta, nu, mu, sigma, gamma, delta, init,
          prob = prob, delay = delay, epsilon = 0.002)
##                                    ^^^^^^^
## default epsilon = 0.05 allows too big leaps => spurious noise
##
str(X)
plot(X)

r <- 10L
Y <- do.call(cbind, replicate(r, simplify = FALSE,
	seir(length.out, beta, nu, mu, sigma, gamma, delta, init,
	     prob = prob, delay = delay, epsilon = 0.002)[, "Z.obs"]))
str(Y) # FIXME: stats:::cbind.ts mangles dimnames
plot(window(Y, start = tsp(Y)[1L] + length(delay) / tsp(Y)[3L]),
     ##        ^^^^^
     ## discards points showing edge effects due to 'delay'
     ##
     plot.type = "single", col = seq_len(r), ylab = "case reports")

})

[Package fastbeta version 0.3.0 Index]