seir.canonical function

Solve the Canonical SEIR Equations

Numerically integrates the canonical SEIR equations, a special case of the more general SEIR equations handled by seir; see Details .

seir.canonical(from = 0, to = from + 1, by = 1,
               R0, ell = (2 * n)/(3 * n + 1),
               m = 1L, n = 1L,
               init = c(1 - p, p), p = .Machine[["double.neg.eps"]],
               weights = rep(c(1, 0), c(1L, m + n - 1L)),
               root = c("none", "peak"), aggregate = FALSE, ...)

## S3 method for class 'seir.canonical'
summary(object, tol = 1e-6, ...)

Arguments

• from, to, by: passed to seq.int in order to generate an increasing, equally spaced vector of time points in units of the generation interval.
• R0: a positive number giving the basic reproduction number.
• ell: a number in $(0,1)$ giving the ratio of the mean latent period and mean generation interval when m is positive. The default value implies that the mean latent period and mean infectious period are equal.
• 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 2 giving initial susceptible and infected proportions.
• p: a number in $(0,1]$ used only to define init when init is unset.
• weights: a numeric vector of length m+n containing non-negative weights, defining the initial distribution of infected individuals among the latent and infectious stages. By default, all infected individuals occupy the first stage.
• root: a character string determining what is returned. "none": the numerical solution as a multiple time series object. "peak": information about the time and state when $Y$ attains a local maximum.
• aggregate: a logical indicating if latent and infectious compartments should be aggregated when root = "none".
• ...: optional arguments passed to lsoda.
• object: an object inheriting from class seir.canonical, typically the value of a call to seir.canonical.
• tol: a positive number giving an upper bound on the relative change (from one time point to the next) in the slope of log prevalence, defining time windows in which growth or decay of prevalence is considered to be exponential.

Returns

If root = "none", a multiple time series object, inheriting from class seir.canonical and transitively from class mts. Beneath the class, it is a length(seq(from, to, by))-by-(1+m+n+1) numeric matrix with columns S, E, I, and Y.

If root = "peak", a numeric vector of length 5 of the form c(tau, S, E, I, Y)

containing the time in units of the generation interval at which $Y$

attains a local maximum and the values at that time of $S$, $sum(E)$, $sum(I)$, and $Y$. Attributes E.full, I.full, and curvature store the corresponding values of $E[i]$, $I[j]$, and $Y''$. If $Y$ does not attain a local maximum between times from and to, then all of the elements are NaN and there are no attributes. If m is zero, then the third element is NaN and there is no attribute E.full.

The method for summary returns a numeric vector of length 2 containing the left and right tail exponents , defined as the asymptotic values of the slope of log prevalence. NaN elements indicate that a tail exponent cannot be approximated from the prevalence time series represented by object, because the time window does not cover enough of the tail, where the meaning of enough is set by tol.

Details

The general SEIR equations handled by seir are

The canonical SEIR equations are obtained by substituting $beta(t) = R0 * gamma$, $\nu(t) = \mu(t) = \delta = 0$, $1/sigma = ell * g$, and $1/gamma = (1 - ell) * g/h$, where $h = (n + 1)/(2 * n)$, then choosing $\tau = t/g$ as an independent variable:

The constraint $\mu(t) = \nu(t)$ implies that the population size $N = S + \sum_{i} E^{i} + \sum_{j} I^{j} + R$ is constant, hence, without loss of generality, seir.canonical assumes $N = 1$

and drops the last equation. The resulting system of $1 + m + n$

equations is augmented with a final equation

due to the usefulness of the solution $Y$ in applications.

See Also

seir.

Examples

if (requireNamespace("deSolve")) withAutoprint({

to <- 100; by <- 0.01; R0 <- 16; ell <- 1/3

peak <- seir.canonical(to = to, by = by, R0 = R0, ell = ell,
                       root = "peak")
peak

to <- 4 * peak[["tau"]] # a more principled endpoint

soln <- seir.canonical(to = to, by = by, R0 = R0, ell = ell)
head(soln)

plot(soln) # dispatching stats:::plot.ts

plot(soln[, "Y"])
abline(v = peak[["tau"]], h = peak[["Y"]],
       lty = 2, lwd = 2, col = "red")

xoff <- function (x, k) x - x[k]
prev <- soln[, "E"] + soln[, "I"]
lamb <- summary(soln)
k <- c(16L, nrow(soln)) # c(1L, nrow(soln)) suffers due to transient
plot(prev, log = "y")
for (i in 1:2)
        lines(prev[k[i]] * exp(lamb[i] * xoff(time(prev), k[i])),
              lty = 2, lwd = 2, col = "red")

})