Version:	0.4.0
Date:	2025-04-28
Title:	Fast Approximation of Time-Varying Infectious Disease Transmission Rates
Description:	A fast method for approximating time-varying infectious disease transmission rates from disease incidence time series and other data, based on a discrete time approximation of an SEIR model, as analyzed in Jagan et al. (2020) <doi:10.1371/journal.pcbi.1008124>.
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
URL:	https://github.com/davidearn/fastbeta
BugReports:	https://github.com/davidearn/fastbeta/issues
Depends:	R (≥ 4.3)
Imports:	grDevices, graphics, stats
Suggests:	adaptivetau, deSolve, tools, utils
BuildResaveData:	no
NeedsCompilation:	yes
Packaged:	2025-04-28 15:31:28 UTC; mikael
Author:	Mikael Jagan [aut, cre]
Maintainer:	Mikael Jagan <jaganmn@mcmaster.ca>
Repository:	CRAN
Date/Publication:	2025-04-28 16:30:05 UTC

R Package fastbeta

Description

An R package for approximating time-varying infectious disease transmission rates from disease incidence time series and other data.

Details

The “main” function is fastbeta.

To render a list of available help topics, use help(package = "fastbeta").

To report a bug or request a change, use bug.report(package = "fastbeta").

Author(s)

Mikael Jagan jaganmn@mcmaster.ca

Richardson-Lucy Deconvolution

Description

Performs a modified Richardson-Lucy iteration for the purpose of estimating incidence from reported incidence or mortality, conditional on a reporting probability and on a distribution of the time to reporting.

Usage

deconvolve(x, prob = 1, delay = 1,
           start, tol = 1, iter.max = 32L, complete = FALSE)

Arguments

Value

A list with elements:

References

Goldstein, E., Dushoff, J., Ma, J., Plotkin, J. B., Earn, D. J. D., & Lipsitch, M. (2020). Reconstructing influenza incidence by deconvolution of daily mortality time series. Proceedings of the National Academy of Sciences U. S. A., 106(51), 21825-21829. doi:10.1073/pnas.0902958106

Examples

set.seed(2L)
n <- 200L
d <- 50L
p <- 0.1
prob <- plogis(rlogis(d + n, location = qlogis(p), scale = 0.1))
delay <- diff(pgamma(0L:(d + 1L), 12, 0.4))

h <- function (x, a = 1, b = 1, c = 0) a * exp(-b * (x - c)^2)
ans <- floor(h(seq(-60, 60, length.out = d + n), a = 1000, b = 0.001))

x0 <- rbinom(d + n, ans, prob)
x <- tabulate(rep(1L:(d + n), x0) +
              sample(0L:d, size = sum(x0), replace = TRUE, prob = delay),
              d + n)[-(1L:d)]

str(D0 <- deconvolve(x, prob, delay, complete = FALSE))
str(D1 <- deconvolve(x, prob, delay, complete =  TRUE))

matplot(-(d - 1L):n,
        cbind(x0, c(rep(NA, d), x), prob * D0[["value"]], p * ans),
        type = c("p", "p", "p", "l"),
        col = c(1L, 1L, 2L, 4L), pch = c(16L, 1L, 16L, NA),
        lty = c(0L, 0L, 0L, 1L), lwd = c(NA, NA, NA, 3L),
        xlab = "Time", ylab = "Count")
legend("topleft", NULL,
       c("actual", "actual+delay", "actual+delay+deconvolution", "p*h"),
       col = c(1L, 1L, 2L, 4L), pch = c(16L, 1L, 16L, NA),
       lty = c(0L, 0L, 0L, 1L), lwd = c(NA, NA, NA, 3L),
       bty = "n")

plot(0L:D1[["iter"]], D1[["chisq"]], xlab = "Iterations", ylab = quote(chi^2))
abline(h = 1, lty = 2L)

Estimate a Time-Varying Infectious Disease Transmission Rate

Description

Generates a discrete approximation of a time-varying infectious disease transmission rate from an equally spaced disease incidence time series and other data.

