seir: Simulate Infectious Disease Time Series

In davidearn/fastbeta: Fast Estimation of Time-Varying Infectious Disease Transmission Rates

Simulate Infectious Disease Time Series

Description

Simulates time series of the susceptible, infected, and removed population sizes and corresponding time series of births, incidence, and reported incidence. Simulations are based on an SIR model with user-defined forcing and a simple model for observation error.

Usage

sir(n, beta, nu, mu, constants, stochastic = TRUE,
    prob = 1, delay = 1, useCompiled = TRUE, ...)

Arguments

n	a positive integer. The simulation uses 0:n as time points, yielding n observation intervals of equal duration.
beta, nu, mu	functions of one or more argument returning transmission, birth, and natural death rates, respectively, at the time point indicated by the first argument. Arguments after the first must be strictly optional. The functions need not be vectorized.
constants	a numeric vector of the form c(S0, I0, R0, gamma, delta), specifying an initial state and rates of removal and loss of immunity, in that order.
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 reported. 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 reported during interval j+i-1, given that it is eventually reported. 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 ode (directly) or ssa.adaptivetau (via its list argument tl.params), depending on stochastic.

Details

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

\begin{aligned} \frac{\text{d} S}{\text{d} t} &= \nu(t) - \beta(t) S I + \delta R - \mu(t) S \\ \frac{\text{d} I}{\text{d} t} &= \beta(t) S I - \gamma I - \mu(t) I \\ \frac{\text{d} R}{\text{d} t} &= \gamma I - \delta R - \mu(t) R \end{aligned}

between times t = 0, \ldots, n, where t is understood to be a unitless measure of time relative to the duration of an observation interval. To generate time series of births and incidence

\begin{aligned} B(t) &= \int_{t - 1}^{t} \nu(s)\,\text{d}s \\ Z(t) &= \int_{t - 1}^{t} \beta(s) S(s) I(s)\,\text{d}s \end{aligned}

the system is augmented with two additional equations describing cumulative incidence and cumulative births (with right hand sides given by the integrands above), and the augmented system with five equations is integrated. Case reports are simulated by scaling incidence by prob and convolving the result with delay.

sir(stochastic = TRUE) works by simulating a Markov process corresponding to the augmented system, as described in the reference. Case reports are 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 an (n+1)-by-(5+r) numeric matrix X, where r is 0 if both prob and delay are missing (indicating no observation error) and 1 otherwise.

Rows correspond to time points 0:n. Columns are named S, I, R, B, Z, and (if r is 1) Z.obs. X[, 1:3] give the state at each time, and X[, 4:6] give the number of births, infections, and infections reported during the observation interval ending at each time.

X[1, 4:6] is NA, and X[2:length(delay), 6] can contain incomplete information if length(delay) >= 2.

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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1063/1.2745299")}

See Also

sir.library.

Examples

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

S0 <- 5e+04
I0 <- 1e+03
R0 <- 1e+06 - S0 - I0
constants <- c(S0 = S0, I0 = I0, R0 = R0, gamma = 0.5, delta = 0)

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

set.seed(0L)
X <- sir(n, beta, nu, mu, constants, prob = prob, delay = delay)
str(X)
plot(X)

r <- 10L
Y <- do.call(cbind, replicate(r, simplify = FALSE,
	sir(n, beta, nu, mu, constants, prob = prob, delay = delay)[, "Z.obs"]))
str(Y) # FIXME: stats:::cbind.ts mangles dimnames
plot(Y, plot.type = "single", col = seq_len(r), ylab = "Case reports")

davidearn/fastbeta documentation built on April 30, 2024, 2:35 a.m.