Renewal models and inference problems

We consider an infectious disease epidemic observed over some time period 1 ≤ s ≤ t in a homogeneous and well-mixed population. While epidemics actually spread on dynamic networks involving stratified contact structures, homogeneous models can provide useful real-time insight into key transmission patterns and are more easily fit and verified with routine surveillance data such as incidence curves [19]. If the incidence or number of newly infected cases at time s is I_s then the set is the incidence curve of the epidemic. We assume that incidence is available on a daily scale so that is a vector of t daily counts but weeks or months could be used instead. A common problem in infectious disease is the inference of the transmissibility of the epidemic given this curve. The renewal model [1, 20] presents a general and popular framework for investigating this problem and its estimates, in some instances, can even approximate those from detailed network models [21].

The renewal model posits that epidemic transmissibility, summarised by the effective or instantaneous reproduction number, R_s, generates the observed incidence as in Eq (1). We assume incidence counts are Poisson distributed with Pois indicating Poisson noise. Here signifies equality in distributions and | means ‘conditioned on’. (1)

While Eq (1) does not directly model how susceptible individuals become infected, these effects are encoded in the reproduction number R_s, which measures the secondary cases generated per effective primary case at s [3]. The quantity , known as the total infectiousness, counts how many effective past cases are still infectious at s i.e., it describes the number of circulating cases that can actively transmit.

The generation time distribution of the epidemic, , controls how past incidence influences Λ_s, with w_u as the probability that a primary case takes between u − 1 and u days to generate a secondary case [1]. We make the standard assumption that is well approximated by the serial interval distribution of the epidemic of interest, which is known [11]. The serial interval is the time between symptom onset of a primary and its secondary case. While onset times are more practically measurable than actual infection times they do not include asymptomatic or subclinical cases. This approximation can limit inference (for example, serial intervals often have larger variances than generation times), but methods are being developed to improve its quality [22].

We focus on inferring the complete set of R_s values, denoted , given the incidence curve with t as the last recorded time. Estimating is important because changes in the values of instantaneous reproduction numbers often signify key transitions in epidemic transmissibility, which might be due to the imposition or relaxation of interventions. Instantaneous reproduction numbers are also the basis of other transmissibility metrics, such as cohort reproduction numbers or growth rates [3]. We define three main inference problems, based on how information in is recruited to construct every R_s estimate. We represent these problems in terms of the posterior distribution their solutions induce over possible R_s values. Estimates are functions of these posteriors.

The first is called filtering, where we sequentially compute the filtering posterior for every s ≤ t [23]. Filtering only uses incidence data up to time s, , for inferring R_s. Solving this problem is fundamental to real-time inference [24]. Filtering solutions are commonly employed for inferring instantaneous reproduction numbers. The second problem, which we call reverse-filtering, is the complement of the first. The reverse-filtering posterior , is important for retrospective or backward-looking estimates and infers R_s from incidence beyond s i.e., [25]. Practical R_s calculations do not involve reverse-filtering. Instead, the information used by r_s is implicit to deriving cohort reproduction numbers, [12]. Later we show that estimates from EpiEstim and the WT method are related to the p_s and r_s distributions, respectively.

The last problem, which is termed smoothing, is our main interest. It asks the question: how can we construct an R_s, at every s ≤ t, that integrates all available incidence information from . To solve this problem we must formulate the smoothing posterior distribution [14]. Functions of this posterior would then yield maximally informed instantaneous reproduction number estimates (and cohort reproduction number estimates by extension). Note that p_s, r_s and q_s depend on the choices of state space model, which describes the dynamics of R_s across time, and observation model, which explains how changes in R_s lead to trends in observed incidence data [18]. These models encode our assumptions about the epidemic of interest and determine how estimates trade off assumptions against data. We next explore how posterior distribution selection determines performance, especially when data are scarce, and examine how EpiEstim and the WT methods fit within this framework.

Inference methods and low incidence

We mostly detail EpiEstim as instantaneous reproduction number, R_s, estimates are the main focus of this work. EpiEstim assumes that the estimate of R_s at time s depends on a rolling past window of data, τ(s) = {s, s − 1, …, s − k + 1}, of size k [11]. Consequently, at time s, only size k subsets of the incidence, , and total infectiousness, , are considered informative about R_s. Their sums over this window are i_τ(s) = ∑_u∈τ(s) I_u and λ_τ(s) = ∑_u∈τ(s) Λ_u. When k = 1, EpiEstim only uses the most recent case count, I_s, (and Λ_s) to estimate R_s. While this maximises flexibility, it usually results in over-fitting and so larger windows are often employed to trade off estimate variance with bias [11, 26].