Usage

fastbeta(series, sigma = 1, gamma = 1, delta = 0,
         m = 1L, n = 1L, init, ...)

Arguments

Details

The algorithm implemented by fastbeta is based on an SEIR model with

• m latent stages (E^{i}, i = 1,\ldots,m);

• n infectious stages (I^{j}, j = 1,\ldots,n);

• time-varying rates \beta, \nu, and \mu of transmission, birth, and natural death; and

• constant rates m \sigma, n \gamma, and \delta of removal from each latent, infectious, and recovered compartment, where removal from the recovered compartment implies return to the susceptible compartment (loss of immunity).

It is derived by linearizing of the system of ordinary differential equations

\begin{alignedat}{10} \text{d} & S &{} / \text{d} t &{} = {}& \delta &R &{} - ( && \lambda(t) &{} + \mu(t)) S &{} + \nu(t) \\ \text{d} & E^{ 1} &{} / \text{d} t &{} = {}& \lambda(t) &S &{} - ( && m \sigma &{} + \mu(t)) E^{ 1} &{} \\ \text{d} & E^{i + 1} &{} / \text{d} t &{} = {}& m \sigma &E^{i} &{} - ( && m \sigma &{} + \mu(t)) E^{i + 1} &{} \\ \text{d} & I^{ 1} &{} / \text{d} t &{} = {}& m \sigma &E^{m} &{} - ( && n \gamma &{} + \mu(t)) I^{ 1} &{} \\ \text{d} & I^{j + 1} &{} / \text{d} t &{} = {}& n \gamma &I^{j} &{} - ( && n \gamma &{} + \mu(t)) I^{j + 1} &{} \\ \text{d} & R &{} / \text{d} t &{} = {}& n \gamma &I^{n} &{} - ( && \delta &{} + \mu(t)) R &{} \end{alignedat} \\ \lambda(t) = \beta(t) \sum_{j} I^{j}

and substituting actual or estimated incidence and births for definite integrals of \lambda S and \nu. This procedure yields a system of linear difference equations from which one recovers a discrete approximation of \beta:

\begin{alignedat}{17} &E_{t + 1}^{ 1} &{} = {}& [(1 - \tfrac{1}{2} ( & m \sigma + \mu_{t})) & E_{t}^{ 1} & & & & & & & &{} + Z_{t + 1} & ] &{} / [1 + \tfrac{1}{2} ( & m \sigma + \mu_{t + 1})] \\ &E_{t + 1}^{i + 1} &{} = {}& [(1 - \tfrac{1}{2} ( & m \sigma + \mu_{t})) & E_{t}^{i + 1} &{} + {}& \tfrac{1}{2} & m \sigma ( & E_{t}^{i} &{} + {}& E_{t + 1}^{i} & ) & & ] &{} / [1 + \tfrac{1}{2} ( & m \sigma + \mu_{t + 1})] \\ &I_{t + 1}^{ 1} &{} = {}& [(1 - \tfrac{1}{2} ( & n \gamma + \mu_{t})) & I_{t}^{ 1} &{} + {}& \tfrac{1}{2} & m \sigma ( & E_{t}^{m} &{} + {}& E_{t + 1}^{m} & ) & & ] &{} / [1 + \tfrac{1}{2} ( & n \gamma + \mu_{t + 1})] \\ &I_{t + 1}^{j + 1} &{} = {}& [(1 - \tfrac{1}{2} ( & n \gamma + \mu_{t})) & I_{t}^{j + 1} &{} + {}& \tfrac{1}{2} & n \gamma ( & I_{t}^{j} &{} + {}& I_{t + 1}^{j} & ) & & ] &{} / [1 + \tfrac{1}{2} ( & n \gamma + \mu_{t + 1})] \\ &R_{t + 1} &{} = {}& [(1 - \tfrac{1}{2} ( & \delta + \mu_{t})) & R_{t} &{} + {}& \tfrac{1}{2} & n \gamma ( & I_{t}^{n} &{} + {}& I_{t + 1}^{n} & ) & & ] &{} / [1 + \tfrac{1}{2} ( & \delta + \mu_{t + 1})] \\ &S_{t + 1} &{} = {}& [(1 - \tfrac{1}{2} ( & \mu_{t})) & S_{t} &{} + {}& \tfrac{1}{2} & \delta ( & R_{t} &{} + {}& R_{t + 1} & ) &{} - Z_{t + 1} &{} + B_{t + 1}] &{} / [1 + \tfrac{1}{2} ( & \mu_{t + 1})] \end{alignedat} \\ \beta_{t} = (Z_{t} + Z_{t + 1}) / (2 S_{t} \sum_{j} I_{t}^{j})

