Version:	0.5.0
Date:	2025-10-15
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-10-15 14:09:35 UTC; mikael
Author:	Mikael Jagan [aut, cre]
Maintainer:	Mikael Jagan <jaganmn@mcmaster.ca>
Repository:	CRAN
Date/Publication:	2025-10-15 14:30:02 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

Combine Time Series Objects

Description

A replacement for the S3 method registered by package stats for generic function cbind and class ts. It sets column names following the rules employed by the internal default method for cbind. It exists to allow users to work around PR#18583, which shows that the method in package stats employs different and rather inconsistent rules. This function must be called directly, as it is not registered as a method for cbind.

Usage

## S3 method for class 'ts'
cbind(..., deparse.level = 1)

Arguments

Value

A “multiple time series” object, inheriting from class mts.

Examples

n <- 3L
x <- matrix(0, n, n, dimnames = list(NULL, LETTERS[seq_len(n)]))
y <- seq_len(n)
tsx <- ts(x)
tsy <- ts(y)
`~` <- identity
for (k in 0L:2L) {
	cat(sprintf("k = %d:\n\n\n", k))
	withAutoprint({
	try(colnames(cbind   (   x,    y, deparse.level = k)))
	try(colnames(cbind   ( tsx,  tsy, deparse.level = k)))
	try(colnames(cbind.ts( tsx,  tsy, deparse.level = k)))
	try(colnames(cbind   (~  x, ~  y, deparse.level = k)))
	try(colnames(cbind   (~tsx, ~tsy, deparse.level = k)))
	try(colnames(cbind.ts(~tsx, ~tsy, deparse.level = k)))
	})
	cat("\n\n")
}
rm(`~`)

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, x.pad = 0, tol = 1, iter.max = 32L, complete = FALSE)

Arguments

Details

Do note that temporal alignment of x (length n) and y=deconvolve(x, ...)$value (length or number of rows d+n) requires padding x on the left, as in cbind(x=c(rep(NA, d), x), y); see the examples.

Value

A list with elements:

subset(value, start == 0) is zero because zero is invariant under the iteration. If delay has l leading zeros and t trailing zeros, then head(value, t) and tail(value, l) are NaN due to divide-by-zero in the iteration. (Conceptually, x and delay provide no information about incidence during those intervals.)

References

Goldstein, E., Dushoff, J., Ma, J., Plotkin, J. B., Earn, D. J. D., & Lipsitch, M. (2009). 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

##
## Example 1: simulation
##

set.seed(2L)
n <- 200L
d <- 50L
p <- 0.1
prob <- plogis(rlogis(d + n, location = qlogis(p), scale = 0.1))
delay <- c(0, diff(pgamma(0L:d, 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, 3),
        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, 3),
       bty = "n")

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


##
## Example 2: application to pneumonia and influenza
##

data(pneumonia, package = "fastbeta")
x <- pneumonia[["series"]][["deaths"]]
delay <- pneumonia[["delay"]][["gpg"]]

n <- length(x)
d <- length(delay) - 1L
r <- 30L

D2 <- deconvolve(x = x, delay = delay, tol = 0, iter.max = r,
                 complete = TRUE)
stopifnot(D2[["iter"]] == r,
          identical(dim(D2[["value"]]), c(d + n, 1L + r)),
          length(D2[["chisq"]]) == 1L + r,
          min(D2[["chisq"]]) < 1)

## Subscript for the first, critical, and last values:
j2 <- c(1L, which.max(D2[["chisq"]] < 1), 1L + r)

matplot(x = seq(from = pneumonia[["series"]][1L, "date"] - d,
                by = 1, length.out = d + n),
        y = cbind(c(rep(NA, d), x), D2[["value"]][, j2]),
        type = "o",
        col = c(1L, 4L, 2L, 3L), pch = 1L, lty = 1L, lwd = 1,
        xlab = "1918", ylab = "deaths")
legend("topleft", NULL,
       c("observed", sprintf("after %d iterations", j2 - 1L)),
       col = c(1L, 4L, 2L, 3L), pch = 1L, lty = 1L, lwd = 1,
       bty = "n")

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.ts(seir.ts02[, c("Z", "B")], mu = a[["mu"]](0))

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.ts(seir.ts02[, c("Z", "B")], mu = a[["mu"]](0))

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.ts(seir.ts02[, c("Z", "B")], mu = a[["mu"]](0))

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

})

Pneumonia and Influenza Mortality in Philadelphia, PA, 1918

Description

Time series of deaths due to pneumonia and influenza in Philadelphia, PA from September 1, 1918 to December 31, 1918, as recorded in the “Special Tables of Mortality” of the U.S. Census Bureau.

Usage

data(pneumonia, package = "fastbeta")

Format

A named list with 2 components, series and delay. series is a data frame with 122 rows and 2 variables:

date

date of the record.

deaths

count of deaths due to influenza and pneumonia.

delay is a data frame with 64 rows and 3 variables:

nday

number of days from infection to death.

goldstein, gpg

probabilities, not summing to 1 due to rounding and truncation; see ‘Source’.

Source

A script generating the pneumonia object is available as system.file("scripts", "pneumonia.R", package = "fastbeta").

series is obtained from Table 2 in the first reference.

delay is obtained from the remaining references. Component goldstein is obtained from Figure 1 in the Supporting Information of Goldstein et al. (2009). Component gpg is obtained from the convolution of two gamma distributions, one for the time from infection to symptom onset fitted to Figure 1 in Moser et al. (1979) and another for the time from symptom onset to death fitted to Chart 2 in Keeton & Cushman (1918).

References

U.S. Census Bureau (1920). Special Tables of Mortality from Influenza and Pneumonia: Indiana, Kansas, and Philadelphia, PA. U.S. Department of Commerce. https://www.census.gov/library/publications/1920/demo/1918-mortality-special-tables.html

Goldstein, E., Dushoff, J., Ma, J., Plotkin, J. B., Earn, D. J. D., & Lipsitch, M. (2009). 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

Moser, M. R., Bender, T. R., Margolis, H. S., Noble, G. R., Kendal, A. P., & Ritter, D. G. (1979). An outbreak of influenza aboard a commercial airliner. Americal Journal of Epidemiology, 110(1), 1-6. doi:10.1093/oxfordjournals.aje.a112781

Keeton, R. W. & Cushman, A. B. (1918). The influenza epidemic in Chicago: the disease as a type of toxemic shock. Journal of the Americal Medical Association. 71(24), 1962-1967.

Examples

data(pneumonia, package = "fastbeta")
str(pneumonia)

plot(deaths ~ date, pneumonia$series, xlab = "1918")

plot(goldstein/sum(goldstein) ~ nday, pneumonia$delay, type = "o",
     lty = 2, pch = 1, xlab = "days", ylab = "probability")
lines(gpg/sum(gpg) ~ nday, pneumonia$delay, type = "o",
      lty = 1, pch = 16)

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.ts(seir.ts01[, c("Z", "B")], mu = a[["mu"]](0))

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.ts, 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)
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, 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 class 'sir.aoi'
summary(object, name = "Y", tol = 1e-6, ...)

Arguments

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.

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.

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)

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)

})

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 script generating the smallpox object from a CSV file accompanying the reference is available as system.file("scripts", "smallpox.R", package = "fastbeta").

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

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