EpiEstim therefore effectively solves a filtering problem, as discussed in the previous section. The filtering posterior distribution produced by EpiEstim is restricted to the informative window τ(s) and denoted . This is parametrised by the shape-scale gamma posterior distribution in Eq (2). (2)

The posterior p_τ(s) results from combining a gamma prior distribution on R_s, , with a Poisson observation likelihood for the incidence (see Eq (1)), as detailed in [11, 13]. We never explicitly include Λ_s terms in our notation (e.g., we could write p_τ(s) as ) since they appear naturally when using Eq (1) and the key difference among our inference problems relate to I_s terms.

The posterior mean estimate from Eq (2) is constructed as , . The variance around this estimate is (a + i_τ(s))(c + λ_τ(s))⁻². The observation model is given by Eq (1) but the state space model of EpiEstim is not explicit. However, if R_τ(s) is the assumed average reproduction number in τ(s) (which is used to estimate R_s) then λ_τ(s) R_τ(s) = ∑_u∈τ(s) Λ_u R_u [13]. Thus, EpiEstim somewhat incorporates a linear moving average state space model, and assumes that the filtering distribution p_s ≈ p_τ(s) by deeming data outside τ(s), , as effectively uninformative [12]. Since p_τ(s) can be computed sequentially across an ongoing epidemic, EpiEstim provides real-time inference.

The WT method takes a complementary approach to EpiEstim, computing transmissibility over a forward-looking window γ(s) = {s, s + 1, …, s + k − 1} [2]. Often k = t − s + 1 i.e., the window extends to the last observed incidence. The WT method uses the observation model of Eq (1) and has an implicit moving average state model R_γ(s) = ∑_u∈γ(s) R_u w_u−s+1, which leads to its cohort reproduction number estimates [3]. As this method effectively uses future information [12], it implicitly involves approximating the reverse-filtering distribution (see previous section) as . We illustrate the information windows employed by the WT method and EpiEstim, as well as the complete filtering and reverse-filtering windows in Fig 1. The goodness of the windowed distributions p_τ(s) and r_γ(s) as approximations to the general posterior distributions p_s and r_s will depend on k and the appropriateness of the state model underlying each method [13].

While EpiEstim and the WT methods are powerful tools for inferring transmissibility in real-time and in retrospect, they have two main and related limitations, which necessarily reduce the reliability of their outputs [12]. First, their performance degrades as s gets close to 1 for EpiEstim and t for the WT method [11]. These edge or censoring effects correspond, at the extreme, to and , which are weakly informed posterior distributions. As a result, at the beginning of the incidence curve EpiEstim can be unreliable (and even unidentifiable). The WT method suffers similarly at the end of that curve [1] (see Fig 1).

The second limitation occurs in phases of the epidemic where incidence is low for a prolonged period of time [17, 26]. In these periods data are sparse and the quality of estimates depend on how well the method of choice mediates between the little available information and its inherent assumptions [15]. We illustrate this with EpiEstim, using the mean and variance of the posterior estimate from Eq (2), , which are defined above. If incidence is small over the window τ(s) then the sum of incidence, i_τ(s), and the total infectiousness λ_τ(s) shrink, meaning that the prior hyperparameters, a and c, strongly influence the resulting estimate mean and variance. This contrasts the data-rich scenario when an epidemic is large, where i_τ(s) and λ_τ(s) overpower a and c.

Further, should a sequence of n ≥ k zero-incidence days occur, then i_τ(s) = 0 and λ_τ(s) → 0 as n increases, with a rate controlled by the serial interval of the epidemic. Here the shape parameter of p_τ(s) is exactly that of the prior distribution . As λ_τ(s) decays, (see Eq (2)), R_s becomes statistically unidentifiable from the window of data and inference is completely prior driven [26, 27]. While lack of data is a fundamental limitation, the point at which we lose inferential power is not fixed, and depends on the window size, k. Studies that formally optimised k for estimate reliability, found that small k is needed to infer sharp changes in transmissibility (e.g., due to lockdowns) [13], indicating that these issues can be acute. Analogous effects occur in the WT method if there are few incident cases across its forward-looking window γ(s) [12].