where we use the notation

X_{t} \sim X(t) : X = S, E^{i}, I^{j}, R, Z, B, \mu, \beta \\ \begin{aligned} Z(t) &= \int_{t - 1}^{t} \lambda(s) S(s) \, \text{d} s \\ B(t) &= \int_{t - 1}^{t} \nu(s) \, \text{d} s \end{aligned}

and it is understood that the independent variable t is a unitless measure of time relative to the spacing of the substituted time series of incidence and births.

Value

A “multiple time series” object, inheriting from class mts, with 1+m+n+1+1 columns (named S, E, I, R, and beta) storing the result of the iteration described in ‘Details’. It is completely parallel to argument series, having the same tsp attribute.

References

Jagan, M., deJonge, M. S., Krylova, O., & Earn, D. J. D. (2020). Fast estimation of time-varying infectious disease transmission rates. PLOS Computational Biology, 16(9), Article e1008124, 1-39. doi:10.1371/journal.pcbi.1008124

Examples

if (requireNamespace("adaptivetau")) withAutoprint({

data(seir.ts02, package = "fastbeta")
a <- attributes(seir.ts02)
str(seir.ts02)
plot(seir.ts02)

## We suppose that we have perfect knowledge of incidence,
## births, and the data-generating parameters
series <- cbind(seir.ts02[, c("Z", "B")], mu = a[["mu"]](0))
colnames(series) <- c("Z", "B", "mu") # FIXME: stats:::cbind.ts mangles dimnames

args <- c(list(series = series),
          a[c("sigma", "gamma", "delta", "m", "n", "init")])
str(args)

X <- do.call(fastbeta, args)
str(X)
plot(X)

plot(X[, "beta"], ylab = "Transmission rate")
lines(a[["beta"]](time(X)), col = "red") # the "truth"

})

Defunct Functions in Package fastbeta

Description

The functions and other objects listed here are no longer part of fastbeta as they are no longer needed.

Usage

## Nothing yet!

Details

These either are stubs reporting that they are defunct or have been removed completely (apart from being documented here).

See Also

Deprecated, base-deprecated, fastbeta-deprecated, fastbeta-notyet.

Deprecated Functions in Package fastbeta

Description

The functions and other objects listed here are provided only for compatibility with older versions of fastbeta and may become defunct as soon as the next release.

Usage

## Nothing yet!

See Also

Defunct, base-defunct, fastbeta-defunct, fastbeta-notyet.

Not Yet Implemented Functions in Package fastbeta

Description

The functions listed here are defined but not yet implemented. Use bug.report(package = "fastbeta") to request an implementation.

Usage

## Nothing yet!

See Also

NotYetImplemented, fastbeta-deprecated, fastbeta-defunct.

Parametric Bootstrapping

Description

A simple wrapper around fastbeta using it to generate a “primary” estimate of a time-varying transmission rate and r bootstrap estimates. Bootstrap estimates are computed for incidence time series simulated using seir, with transmission rate defined as the linear interpolant of the primary estimate.

Usage

fastbeta.bootstrap(r,
                   series, sigma = 1, gamma = 1, delta = 0,
                   m = 1L, n = 1L, init, ...)

Arguments

Value

A “multiple time series” object, inheriting from class mts, with 1+r columns storing the one primary and r bootstrap estimates. It is completely parallel to argument series, having the same tsp attribute.

