seir | R Documentation |
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.
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, ...)
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 |
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 |
stochastic |
a logical indicating if the simulation should be stochastic; see ‘Details’. |
prob |
a numeric vector of length |
delay |
a numeric vector of positive length such that |
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 |
... |
optional arguments passed to |
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
.
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.
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")}
seir.auxiliary
, seir.library
.
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")
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.