These prior-driven scenarios are realistic for epidemics in waning or tail phases, and can precede either elimination (i.e., epidemic extinction) or resurgence [28]. While some estimate degradation is guaranteed for any R_s inference method when faced with either edge-effects or low incidence, robustness can still be improved. Edge-effects can be largely overcome by constructing the smoothed posterior distribution for estimating the instantaneous reproduction number R_s, denoted q_s. Solving the smoothing problem melds the advantages of the opposite looking windows of EpiEstim and the WT method, removing the vulnerability near the ends of the incidence curve . This follows as q₁ = r₁ and q_t = p_t (see Fig 1). Further, by maximising the information used for inferring every R_s and by minimising our state model assumptions, we can ameliorate the impact of low incidence. We next develop a method, termed EpiFilter, to realise these improvements.

Bayesian (forward) recursive filtering

We reformulate the inference problem of estimating instantaneous reproduction numbers R_s from past incidence as an optimal Markov state filtering problem. Filtering describes a general class of engineering problems aimed at optimally, usually in a MSE sense, inferring some hidden state in real time from noisy observations [14, 18]. Given some functions f_s and g_s, which describe the state (R_s in our case) space dynamics and the process of generating noisy observations (the I_s here), the filtering problem tries to construct the posterior distribution p_s (see previous section) [23], which EpiEstim approximates. The conditional mean estimate leads to the minimum MSE of [23], which depends on all the past information.

The famed Kalman filter [24] was the genesis of these methods. Here we focus on Bayesian recursive filters for models with noisy count observations. These generalise the Kalman filter [23] and have been successfully applied to similar problems in phylodynamics and computational biology [29–31]. We reconsider our renewal model inference problem within this engineering state-observation framework, as described in Eq (3) [14]. (3)

Here R_s is the hidden Markov state that we wish to infer. It dynamically depends on the previous state R_s−1 and a noise term ϵ_s−1 via f_s. Observation I_s is then elicited due to R_s and a noise term ν_s, according to g_s [32].

We develop our filter under two very mild assumptions. First, we define some closed space, , over which R_s is valid. For a given resolution m, extrema R_min and R_max, and grid size δ = m⁻¹(R_max − R_min) then . This means the instantaneous reproduction number R_s must take a discrete value in , the i^th element of which is denoted . We formalise this notion in Eq (4a). (4a) (4b)

This is not restrictive since we can compute our filter for large m if needed and usually we are only interested in R_s on a coarse scale (e.g., policymakers may only want to know if R_s ≤ 1 or not). Other approaches, which depend on MCMC or related sampling methods (e.g., [9] and some implementations of EpiEstim), all implicitly assume some discretisation [29]. In Fig A of the S1 Appendix we show that convergence often occurs at small m.

Second, we propose a linear model for f_s, as defined in Eq (4b). There ϵ_s−1 is a standard white noise term i.e., with Norm signifying a normal distribution and η as a free parameter. We assume that a noisy linear projection of states over consecutive time-points provides a good approximation of the state trajectory. Not only is this assumption standard in engineering [23] and epidemiology [33] but it is also more flexible than the state model inherent to EpiEstim and the WT method. We scale the noise of this projection by a fraction, η < 1, of the magnitude of . This parameter controls the correlation among successive instantaneous reproduction numbers (and hence the state noise) but is set to ensure that R_s is a-priori non-negative.

Our observation model, g_s, is implicit and leads to the probability law in Eq (1). As a result, both our observations and state models are discrete (see Fig 1 for summaries). Because the state model governing R_s in Eq (4b) is stochastic, Eq (1) actually describes an over-dispersed (doubly stochastic) Poisson incidence curve. Consequently, η, allows us to partially model some heterogeneity in transmissibility, and may increase robustness to violations of the well-mixed assumption inherent to renewal models [21]. We can optimise our choice of η value by minimising the square error of the one-step-ahead predictions that result from our observation model [13].

We now define the Bayesian recursive filtering procedure, which is a main contribution of this work, and can be solved exactly, in real-time and with minimal computational effort. We adapt general recursive filtering equations from [14, 18, 25, 32], which are valid for various types of observation and state models, to our renewal model inference problem. The proof of the equations we employ can be found in these works. While we solve discrete, univariate problems (our state model is one dimensional), extensions to continuous-time, multivariate problems also exist [23, 30]. These recursive equations can also be approximately solved using particle filters [14, 32].

Recursive filtering involves two steps: prediction and correction. The first, given in Eq (5a), constructs a sequential prior predictive distribution, . This is informed by past incidence data and the last state R_s−1. The second step then corrects or updates this prior prediction into a posterior filtering distribution, p_s, which constrains p_s using the latest observation, I_s, according to Eq (5b). (5a) (5b)