Examples

if (requireNamespace("adaptivetau")) withAutoprint({

data(seir.ts02, package = "fastbeta")
a <- attributes(seir.ts02)
str(seir.ts02)
plot(seir.ts02)

## We suppose that we have perfect knowledge of incidence,
## births, and the data-generating parameters
series <- cbind(seir.ts02[, c("Z", "B")], mu = a[["mu"]](0))
colnames(series) <- c("Z", "B", "mu") # FIXME: stats:::cbind.ts mangles dimnames

args <- c(list(r = 100L, series = series),
          a[c("sigma", "gamma", "delta", "m", "n", "init")])
str(args)

R <- do.call(fastbeta.bootstrap, args)
str(R)
plot(R)
plot(R, level = 0.95)

})

Calculate Coefficient Matrix for Iteration Step

Description

Calculates the coefficient matrix corresponding to one step of the iteration carried out by fastbeta:

y <- c(1, E, I, R, S)
for (pos in seq_len(nrow(series) - 1L)) {
    L <- fastbeta.matrix(pos, series, ...)
    y <- L %*% y
}

Usage

fastbeta.matrix(pos,
                series, sigma = 1, gamma = 1, delta = 0,
                m = 1L, n = 1L)

Arguments

Value

A lower triangular matrix of size 1+m+n+1+1.

Examples

if (requireNamespace("adaptivetau")) withAutoprint({

data(seir.ts02, package = "fastbeta")
a <- attributes(seir.ts02); p <- length(a[["init"]])
str(seir.ts02)
plot(seir.ts02)

## We suppose that we have perfect knowledge of incidence,
## births, and the data-generating parameters
series <- cbind(seir.ts02[, c("Z", "B")], mu = a[["mu"]](0))
colnames(series) <- c("Z", "B", "mu") # FIXME: stats:::cbind.ts mangles dimnames

args <- c(list(series = series),
          a[c("sigma", "gamma", "delta", "init", "m", "n")])
str(args)

X <- unclass(do.call(fastbeta, args))[, seq_len(p)]
colnames(X)
Y <- Y. <- cbind(1, X[, c(2L:p, 1L)], deparse.level = 2L)
colnames(Y)

args <- c(list(pos = 1L, series = series),
          a[c("sigma", "gamma", "delta", "m", "n")])
str(args)

L <- do.call(fastbeta.matrix, args)
str(L)
symnum(L != 0)

for (pos in seq_len(nrow(series) - 1L)) {
    args[["pos"]] <- pos
    L. <- do.call(fastbeta.matrix, args)
    Y.[pos + 1L, ] <- L. %*% Y.[pos, ]
}
stopifnot(all.equal(Y, Y.))

})

Peak to Peak Iteration

Description

Approximates the state of an SEIR model at a reference time from an equally spaced, T-periodic incidence time series and other data. The algorithm relies on a strong assumption: that the incidence time series was generated by the asymptotic dynamics of an SEIR model admitting a locally stable, T-periodic attractor. Hence do interpret with care.

Usage

ptpi(series, sigma = 1, gamma = 1, delta = 0,
     m = 1L, n = 1L, init,
     start = tsp(series)[1L], end = tsp(series)[2L],
     tol = 1e-03, iter.max = 32L,
     backcalc = FALSE, complete = FALSE, ...)

Arguments

Details

ptpi can be understood as an iterative application of fastbeta to a subset of series. The basic algorithm can be expressed in R code as:

w <- window(series, start, end); i <- nrow(s); j <- seq_along(init)
diff <- Inf; iter <- 0L
while (diff > tol && iter < iter.max) {
    init. <- init
    init <- fastbeta(w, sigma, gamma, delta, m, n, init)[i, j]
    diff <- sqrt(sum((init - init.)^2) / sum(init.^2))
    iter <- iter + 1L
}
value <- init

Back-calculation involves solving a linear system of equations; the back-calculated result can mislead if the system is ill-conditioned.

Value

A list with elements:

