fastbeta (version 0.5.1)

sir.aoi: 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(),
        root = NULL, root.max = 1L, root.break = TRUE,
        aggregate = FALSE, skip.Y = FALSE, ...)
# S3 method for sir.aoi
summary(object, name = "Y", tol = 1e-6, ...)

Value

A “multiple time series” object, inheriting from class sir.aoi and transitively from class mts, storing the numerical solution. Beneath the class, it is a T-by-(1+n+1) numeric matrix of the form cbind(S, I, Y), T <= length(seq(from, to, by)).

If root is a function, then an attribute root.info of the form list(tau, state = cbind(S, I, Y)) stores the first K

roots of that function and the state of the system at each root, K <= root.max.

If aggregate = TRUE, then infected compartments are aggregated so that the number of columns named I is 1 rather than n. This column stores prevalence, the proportion of the population that is infected. For convenience, there are 5 additional columns named I.E, I.I, foi, inc, and crv. These store the non-infectious and infectious components of prevalence (so that I.E + I.I = I), the force of infection, incidence (so that foi * S = inc), and the curvature of Y.

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.

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).

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. Roots of this function in the interval from from to to are sought alongside the numerical solution.

root.max

a positive integer giving a stopping condition for the root finder. Root finding continues until the count of roots found is root.max.

root.break

a logical indicating if the solver should stop once root.max roots are found. If TRUE, then the numerical solution ends at the last time point less than or equal to the last root.

aggregate

a logical indicating if infected compartments should be aggregated.

skip.Y

a logical indicating if solution of the equation for \(Y\) should not be attempted, e.g., because that equation seems stiffer than the rest.

...

optional arguments passed to the solver, function lsoda in package deSolve.

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). 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} F I_{j}) H(S) \\ \text{d} & I_{ 1} &{} / \text{d} t &{} = (\textstyle\sum_{j} \beta_{j} F I_{j}) 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}) F I_{j}) H(S) \\ \text{d} & I_{ 1} &{} / \text{d} \tau &{} = (\textstyle\sum_{j} (\mathcal{R}_{0,j}/\ell_{j}) F I_{j}) 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}) F I_{j}) (H(S) - 1/(\textstyle\sum_{j} \mathcal{R}_{0, j} F)) $$ due to the usefulness of the solution \(Y\) in applications.

Examples

# \dontshow{
## for R_DEFAULT_PACKAGES=NULL
library(   stats, pos = "package:base", verbose = FALSE)
library(graphics, pos = "package:base", verbose = FALSE)
# }
if (requireNamespace("deSolve")) withAutoprint({

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

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

info.peak <- attr(soln.peak, "root.info")
to <- 4 * info.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::lsoda

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 = info.peak[["tau"]], h = info.peak[["state"]][, "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, ...)
	attr(sir.aoi(to = to, by = by, R0 = R0, ell = ell,
	             root = root, aggregate = TRUE, ...),
	     "root.info")[["tau"]]
Ymax <- info.peak[["state"]][, "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 more than once
L <- list(function (Y, dY)             Y - Ymax * 0.5,
          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, root.max = 2L)

})