Here is the state model from Eq (4b), is the observation model from Eq (1) and the constant of proportionality for Eq (5b) is simply a normalising factor.

Solving Eq (5) iteratively and simultaneously over the grid of leads to our novel real-time estimate of the time-varying effective reproduction number. We initialise this process with a uniform prior distribution over for p₁ and note that p_s and p_s are m element vectors that sum to 1, with i^th term corresponding to when . Eq (5) forms the first half of EpiFilter, is flexible and can be adapted to many related problems [14]. A key difference between EpiFilter and the EpiEstim-type methods [11, 13] is that the latter approximate the distributions and with and , respectively. Estimators based on these approximations can be suboptimal, especially when data (i.e., cases) are scarce.

Bayesian (backward) recursive smoothing

While Eq (5) provides a complete real-time solution to the filtering problem, it is necessarily limited at the starting edge of the incidence curve, where past data are sparse or unavailable. Further, because it does not update past estimates as new data accumulate, it cannot provide optimal retrospective estimates. Here we develop the second half of EpiFilter, which involves solving the optimal smoothing problem i.e., computing the smoothing posterior distribution , which provides maximally informed estimates of the instantaneous reproduction number R_s, given the complete incidence curve . To our knowledge, smoothing has not yet been explicitly considered in infectious disease epidemiology (either exactly or approximately).

We specialise the general methodology from [14, 25] to obtain the recursive smoother of Eq (6a). This equation uses the filtering distribution, and the predictive distributions , which we obtain from Eq (5). Our state model means that . (6a) (6b)

We realise Eq (6a) exactly by taking a forward-backward algorithmic approach (this is the backward pass whereas Eq (5) is the forward one). We solve this equation by noting that q_t = p_t and iterating backwards in time to obtain the first smoothing distribution q₁. The integrals become sums over our grid and distributions are m element vectors. Eq (6a) sequentially updates our earlier filtering solutions to include future data, forms the second half of EpiFilter and can also be approximately solved using particle smoothers [14].

This approach neatly links the filtering and smoothing distributions p_s and q_s. If we assume that the reverse-filtering distribution r_s is reasonably approximated by then we can also apply the two-filter smoothing solution of [25] to get Eq (6b). If either future () or past () incidence is uninformative then either r_s or p_s will reduce to the prior , leading to q_s ∝ p_s or q_s ∝ r_s, respectively. The end and beginning of the epidemic provide important examples of each of these scenarios. Consequently, Eq (6b) shows how smoothing connects EpiEstim and the WT methods, and explains why EpiFilter, which can be used for both real-time and retrospective inference, better overcomes edge-effects and periods of low data.

Further, the smoothed posterior q_s yields the conditional mean estimate , which is known to significantly improve on the MSE of the filtered equivalent (see previous section) [23]. While filtering provides the minimum MSE estimator of every instantaneous reproduction number R_s given past knowledge, smoothing provides the minimum given all knowledge. This relationship is formal, with filtered and smoothed MSE values mapping to the amount of mutual information that provides about [34]. Extracting the maximum information from the incidence curve should engender estimates that are more robust and statistically efficient in periods of low incidence. We summarise the EpiFilter algorithm in Fig 1.

While our main interest is on optimised and rigorous real-time and retrospective estimates of transmissibility, which are completely defined by the smoothing distribution q_s, we may also want to predict future incidence, for informing epidemic preparedness plans and for validating past R_s estimates [13, 35]. We compute the filtered one-step-ahead posterior predictive distribution as in Eq (7) (integrals are over ) [14]. (7)

We assume, as in [36], that . Replacing p_s with q_s yields the smoothed equivalent of Eq (7). We will use Eq (7) to compare EpiFilter against APEestim, which is the prediction-optimised version of EpiEstim, developed in [13]. Since predictions depend strongly on p_s or q_s, optimising these distributions can be important, for example, when forecasting second waves of infection.

Last, we comment on how we represent uncertainty in our estimates and predictions. Often we will provide 95% equal tailed Bayesian credible intervals. These are computed directly from the 2.5^th and 97.5^th quantiles of the relevant posterior distribution and used to assess performance statistics such as coverage of true values.

Improved estimation at low incidence