References

Jagan, M., deJonge, M. S., Krylova, O., & Earn, D. J. D. (2020). Fast estimation of time-varying infectious disease transmission rates. PLOS Computational Biology, 16(9), Article e1008124, 1-39. doi:10.1371/journal.pcbi.1008124

Examples

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

data(seir.ts01, package = "fastbeta")
a <- attributes(seir.ts01); p <- length(a[["init"]])
str(seir.ts01)
plot(seir.ts01)

## We suppose that we have perfect knowledge of incidence,
## births, and the data-generating parameters, except for
## the initial state, which we "guess"
series <- cbind(seir.ts01[, c("Z", "B")], mu = a[["mu"]](0))
colnames(series) <- c("Z", "B", "mu") # FIXME: stats:::cbind.ts mangles dimnames

plot(series[, "Z"])
start <- 23; end <- 231
abline(v = c(start, end), lty = 2)

set.seed(0L)
args <- c(list(series = series),
          a[c("sigma", "gamma", "delta", "m", "n", "init")],
          list(start = start, end = end, complete = TRUE))
init <- seir.ts01[which.min(abs(time(seir.ts01) - start)), seq_len(p)]
args[["init"]] <- init * rlnorm(p, 0, 0.1)
str(args)

L <- do.call(ptpi, args)
str(L)

S <- L[["x"]][, "S", ]
plot(S, plot.type = "single")
lines(seir.ts01[, "S"], col = "red", lwd = 4) # the "truth"
abline(h = L[["value"]]["S"], v = start, col = "blue", lwd = 4, lty = 2)

## Relative error
L[["value"]] / init - 1

})

Simulate Infectious Disease Time Series

Description

Simulates incidence time series based on an SEIR model with user-defined forcing and a simple model for observation error.

Note that simulation code depends on availability of suggested packages adaptivetau and deSolve. If the dependency cannot be loaded then an error is signaled.

Usage

seir(length.out = 1L,
     beta, nu = function (t) 0, mu = function (t) 0,
     sigma = 1, gamma = 1, delta = 0,
     m = 1L, n = 1L, init,
     stochastic = TRUE, prob = 1, delay = 1,
     aggregate = FALSE, useCompiled = TRUE, ...)

## A basic wrapper for the m=0L case:

 sir(length.out = 1L,
     beta, nu = function (t) 0, mu = function (t) 0,
     gamma = 1, delta = 0,
     n = 1L, init,
     stochastic = TRUE, prob = 1, delay = 1,
     aggregate = FALSE, useCompiled = TRUE, ...)

Arguments

Details

Simulations are based on an SEIR model with

• m latent stages (E^{i}, i = 1,\ldots,m);

• n infectious stages (I^{j}, j = 1,\ldots,n);

• time-varying rates \beta, \nu, and \mu of transmission, birth, and natural death; and

• constant rates m \sigma, n \gamma, and \delta of removal from each latent, infectious, and recovered compartment, where removal from the recovered compartment implies return to the susceptible compartment (loss of immunity).

seir(stochastic = FALSE) works by numerically integrating the system of ordinary differential equations

\begin{alignedat}{10} \text{d} & S &{} / \text{d} t &{} = {}& \delta &R &{} - ( && \lambda(t) &{} + \mu(t)) S &{} + \nu(t) \\ \text{d} & E^{ 1} &{} / \text{d} t &{} = {}& \lambda(t) &S &{} - ( && m \sigma &{} + \mu(t)) E^{ 1} &{} \\ \text{d} & E^{i + 1} &{} / \text{d} t &{} = {}& m \sigma &E^{i} &{} - ( && m \sigma &{} + \mu(t)) E^{i + 1} &{} \\ \text{d} & I^{ 1} &{} / \text{d} t &{} = {}& m \sigma &E^{m} &{} - ( && n \gamma &{} + \mu(t)) I^{ 1} &{} \\ \text{d} & I^{j + 1} &{} / \text{d} t &{} = {}& n \gamma &I^{j} &{} - ( && n \gamma &{} + \mu(t)) I^{j + 1} &{} \\ \text{d} & R &{} / \text{d} t &{} = {}& n \gamma &I^{n} &{} - ( && \delta &{} + \mu(t)) R &{} \end{alignedat} \\ \lambda(t) = \beta(t) \sum_{j} I^{j}

