Constraining \(R(t)\) over time

Overview

Using the renewal equation and given that \(I(t)\) and \(ω\) are known, \(R(t)\) can be solved for algebraically starting with \(R(t=1)\) and iterating forwards in time. However, this will produce highly volatile estimates of \(R(t)\) that recover the incidence curve directly.

Solving directly for \(R(t)\) at every timestep is undesirable for several reasons:

• observed infections are the result of a noisy process with random day-to-day variations, which do not reflect an underlying change in infectivity;
• real-world infection data are rarely complete, especially in an emerging epidemic, meaning that a certain amount of uncertainty must be incorporated into any estimation framework;
• infection incidences, \(I(t)\), are the data of interest but are cannot be observed directly, so many calculations instead use the observed reported cases, \(C(t)\), which requires some additional processing to incorporate into calculations of \(R(t)\).

Therefore a variety of methods for estimating \(R(t)\) exist that impose structure on how \(R(t)\) varies with time.

Each method confers various assumptions and implications for resulting estimates of \(R(t)\), and new methods represent a large area of innovation with regards to real-time infectious disease modeling. With these constraints, we can make inference from sampled case-count data as a signal of unobserved infections in the larger unobserved population.

Fixed sliding windows

A straightforward method of imposing structure on \(R(t)\) over time involves constraining \(R(t)\) to be drawn from the same distribution within moving time subsets, called sliding windows.

We add the prefix of “fixed-size” to distinguish from methods that may adapt the size of the sliding window over time.

Consider the scenario where \(I(t)\) are drawn from a series of Poisson distributions and where \(R(t)\) are drawn from a series of Gamma distributions. Using a sliding window size, \(τ\), of 5 days:

• \(R(t)\) on days 2 to 6 are assumed to be drawn from the Gamma distribution with parameters \(a_1\) and \(b_1\)

• \(R(t)\) on days 3 to 7 are drawn from a Gamma distribution with parameters \(a_2\) and \(b_2\)

In the above scenario, days 3 through 6 are in both windows and thus will be values that could be reasonably drawn from Gamma distributions with either \(a_1\) and \(b_1\) or \(a_2\) and \(b_2\).

Using an assumption of Gamma distributions for the prior distribution of \(ω\) and \(R(t)\), Cori et al. (2013) analytically derived a posterior distribution \(R(t)\) using fixed-size sliding windows. This distribution has the following directly calculated (rather than inferred) mean and coefficient of variation of \(R(t)\):

\[ E[R(t)]= \frac{a + \sum_{i = \max(1, t - τ)}^{t} I(t)} { 1/b + \sum_{i = \max(1, t - τ)}^{t} Λ(t)} \]

\[ C.V. [R(t)]=\frac{1}{a + \sum_{i = \max(1, t - τ)}^{t} I(t)} \]

Thus, sliding windows with larger \(τ\) improve the stability of the estimate of \(R(t)\) (as compared to smaller \(τ\)) because the coefficient of variation of \(R(t)\) decreases as number of infections increases (see Web Appendix 1 of Cori et al. (2013)). Fixed sliding windows are a key feature of the EpiEstim package.

Drawbacks of fixed sliding windows

There are limitations of the fixed sliding window approach, articulated well in Gostic et al. (2020) and summarized here:

• There is no posterior distribution for the expected value of incidence.

• In the fixed size sliding window approach, \(τ\) must be explicitly defined prior to inference. Shorter \(τ\) will lead to quicker response but more variable estimates of \(R(t)\), which increases the risk of over-fitting. At the extreme, if the \(τ\) is set to 1 day, the resulting \(R(t)\) will recover exactly the input case data. The default recommendation for \(τ\) is one week (7 days)

• In addition, there is debate in the literature about where in time the estimate of \(R(t)\) for each window should go: Gostic et al. (2020) recommends using the midpoint of each sliding window rather than time \(t\).

• The choice of both \(τ\) and the location of the estimate of \(R(t)\) within each window results in gaps in predictions for \(R(t)\):

  • at the beginning of the modeling period, if not enough cases have occurred, Web Index 4 of Cori et al. (2013) suggests that calculation of R(t) should start at least one serial interval’s worth of time after the 12th case is recorded

  • at the end of the modeling period, to account for reporting delays or time between the midpoint of τ and the end of τ.

Modifications of fixed sliding windows in EpiEstim

Several packages modify the functionality of EpiEstim. These modifications include:

• fitting \(τ\) to minimize the Accumulated Prediction Error, as in the package APEestim

• fitting \(R(t)\) via harmonizing the EpiEstim method and methods from the case reproduction number in a variational method, as in the EpiInvert package.

• and dissaggregation of data into smaller time units, as in the ern package.

Random walk

Another method of constraining how \(R(t)\) evolves in time is to define the relationship between \(R(t)\), infections, and time in a random walk or auto-regressive framework.

In this framework, there are latent or unobserved variables, e.g., \(R(t)\), that depend on observed variables, e.g., \(I(t)\) via the renewal equation, and the evolution of the unobserved variables through time can be parameterized.