The reliable estimation of time-varying effective reproduction numbers, R_s, at low incidence, I_s, is a key challenge limiting our understanding of transmission [17]. Periods with small counts of new cases contain little information and so present necessary statistical difficulties [26]. Uncovering the transmissibility patterns in these periods can evidence intervention relaxation strategies, indicate potential elimination and provide insight into mechanisms of spread [28, 50]. Here we compare EpiFilter, which allows exact inference over a state grid , with EpiEstim and APEestim at these data-poor settings. We use EpiEstim with weekly and monthly windows. Weekly windows are the default recommendation in [11]. We include monthly ones as long windows can improve robustness at low incidence [28]. APEestim, developed in [13], optimises the window choice of EpiEstim to minimise one-step-ahead prediction errors. The assumptions and choices inherent in estimation methods become important and visible when data are scarce, and can bias inference or support spurious predictions [15].

We generate epidemics using the Poisson noise model of Eq (1) under the serial interval distribution of Ebola virus disease described in [37]. We examine three diverse scenarios in Fig 2. These describe (A) rapidly controlled epidemics (R_s step changes from 2 to 0.5 at s = 100), (B) small outbreaks with exponentially rising and falling R_s (change-point at s = 30 and rates of 0.02 and -0.008 per time unit) and (C) medium outbreaks that are initially controlled (R_s changes from 4 to 0.6 at s = 40) then resurge (R_s rebounds to 2 at s = 80) into large epidemics, before finally being suppressed (R_s = 0.2 from s = 150). These scenarios are similar to ones investigated in [12, 28] and describe various epidemic dynamics that culminate in elimination. We simulate 200 incidence curves for each scenario and apply APEestim (optimal window k*), EpiEstim (k = 7 and 31) and EpiFilter (state noise η = 0.1) to estimate R_s and sequentially predict I_s (given data up to s − 1) for each curve.

For the first two methods we compute mean one-step-ahead incidence predictions () as in [13] and instantaneous reproduction number estimates () from Eq (2) with τ(s) delimiting the window times used. We obtain smoothed EpiFilter estimates () and filtered predictions () from Eqs (6) and (7). We could use smoothed predictions (), which implicitly include data beyond s − 1, to assess fit but as we want to test model adequacy is more appropriate [38]. Our main focus is on real-time performance so we do not investigate the WT method. See [11, 12] for comparisons of the WT method and EpiEstim. Fig 2A–2C provide representative single runs (we only show k = 31 for EpiEstim as k = 7 is often similar to APEestim) while Fig 2D presents reproduction number MSE and one-step-ahead predicted incidence MSE (PMSE) distributions from the 200 runs.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Small or waning epidemics.

We compare reproduction number estimates ( or ) and one-step-ahead incidence predictions ( or ) from APEestim with optimal window k*, EpiEstim with window k and EpiFilter with state noise η. We simulate 200 epidemics with low daily case numbers or long tails (long sequences of zero cases) using the standard renewal model (Eq (1)) for three scenarios, representative examples of which are given in A-C. The true R_s and I_s are in black. All mean estimates or predictions are in red and blue with 95% credible intervals. APEestim and EpiEstim use a Gam(1, 2) prior distribution and EpiFilter a grid with m = 2000, R_min = 0.01 and R_max = 10. In D we provide statistics of the MSE of these estimates (relative to R_s) and the PMSE of these predictions (relative to I_s) for all 200 runs. We find that EpiFilter is more robust to small incidence (better uncertainty), whereas the other approaches can quickly decay to their prior distribution. It achieves significantly smaller MSE (2–10 fold reductions) and comparable PMSE to APEestim (which is optimised for prediction).

https://doi.org/10.1371/journal.pcbi.1009347.g002

In the Methods we showed that since EpiEstim-type methods group data over some window into the past, they revert to their prior distribution as the total cases in this window becomes small. During low incidence periods, e.g., when the epidemic is waning, using long windows makes sense [28]. However, this reduces predictive accuracy as fluctuations in R_s are underfit. This is the use-case for APEestim, which optimises for prediction. The consequences of these trade-offs are made clear in Fig 2A–2C where APEestim provides the best one-step-ahead I_s predictions but quickly reverts to its prior (seen as wide credible intervals) when cases are few. Longer window EpiEstim improves on this destabilisation but cannot track salient fluctuations in R_s. Fig 2D verifies these trends. APEstim has the smallest PMSE but often worse MSE than the larger k EpiEstim choices.

Interestingly, EpiFilter is able to jointly optimise both instantaneous reproduction number tracking and incidence predictions (computed via Eq (7)). In Fig 2A–2C we see that EpiFilter maintains stable and accurate estimates of R_s and only slowly reverts to its prior (which has the same support as that of EpiEstim and APEesim). However, unlike long window methods it does not sacrifice prediction fidelity. The improvement in MSE shown in Fig 2D is stark (the numerical reduction in MSE when compared to the next best method is on average at least 2-fold and often 10-fold). The PMSE, while larger than that of APEestim (which is optimised for predictions) is still good. These points are reinforced in Fig 3, which expands on the statistics of scenarios A-C.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Temporal statistics of small or waning epidemics.