where it is understood that the independent variable t is a unitless measure of time relative to an observation interval. To get time series of incidence and births, the system is augmented with two equations describing cumulative incidence and births

\begin{aligned} \text{d} Z / \text{dt} &{} = \lambda(t) S \\ \text{d} B / \text{dt} &{} = \nu(t) \end{aligned}

and the augmented system is numerically integrated. Observed incidence is simulated from incidence by scaling the latter by prob and convolving the result with delay.

seir(stochastic = TRUE) works by simulating a Markov process corresponding to the augmented system, as described in the reference. Observed incidence is simulated from incidence by binning binomial samples taken with probabilities prob over future observation intervals according to multinomial samples taken with probabilities delay.

Value

A “multiple time series” object, inheriting from class mts. Beneath the class, it is a length.out-by-(1+m+n+1+2) numeric matrix with columns S, E, I, R, Z, and B, where Z and B specify incidence and births as the number of infections and births since the previous time point.

If prob or delay is not missing, then there is an additional column Z.obs specifying observed incidence as the number of infections observed since the previous time point. The first length(delay) elements of this column contain partial counts.

References

Cao, Y., Gillespie, D. T., & Petzold, L. R. (2007). Adaptive explicit-implicit tau-leaping method with automatic tau selection. Journal of Chemical Physics, 126(22), Article 224101, 1-9. doi:10.1063/1.2745299

See Also

seir.auxiliary, seir.library.

Examples

if (requireNamespace("adaptivetau")) withAutoprint({

beta <- function (t, a = 1e-01, b = 1e-05) b * (1 + a * sinpi(t / 26))
nu   <- function (t) 1e+03
mu   <- function (t) 1e-03

sigma <- 0.5
gamma <- 0.5
delta <- 0

init <- c(S = 50200, E = 1895, I = 1892, R = 946011)

length.out <- 250L
prob <- 0.1
delay <- diff(pgamma(0:8, 2.5))

set.seed(0L)
X <- seir(length.out, beta, nu, mu, sigma, gamma, delta, init = init,
          prob = prob, delay = delay, epsilon = 0.002)
##                                              ^^^^^
## default epsilon = 0.05 allows too big leaps => spurious noise
##
str(X)
plot(X)

r <- 10L
Y <- do.call(cbind, replicate(r, simplify = FALSE,
	seir(length.out, beta, nu, mu, sigma, gamma, delta, init = init,
	     prob = prob, delay = delay, epsilon = 0.002)[, "Z.obs"]))
str(Y) # FIXME: stats:::cbind.ts mangles dimnames
plot(window(Y, start = tsp(Y)[1L] + length(delay) / tsp(Y)[3L]),
     ##        ^^^^^
     ## discards points showing edge effects due to 'delay'
     ##
     plot.type = "single", col = seq_len(r), ylab = "Case reports")

})

Auxiliary Functions for the SEIR Model without Forcing

Description

Calculate the basic reproduction number, endemic equilibrium, and Jacobian matrix of the SEIR model without forcing.

Usage

seir.R0      (beta, nu = 0, mu = 0, sigma = 1, gamma = 1, delta = 0,
              m = 1L, n = 1L, N = 1)
seir.ee      (beta, nu = 0, mu = 0, sigma = 1, gamma = 1, delta = 0,
              m = 1L, n = 1L, N = 1)
seir.jacobian(beta, nu = 0, mu = 0, sigma = 1, gamma = 1, delta = 0,
              m = 1L, n = 1L)

Arguments

Details

If \mu, \nu = 0, then the basic reproduction number is computed as

\mathcal{R}_{0} = N \beta / \gamma

and the endemic equilibrium is computed as

