ptpi function

Peak to Peak Iteration

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.

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

• series: a multiple time series object, inheriting from class mts, with three columns storing (parallel , equally spaced) time series of incidence, births, and the per capita natural mortality rate, in that order.
• sigma, gamma, delta: non-negative numbers. m*sigma, n*gamma, and delta are the rates of removal from each latent, infectious, and recovered compartment.
• m: a non-negative integer indicating a number of latent stages.
• n: a positive integer indicating a number of infectious stages.
• init: a numeric vector of length 1+m+n+1 giving an initial guess for the state at time start.
• start, end: start and end times for the iteration, whose difference should be approximately equal to an integer number of periods. One often chooses the time of the first peak in the incidence time series and the time of the last peak in phase with the first.
• tol: a tolerance indicating a stopping condition; see Details .
• iter.max: the maximum number of iterations.
• backcalc: a logical indicating if the state at time tsp(series)[1] should be back-calculated from the state at time start if that is later.
• complete: a logical indicating if intermediate states should be recorded in an array. Useful mainly for didactic or diagnostic purposes.
• ...: optional arguments passed to deconvolve, if the first column of series represents observed incidence rather than actual or estimated incidence.

Returns

A list with elements: - value: an approximation of the state at time start or at time tsp(series)[1L], depending on backcalc.

• diff: the relative difference between the last two approximations.

• iter: the number of iterations performed.

• x: if complete = TRUE, then a multiple time series

  object, inheriting from class mts, with dimensions c(nrow(w), length(value), iter), where w = window(series, start, end). x[, , k] contains the state at each time(w) in iteration k.

Details

ptpi can be understood as an iterative application of fastbeta to a subset of series. The basic algorithm can be expressed in 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.

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. tools:::Rd_expr_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

})