sir.aoi: Solve the SIR Equations Structured by Age of Infection

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

Solve the SIR Equations Structured by Age of Infection

Description

Numerically integrates the SIR equations with rates of transmission and recovery structured by age of infection.

Usage

sir.aoi(from = 0, to = from + 1, by = 1,
        R0, ell = 1, eps = 0, n = max(length(R0), length(ell)),
        init = c(1 - init.infected, init.infected),
        init.infected = .Machine[["double.neg.eps"]],
        weights = rep(c(1, 0), c(1L, n - 1L)),
        F = function (x) 1, Fargs = list(),
        H = identity, Hargs = list(),
        method = c("lsoda", "lsode", "vode", "daspk", "radau"),
        root = NULL, aggregate = FALSE, ...)

## S3 method for class 'sir.aoi'
summary(object, name = "Y", 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 mean time spent infectious.
R0	a numeric vector of length n such that sum(R0) is the basic reproduction number and R0[j] is the contribution of infected compartment j. Otherwise, a numeric vector of length 1, handled as equivalent to rep(R0/n, n).
ell	a numeric vector of length n such that ell[j] is the ratio of the mean time spent in infected compartment j and the mean time spent infectious; internally, ell/sum(ell[R0 > 0]) is used, hence ell is determined only up to a positive factor. Otherwise (and by default), a numeric vector of length 1, handled as equivalent to rep(1, n).
eps	a non-negative number giving the the ratio of the mean time spent infectious and the mean life expectancy; eps = 0 implies that life expectancy is infinite (that there are no deaths).
n	a positive integer giving the number of infected compartments. Setting n and thus overriding the default expression is necessary only if the lengths of R0 and ell are both 1.
init	a numeric vector of length 2 giving initial susceptible and infected proportions.
init.infected	a number in (0, 1] used only to define the default expression for init; see ‘Usage’.
weights	a numeric vector of length n containing non-negative weights, defining the initial distribution of infected individuals among the infected compartments. By default, all infected individuals occupy the first compartment.
F, H	functions returning a numeric vector of length 1 or of length equal that of the first formal argument. The body must be symbolically differentiable with respect to the first formal argument; see D.
Fargs, Hargs	lists of arguments passed to F or H. In typical usage, F and H define parametric families of functions of one variable, and Fargs and Hargs supply parameter values. For example: H = function(x, h) x^h, Hargs = list(h = 0.996).
method	a character string indicating a numerical integration method as the name of a solver in package deSolve. Otherwise, an object of class rkMethod specifying a Runge-Kutta method, obtained as the return value of a call to rkMethod.
root	a function returning a numeric vector of length 1, with formal arguments (tau, S, I, Y, dS, dI, dY, R0, ell) (or a subset); otherwise, NULL. Rootfinding is supported by three solvers: "lsoda", "lsode", and "radau".
aggregate	a logical indicating if infected compartments should be aggregated.
...	optional arguments passed to the solver indicated by method, used to set control parameters.
object	an R object inheriting from class sir.aoi, typically the value of a call to sir.aoi.
name	a character string in colnames(object), not "S". Tail exponents of V are approximated, where, for example, V = Y if name = "Y" and V = \sum_{j} I_{j} (prevalence) if name = "I".
tol	a positive number giving an upper bound on the relative change (from one time point to the next) in the slope of \log(V), defining time windows in which growth or decay of V is considered to be exponential.

Details

The SIR equations with rates of transmission and recovery structured by age of infection are

\begin{alignedat}{4} \text{d} & S &{} / \text{d} t &{} = \mu (1 - S) - (\textstyle\sum_{j} \beta_{j} I_{j}) F H(S) \\ \text{d} & I_{ 1} &{} / \text{d} t &{} = (\textstyle\sum_{j} \beta_{j} I_{j}) F H(S) - (\gamma_{1} + \mu) I_{1} \\ \text{d} & I_{j + 1} &{} / \text{d} t &{} = \gamma_{j} I_{j} - (\gamma_{j + 1} + \mu) I_{j + 1} \\ \text{d} & R &{} / \text{d} t &{} = \gamma_{n} I_{n} - \mu R \end{alignedat}

where S, I_{j}, R \ge 0, S + \sum_{j} I_{j} + R = 1, F is a forcing function, and H is a susceptible heterogeneity function. In general, F and H are nonlinear. In the standard SIR equations, F is 1 and H is the identity function.

Nondimensionalization using parameters \mathcal{R}_{0,j} = \beta_{j}/(\gamma_{j} + \mu), \ell_{j} = (1/(\gamma_{j} + \mu))/t_{+}, and \varepsilon = t_{+}/(1/\mu) and independent variable \tau = t/t_{+}, where t_{+} = \sum_{j:\mathcal{R}_{0,j} > 0} 1/(\gamma_{j} + \mu) designates as a natural time unit the mean time spent infectious, gives

\begin{alignedat}{4} \text{d} & S &{} / \text{d} \tau &{} = \varepsilon (1 - S) - (\textstyle\sum_{j} (\mathcal{R}_{0,j}/\ell_{j}) I_{j}) F H(S) \\ \text{d} & I_{ 1} &{} / \text{d} \tau &{} = (\textstyle\sum_{j} (\mathcal{R}_{0,j}/\ell_{j}) I_{j}) F H(S) - (1/\ell_{1} + \varepsilon) I_{1} \\ \text{d} & I_{j + 1} &{} / \text{d} \tau &{} = (1/\ell_{j}) I_{j} - (1/\ell_{j+1} + \varepsilon) I_{j + 1} \\ \text{d} & R &{} / \text{d} \tau &{} = (1/\ell_{n}) I_{n} - \varepsilon R \end{alignedat}

sir.aoi works with the nondimensional equations, dropping the last equation (which is redundant given R = 1 - S - \sum_{j} I_{j}) and augments the resulting system of 1 + n equations with a new equation

\text{d}Y/\text{d}\tau = (\textstyle\sum_{j} (\mathcal{R}_{0,j}/\ell_{j}) I_{j}) (H(S) - 1/(\textstyle\sum_{j} \mathcal{R}_{0, j}))

due to the usefulness of the solution Y in applications.

Value

root = NULL	a “multiple time series” object, inheriting from class sir.aoi and transitively from class mts. Beneath the class, it is a length(seq(from, to, by))-by-(1+n+1) numeric matrix of the form cbind(S, I, Y).
root = function (tau, S, I, Y, dS, dI, dY, R0, ell)	a numeric vector of length 1+1+n+1 of the form c(tau, S, I, Y) storing the root of the function root in units of the mean time spent infectious and the state at that time. Attribute curvature stores the curvature of Y at the root. If a root is not found between times from and to, then the value is NULL.

If aggregate = TRUE, then infected compartments are aggregated so that the number of columns (elements, if root is a function) named I is 1 rather than n. This column (element) stores prevalence, the proportion of the population that is infected. For convenience, there are 4 additional columns (elements) named I.E, I.I, foi, and inc. These store the non-infectious and infectious components of prevalence (so that I.E + I.I = I), the force of infection, and incidence.

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

Note

sir.aoi is not a special case of sir nor a generalization. The two functions were developed independently and for different purposes: sir.aoi to validate analytical results concerning the SIR equations as formulated here, sir to simulate incidence time series suitable for testing fastbeta.

Examples

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

to <- 100; by <- 0.01; R0 <- c(0, 16); ell <- c(0.5, 1)

peak <- sir.aoi(to = to, by = by, R0 = R0, ell = ell,
                root = function (S, R0) sum(R0) * S - 1,
                aggregate = TRUE)
peak

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

soln <- sir.aoi(to = to, by = by, R0 = R0, ell = ell,
                aggregate = TRUE, atol = 0, rtol = 1e-12)
##                                ^^^^      ^^^^
## 'atol', 'rtol', ... are passed to deSolve::ode

head(soln)
tail(soln)

plot(soln) # dispatching stats:::plot.ts
plot(soln, log = "y")

(lamb <- summary(soln)) # left and right tail exponents

xoff <- function (x, k) x - x[k]
k <- c(16L, nrow(soln)) # c(1L, nrow(soln)) worse due to transient

plot(soln[, "Y"], log = "y", ylab = "Y")
abline(v = peak[["tau"]], h = peak[["Y"]],
       lty = 2, lwd = 2, col = "red")
for (i in 1:2)
	lines(soln[k[i], "Y"] * exp(lamb[i] * xoff(time(soln), k[i])),
	      lty = 2, lwd = 2, col = "red")

wrap <-
function (root)
	sir.aoi(to = to, by = by, R0 = R0, ell = ell,
	        root = root, aggregate = TRUE)
Ymax <- peak[["Y"]]

## NB: want *simple* roots, not *multiple* roots
L <- list(function (Y) (Y - Ymax * 0.5)  ,
          function (Y) (Y - Ymax * 0.5)^2,
          function (Y) (Y - Ymax      )  ,
          function (Y) (Y - Ymax      )^2)
lapply(L, wrap)

## NB: critical values can be attained twice
L <- list(function (Y, dY) if (dY > 0) Y - Ymax * 0.5 else 1,
          function (Y, dY) if (dY < 0) Y - Ymax * 0.5 else 1)
lapply(L, wrap)

})

davidearn/fastbeta documentation built on July 4, 2025, 6:28 p.m.