We expand on the results from Fig 2D by decomposing the MSE and one-step-ahead PMSE statistics across the 200 simulated trajectories for every scenario in Fig 2. We do not consider the k = 31 EpiEstim example given its poor performance. We present the k = 7 case, which is the generally recommended EpiEstim setting. We observe that EpiFilter significantly improves on MSE throughout the epidemic trajectory (and not only in periods of low incidence) while maintaining comparable prediction accuracies. Coverage statistics for these scenarios, which are given in Fig C of the S1 Appendix, confirm that EpiFilter also consistently contains the true R_s and I_s values within its credible intervals.

https://doi.org/10.1371/journal.pcbi.1009347.g003

There the MSE and PMSE computed over the 200 replicate simulations are plotted with time. We observe a clear and significant reduction in MSE at most time points (both at low and large incidence) for EpiFilter, with similar PMSE performance, as compared to EpiEstim and APEestim. This confirms the widespread benefits of smoothing solutions. Moreover, the coverage of both the true R_s and I_s of EpiFilter (i.e., the probability that R_s or I_s is contained within estimated or predicted 95% equal tailed credible intervals) is more consistent than all other approaches. We present this in Fig C of the S1 Appendix. Thus, EpiFilter combines the advantages of APEestim and long-window EpiEstim, strongly minimises estimation error and, further, is able to reliably detect transmission change-points automatically.

Improved estimation between epidemic waves

Maintaining robust instantaneous reproduction number, R_s, estimation when incidence, I_s, becomes small is not just statistically important. Two possible outcomes may follow periods of small I_s: either the epidemic goes extinct (elimination occurs, as in the previous section), or an additional wave of infection surfaces (resurgence) e.g., due to imports or unmonitored local transmission. Predicting which outcome is likely, in real-time, is of global concern as countries aim to relax interventions during the ongoing COVID-19 pandemic, while also minimising the risks of further cases [17, 39]. As changes in instantaneous reproduction numbers signal variations in transmission and hence incidence, reliably identifying and inferring R_s trends in the trough preceding potential new peaks can be crucial for preparedness, providing evidence for timely and effective epidemic interventions [10].

Reliable estimation of R_s between epidemic waves depends on the prior assumptions of the inference method used and on how that method relies on those assumptions when data are scarce [40, 41]. Here we examine this dependence and investigate cases where resurgence follows a low-incidence period. As in the above section, we compare EpiFilter (η = 0.1) with APEestim (optimal window k*) and EpiEstim (weekly, k = 7, and monthly, k = 31, windows) over 200 simulated epidemics under the serial interval of Ebola virus. We explore scenarios depicting (A) epidemics that are initially controlled (R_s falls from 2.5 to 0.5 at s = 70) but which resurge just as quickly (R_s returns to 2.5 from s = 230), (B) periodic or seasonal transmission (R_s is sinusoidal with magnitude 1.3 ±1.2 and period of 120 time units) and (C) outbreaks with exponentially rising and then falling transmissibility (change-points at s = 40 and 190 and exponent rates 0.03, -0.015 and 0.02).

These examples are similar to some in [13] and describe diverse epidemics with multiple peaks and troughs. We provide representative runs of each scenario in A-C and collect MSE (for R_s) and PMSE (relative to I_s) distributions over the 200 runs of every scenario in Fig 4D. We also provide these statistics across time in Fig 5. We observe a similar pattern in performance among the methods as in our previous analyses of Figs 2 and 3. APEestim is best able to predict upcoming incidence and achieves the lowest PMSE as expected. While the monthly window EpiEstim is less useful in these cases (since prior reversion does not occur as often, though could be an issue for more extended troughs), the weekly window version loses predictive fidelity for minor MSE improvements.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Epidemics with multiple waves.

