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 |
R0 |
a numeric vector of length |
ell |
a numeric vector of length |
n |
a positive integer giving the number of infected compartments.
Setting |
init |
a numeric vector of length 2 giving initial susceptible and infected proportions. |
init.infected |
a number in |
weights |
a numeric vector of length |
root |
a function returning a numeric vector of length 1, with formal
arguments |
aggregate |
a logical indicating if infected compartments should be aggregated. |
... |
optional arguments passed to |
object |
an R object inheriting from class |
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
|
root = function (tau , S , I , Y , dS , dI , dY , R0 , ell) |
a numeric vector of length |
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)
})