\begin{bmatrix} S^{\hphantom{1}} \\ E^{i} \\ I^{j} \\ R^{\hphantom{1}} \end{bmatrix} = \begin{bmatrix} \gamma / \beta \\ w \delta / (m \sigma) \\ w \delta / (n \gamma) \\ w \end{bmatrix}

where w is chosen so that the sum is N.

If \mu, \nu > 0, then the basic reproduction number is computed as

\mathcal{R}_{0} = \nu \beta a^{-m} (1 - b^{-n}) / \mu^{2}

and the endemic equilibrium is computed as

\begin{bmatrix} S^{\hphantom{1}} \\ E^{i} \\ I^{j} \\ R^{\hphantom{1}} \end{bmatrix} = \begin{bmatrix} \mu a^{m} / (\beta (1 - b^{-n})) \\ w a^{m - i} b^{n} (\delta + \mu) / (m \sigma) \\ w b^{n - j} (\delta + \mu) / (n \gamma) \\ w \end{bmatrix}

where w is chosen so that the sum is \nu / \mu, the population size at equilibrium, and a = 1 + \mu / (m \sigma) and b = 1 + \mu / (n \gamma).

Currently, none of the functions documented here are vectorized. Arguments must have length 1.

Value

seir.R0 returns a numeric vector of length 1. seir.ee returns a numeric vector of length 1+m+n+1. seir.jacobian returns a function of one argument x (which must be a numeric vector of length 1+m+n+1) whose return value is a square numeric matrix of size length(x).

See Also

seir, for the system of ordinary differential equations on which these computations are predicated.

Often Used Simulations

Description

Infectious disease time series simulated using seir, for use primarily in examples, tests, and vignettes. Users should not rely on simulation details, which may change between package versions.

Note that simulation code depends on availability of suggested packages adaptivetau and deSolve. If the dependency cannot be loaded then the value of the data set is NULL.

Usage

## if (requireNamespace("deSolve"))
data(seir.ts01, package = "fastbeta")
## else ...

## if (requireNamespace("adaptivetau"))
data(seir.ts02, package = "fastbeta")
## else ...

Format

A “multiple time series” object, inheriting from class mts, always a subset of the result of a call to seir, discarding transient behaviour. Simulation parameters may be preserved as attributes.

Source

Scripts sourced by data to reproduce the simulations are located in subdirectory ‘data’ of the fastbeta installation; see, e.g. system.file("data", "seir.ts01.R", package = "fastbeta").

See Also

seir.

Examples

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

data(seir.ts01, package = "fastbeta")
str(seir.ts01)
plot(seir.ts01)

})

if (requireNamespace("adaptivetau")) withAutoprint({

data(seir.ts02, package = "fastbeta")
str(seir.ts02)
plot(seir.ts02)

})

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

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

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)

})

Smallpox Mortality in London, England, 1661-1930

Description

Time series of deaths due to smallpox, deaths due to all causes, and births in London, England, from 1661 to 1930, as recorded in the London Bills of Mortality and the Registrar General's Weekly Returns.

Usage

data(smallpox, package = "fastbeta")

Format

A data frame with 13923 observations of 5 variables:

from

start date of the record.

nday

length of the record, which is the number of days (typically 7) over which deaths and births were counted.

smallpox

count of deaths due to smallpox.

allcauses

count of deaths due to all causes.

births

count of births.

Source

A precise description of the data set and its correspondence to the original source documents is provided in the reference.

A script generating the smallpox data frame from a CSV file accompanying the reference is available as system.file("scripts", "smallpox.R", package = "fastbeta").

References

Krylova, O. & Earn, D. J. D. (2020). Patterns of smallpox mortality in London, England, over three centuries. PLOS Biology, 18(12), Article e3000506, 1-27. doi:10.1371/journal.pbio.3000506

Examples

data(smallpox, package = "fastbeta")
str(smallpox)
table(smallpox[["nday"]]) # not all 7 days, hence:
plot(7 * smallpox / as.double(nday) ~ from, smallpox, type = "l")