sir.aoi {fastbeta}R Documentation

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, 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)),
        root = NULL, aggregate = FALSE, ...)

## S3 method for class 'sir.aoi'
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 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).

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.

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.

aggregate

a logical indicating if infected compartments should be aggregated.

...

optional arguments passed to lsoda.

object

an R object inheriting from class sir.aoi, typically the value of a call to sir.aoi.

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.

Details

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

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

where N = S + \sum_{j} I_{j} + R is the (constant, positive) population size. Nondimensionalization using parameters N = 1, \mathcal{R}_{0,j} = \beta_{j}/\gamma_{j}, and \ell_{j} = (1/\gamma_{j})/\sum_{j:\mathcal{R}_{0,j} > 0} (1/\gamma_{j}) and time unit \tau = t/\sum_{j:\mathcal{R}_{0,j} > 0} (1/\gamma_{j}), gives

\begin{alignedat}{4} \text{d} & S &{} / \text{d} \tau &{} = -\textstyle\sum_{j} (\mathcal{R}_{0,j}/\ell_{j}) S I_{j} \\ \text{d} & I_{ 1} &{} / \text{d} \tau &{} = \textstyle\sum_{j} (\mathcal{R}_{0,j}/\ell_{j}) S I_{j} - (1/\ell_{1}) I_{1} \\ \text{d} & I_{j + 1} &{} / \text{d} \tau &{} = (1/\ell_{j}) I_{j} - (1/\ell_{j+1}) I_{j + 1} \\ \text{d} & R &{} / \text{d} \tau &{} = (1/\ell_{n}) I_{n} \\ \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 = (\sum_{j} \mathcal{R}_{0, j} S - 1) \sum_{j:\mathcal{R}_{0,j} > 0} I_{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 or element stores prevalence, the proportion of the population that is infected. For convenience, there are 2 additional columns (elements) named I.E and I.I. These store the non-infectious and infectious components of prevalence, as indicated by sign(R0), hence I.E + I.I = I.

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.

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)
head(soln)

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

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

xoff <- function (x, k) x - x[k]
lamb <- summary(soln)
k <- c(16L, nrow(soln)) # c(1L, nrow(soln)) suffers due to transient
plot(soln[, "I"], log = "y", ylab = "Prevalence")
for (i in 1:2)
	lines(soln[k[i], "I"] * 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
F <- list(function (Y) (Y - Ymax * 0.5)  ,
          function (Y) (Y - Ymax * 0.5)^2,
          function (Y) (Y - Ymax      )  ,
          function (Y) (Y - Ymax      )^2)
lapply(F, wrap)

## NB: critical values can be attained twice
F <- 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(F, wrap)

})

[Package fastbeta version 0.4.0 Index]