The auto-regressive component means that the current value of \(R(t)\) is correlated via some mechanism with \(R(t-1)\) (and potentially other past values).

The packages EpiNow2, epinowcast, and WhiteLabRt contain an implementations of a random walk procedure that look generally as follows:

\[ f(R(t))=f(R(t-1))+N(0,σ_R) \]

where \(σ_R\) has a user-defined prior, e.g., with hyperparameters \(ρ\) and \(φ\): \[ σ_R \sim HalfNormal(ρ,φ) \]

The random walk implies that adjacent \(R(t)\) values may be drawn from similar or even the same distribution, and would be correlated in time based on previous values.

The function \(f\) can be a transformation of \(R(t)\), e.g. in log space to correct for the skewness of \(R(t)\). This can help provide a variable that is more Gaussian, provide a variable that obeys the properties that we expect from \(R(t)\) (i.e., is non-negative), and aid in interpretability.

Filtering

Filtering is another way that \(R(t)\) is constrained in time. Filtering is similar to a random walk constraint, in that a hyperparameter controls the amount of difference between \(R(t)\) in adjacent time-steps, but with a different functional form to a random walk.

One way that a filter could be implemented is in a Hidden Markov Model (Rabiner and Juang (1986)). A simple forward-looking linear filter for \(R(t)\) in an Hidden Markov Model might look as follows, with a tuning parameter (\(η\)) to influence the amount that \(R(t)\) can vary between time-steps and a standard white noise component (\(ϵ\)):

\[ R(t)=R(t-1)+(η√(R(t-1))) * ϵ(t-1) \]

The package EpiFilter implements a two-stage filtering and smoothing method for estimating \(R(t)\). A key innovation of EpiFilter is that the states of historical \(R(t)\) are constrained to a predefined set of values; this dramatically reduces calculation time. The smoothing stage refines estimates of \(R(t)\) by incorporating future incidence, in this way using all available data in estimates of historical \(R(t)\). These modeling steps help avoid \(R(t)\) instability when infections are low and instability at the beginning and (more importantly) the end of the modeling period.

Another way that filtering can be implemented is across the entire R(t) time-series, as in the packages RtEstim and EpiLPS.

Notably, \(R(t)\) estimated in this way thus contains information about past and pending infections, e.g., for \(R(t=i)\), the smoothing step will affect \(R(t=i)\) using information from \(0<i≤t_{max}\). This complicates comparisons to outputs from other methods that only use historical information to estimate \(R(t)\), e.g., estimates for \(R(t=i)\) containing only information from \(t<i\).

Gausian Process models

Gaussian Process models (Schulz et al. (2018)) are yet another flexible method of constraining the evolution of \(R(t)\) in time than the methods discussed thus far.

In Gaussian Process modeling, a family of basis functions are fit to available data, permitting inference about continuous processes without needing to a priori define where inflection points occur.

The core of Gaussian Process operations is a kernel, which is used to assess the similarity between input vectors, say \(x\) and \(x'\).

There are many options for potential kernels, and each contains different hyperparameters that are used to control the amount of smoothing that is enforced, as well as other factors. One such choice is the squared exponential kernel:

\[ k(x,x') = α^2 \exp ( \frac {-(x -x')^2} {2l^2} ) \]

In this kernel, the hyperparameters are the length scale, \(l\), which controls the smoothness of the model, and the magnitude, α, which controls the range of values used in the fitting process. These parameters can be given prior distributions and fit using optimization.

EpiNow2 contains options to use Gaussian Process models to control how R(t) evolves in time. As one example, the relationship between first difference values of R(t) can be constrained using a zero-mean Gaussian Process model with the above kernel as the covariance function:

\[ \log R(t) = \log R(t-1) + GP(0, k(R(t), R(t')) ) \]

The advantage of Gaussian Process models is that R(t) is enforced to change smoothly in time . Limitations include complexity and computational time (Riutort-Mayol et al. (2020)), which in general means Gaussian Process runtimes and required computational resources are considerable as compared to other methods.

Cori A, Ferguson NM, Fraser C, Cauchemez S. 2013. A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics. American Journal of Epidemiology 178:1505–1512; doi:10.1093/aje/kwt133.

Gostic KM, McGough L, Baskerville EB, Abbott S, Joshi K, Tedijanto C, et al. 2020. Practical considerations for measuring the effective reproductive number, Rt. PLOS Computational Biology 16:e1008409; doi:10.1371/journal.pcbi.1008409.

Rabiner L, Juang B. 1986. An introduction to hidden markov models. ieee assp magazine 3: 4–16.

Riutort-Mayol G, Bürkner P, Andersen M, Solin A, Vehtari A. 2020. Practical hilbert space approximate bayesian gaussian processes for probabilistic programming. ArXiv prepr. arXiv preprint arXiv:200411408.

Schulz E, Speekenbrink M, Krause A. 2018. A tutorial on gaussian process regression: Modelling, exploring, and exploiting functions. Journal of mathematical psychology 85: 1–16.