We compare reproduction number estimates ( or ) and one-step-ahead incidence predictions ( or ) from APEestim with optimal window k*, EpiEstim with window k and EpiFilter with state noise η. We simulate 200 epidemics with multiple waves of infection using the standard renewal model (Eq (1)) for three scenarios, representative examples of which are given in A-C. The true R_s and I_s are in black. All mean estimates or predictions are in red and blue with 95% equal tailed credible intervals. APEestim and EpiEstim use a Gam(1, 2) prior distribution and EpiFilter a grid with m = 2000, R_min = 0.01 and R_max = 10. In D we provide statistics of the MSE of these estimates (relative to R_s) and the PMSE of these predictions (relative to I_s) for all 200 runs. We find EpiFilter is best able to negotiate troughs between epidemic peaks and hence infer resurging infectious dynamics, achieving significantly smaller MSE (2–10 fold reductions) and comparable PMSE to APEestim (which is optimised for prediction).

https://doi.org/10.1371/journal.pcbi.1009347.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Temporal statistics of epidemics with multiple waves.

We expand on the results from Fig 4D by decomposing the MSE and one-step-ahead PMSE statistics across the 200 simulated trajectories for every scenario in Fig 4. We do not consider the k = 31 EpiEstim example given its poor performance. We present the k = 7 case, which is the generally recommended EpiEstim setting. We observe that EpiFilter significantly improves on MSE throughout the resurgent epidemic trajectory (whether incidence is small or large) while maintaining comparable prediction accuracies. Coverage statistics for these scenarios, which are given in Fig D of the S1 Appendix, confirm that EpiFilter consistently contains the true R_s and I_s values within its credible intervals.

https://doi.org/10.1371/journal.pcbi.1009347.g005

EpiFilter once again combines the advantages of the other approaches. For every scenario in Fig 4 it provides accurate tracking of changes in R_s with stable credible intervals and a MSE that is at least and sometimes even that of the next best method. The improvement in MSE is clear and maintained at almost every time point of each scenario, regardless of whether incidence is low or high. Concurrently, the incidence PMSE of EpiFilter rivals that of APEestim and it attains the most consistent coverage of the true R_s and I_s values, as shown in Fig D of the S1 Appendix. EpiFilter is therefore a proficient tool for inferring resurgence. We also find that the η = 0.1 parameter value seems to be an all-purpose heuristic, meaning that usage of EpiFilter can be simpler than EpiEstim and other window or change-point based methods. The improvements of EpiFilter likely result from its minimal assumptions (see Eq (4)) and its increased information extraction. We next test our method on empirical data.

COVID-19 in New Zealand and H1N1 influenza in the USA

The previous sections confirmed EpiFilter as a powerful inference and prediction tool, especially in data-poor conditions, using simulated epidemics. We now confront our method with empirical data from the 1918 H1N1 influenza pandemic in Baltimore (USA) [42] and the ongoing COVID-19 pandemic in New Zealand (up to 17 August 2020) [43]. The H1N1 dataset has been well-studied and so we first use this to benchmark EpiFilter. We clean this dataset by applying a 5-day moving average filter as recommended in [42]. Previous work [11] analysed this dataset with EpiEstim and found that sensible instantaneous reproduction number, R_s, estimates are obtained when a weekly window (k = 7) is applied. However, more recent work, using APEestim [13], showed that while k = 7 provides stable estimates for this epidemic, it is a poor predictor of the incidence data. Instead, an optimised window of 2 days (k* = 2) yields good predictions but the resulting R_s estimates are noisy.

We reproduce the instantaneous reproduction number estimates () and incidence predictions () from both studies in Fig 6 and compare them against EpiFilter with η = 0.1 ( and ). Top and middle rows of Fig 6 illustrate the aforementioned trade-off between estimate stability and prediction accuracy. The bottom row confirms the power of EpiFilter. Our R_s estimates are of comparable stability to those of EpiEstim at k = 7, yet our prediction fidelity matches that of APEestim. Our improved inference again benefits from using more information (i.e., the backward pass in Fig 1) and making less restrictive prior assumptions. We see the latter from the R_s credible intervals over 40 ≤ s ≤ 60. There EpiEstim seems overconfident, and this results in a rigid overestimation of incidence. However, EpiFilter mediates its estimate uncertainty to a level similar to APEestim.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. H1N1 influenza transmission in Baltimore (1918).

We compare APEestim (top), EpiEstim with recommended weekly window (middle) (both with Gam(1, 2) prior distribution) and EpiFilter (with m = 2000, η = 0.1, R_min = 0.01 and R_max = 10) on the H1N1 influenza dataset from [42]. We use a 5-day moving average filter, as in [42], to ameliorate known sampling biases. Estimates of reproduction numbers, R_s, and corresponding 95% equal tailed credible intervals are in red. One-step-ahead predictions of incidence, I_s, (with 95% credible intervals) are in blue with the actual incidence in black. We find that EpiFilter combines the benefits of APEestim and EpiEstim, achieving both good estimates and predictions.

