seir: 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 = function (t) 0, mu = function (t) 0,
     sigma = 1, gamma = 1, delta = 0,
     m = 1L, n = 1L, init,
     stochastic = TRUE, prob = 1, delay = 1,
     aggregate = FALSE, useCompiled = TRUE, ...)
## A basic wrapper for the m=0L case:
 sir(length.out = 1L,
     beta, nu = function (t) 0, mu = function (t) 0,
     gamma = 1, delta = 0,
     n = 1L, init,
     stochastic = TRUE, prob = 1, delay = 1,
     aggregate = FALSE, useCompiled = TRUE, ...)

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.

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. m*sigma, n*gamma, and delta are the rates of removal from each latent, infectious, and recovered compartment.
m: a non-negative integer indicating a number of latent stages.
n: a positive integer indicating a number of infectious stages.
init: a numeric vector of length 1+m+n+1 giving an initial state with compartments ordered as $(S, E, I, R)$.
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.
aggregate: a logical indicating if latent and infectious compartments should be aggregated.
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.

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. tools:::Rd_expr_doi("10.1063/1.2745299")

Examples

Run this code

# \dontshow{
## for R_DEFAULT_PACKAGES=NULL
library(stats, pos = "package:base", verbose = FALSE)
library(utils, pos = "package:base", verbose = FALSE)
# }
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 = 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.ts, replicate(r, simplify = FALSE,
	seir(length.out, beta, nu, mu, sigma, gamma, delta, init = init,
	     prob = prob, delay = delay, epsilon = 0.002)[, "Z.obs"]))
str(Y)
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")

})

Run the code above in your browser using DataLab