https://doi.org/10.1371/journal.pcbi.1009347.g006

We explore COVID-19 transmission patterns in New Zealand using incidence data up to 17 August 2020 from [43]. New Zealand presents an insightful case study because officials combined swift lockdowns with intensive testing to achieve and sustain very low incidence levels that eventually led to local elimination of COVID-19 [44]. However, an upsurge in cases in early August inspired concerns about a second wave (which led to new interventions and is why we do not consider data beyond 17 August). Here we investigate the time-varying transmission in New Zealand to see if this uptick suggests that the epidemic was resurfacing in mid-August. We believe smoothing can confer important inferential advantages in exactly these types of low incidence scenarios.

We make the common assumptions that case under-reporting is constant [11], which seems reasonable given the intensive surveillance employed by New Zealand [45]. We ignore reporting delays, which are known to be small [46] and use the COVID-19 serial interval distribution from [47]. We do not explicitly distinguish imported from local cases in our analysis. The latter could bias our study [28, 39] but our focus is on demonstrating differences between filtering and smoothing on instantaneous reproduction number, R_s, trends and not on providing detailed R_s estimates during this period. We plot the results of our exploration in Fig 7.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. COVID-19 transmission in New Zealand.

We compute smoothed and filtered reproduction number estimates, (red) and (blue) respectively, from the COVID-19 incidence curve for New Zealand (available at [43]) in the left panels. We use EpiFilter with m = 2000, η = 0.1, R_min = 0.01 and R_max = 10 with a uniform prior distribution over the grid . The top of 7A shows conditional mean estimates and 95% credible intervals for (red) and (blue). Vertical lines indicate the start and end of lockdown, a major intervention that was employed to halt transmission. The additional ‘future’ information used in smoothing has a notable effect. The bottom of 7A provides smoothed one-step-ahead predictions (blue, with 95% credible intervals) of the actual reported cases I_s (black). The inset gives the estimated probability of R_s ≤ 1. We observe a clear trend of subcritical transmission that eventually seeds a second wave by August. In 7B we compare EpiFilter with EpiEstim (using weekly windows) and APEestim (both with Gam(1, 2) priors) with all left subfigures presenting R_s estimates and right ones providing filtered I_s predictions. We observe that both APEestim and EpiEstim lead to largely unusable estimates that mask transmission trends, in sharp contrast to EpiFilter.

https://doi.org/10.1371/journal.pcbi.1009347.g007

We apply EpiFilter and obtain filtered (, blue) and smoothed (, red) conditional mean estimates together with their 95% confidence intervals. These are in the left panels of Fig 7 and computed from Eqs (5) and (6) respectively. The times of lockdown and release are included for reference. Interestingly, we see a notable difference in the quality of inference between and . The former, as expected, is unreliable at the beginning of the incidence curve and features wider uncertainty and noisier trends. The smoothed , by using both forward and backward-looking data largely overcomes these issues and clarifies transmission dynamics. In the right panels we compare EpiFilter with EpiEstim (weekly windows) and APEestim. The difference is striking. Neither APEestim nor EpiEstim recovers a clear trend and both are appreciably worse than even the filtered estimate .

Our smoothed analysis suggests that R_s has resurged and supports the re-implementation of measures around 14 August. We recover suppression of the initial wave in April, likely associated with the implementation of key interventions, including lockdown [46]. Following this, we infer a prolonged period of subcritical transmission, where . Low and then zero incidence over this period strikingly destabilises both APEestim and EpiEstim and conceals the rise in R_s that occurs after July. EpiFilter signals this upsurge, which continually grows until August, where a second wave becomes likely. We also provide one-step-ahead predictions (which are from the smoothed hence the notation ) and their equal tailed 95% credible intervals against the reported incidence from [43]. These verify that reasonably describes the data (we also present when comparing methods). We find that around July, further supporting this resurgence hypothesis.

While the analyses of the H1N1 and COVID-19 data above illustrate the advantages of EpiFilter, these estimates can be further improved. We used case data and not actual infection times (which relates to the approximation of the generation time by the serial interval distribution) and did not account for introductions, case ascertainment fractions and reporting delays. For more practical analyses it may be necessary to first compensate for these biases to obtain the best possible incidence curve. Methods from [6, 12, 48] can be applied if relevant incubation period, reporting distributions and contact tracing data are available to diagnose and correct for these issues. The resulting pre-processed incidence curve can then be input to EpiFilter to obtain more realistic R_s estimates.