Improved estimation of time-varying reproduction numbers at low case incidence and between epidemic waves

INTRODUCTION

During an unfolding epidemic, one of the most commonly available and useful types of surveillance data is the daily (or weekly) number of new infected cases. This time-series of case counts, also known as the incidence curve, not only measures the epidemic size and burden, but also provides information about trends or changes in its transmissibility [1], [2]. These trends are captured by the time-varying effective reproduction number, denoted R_s at time s, which defines how the number of secondary cases generated per primary case varies across the duration of the outbreak [3]. When R_s > 1 we can expect and prepare for growing incidence, whereas R_s < 1 signifies that the epidemic is being controlled [4].

Inferring changes in R_s given an observed incidence curve is therefore crucial, both to understanding transmissibility and to forecasting upcoming case loads, especially in an ongoing epidemic, where it can help inform policymaking and intervention choices [1], [5]. Real-time and retrospective R_s estimates have been used to characterise rates and patterns of spread in various diseases such as malaria [6] and Ebola virus disease [7] and have already proven valuable across the COVID-19 pandemic, providing updating synopses of global transmission [8] or evidencing the impact of implemented control actions (e.g. lockdowns and social distancing) [9], [10].

Many of the studies that aim to infer R_s either apply the Wallinga-Teunis (WT) method [2] or the Cori et al method, known as EpiEstim [11]. Both methods take complementary viewpoints on how incidence data inform on R_s and hence have diverse use-cases. The WT method reconstructs the average number of new cases caused by infectious individuals at s and so requires incidence data beyond time s for its estimate, which is called the case reproduction number. It is better suited for retrospective analyses [12]. Alternatively, EpiEstim infers how past infections propagate to form the incidence observed at s, only requiring data prior to time s. EpiEstim is preferred for real-time investigations and its R_s estimate is termed the instantaneous reproduction number [3].

While both methods provide useful and important estimators of transmission, they are not perfect. Two main limitations exist. First, each suffers from data censoring or edge-effects [3]. Because the WT method is forward-looking i.e. depends on data later than s, its estimates are right censored when s is close to the last observed time point [12]. In contrast, EpiEstim looks backward in time and suffers edge-effects when s is near the first observed time point [11]. Estimates in the vicinity of the start and end of the incidence time-series are therefore unreliable under EpiEstim and the WT method, respectively.

Techniques have been proposed to limit this unreliability [5], [14], but the problem is intrinsic, and inevitable near the actual start and end times of an epidemic, where there is necessarily no or scarce data. This leads into the second key limitation: inference in periods of small incidence. This presents a fundamental challenge for any R_s estimation approach and effectively creates additional edge-effects. When few or no case counts are available to constrain inference, methods are largely driven by their inherent prior distributions and assumptions, which can result in misleading and unreliable estimates [11], [15]. Understanding how to best mediate the trade-off between prior assumptions and data when incidence is small is of both statistical and epidemiological importance.

Following a period of low incidence, two outcomes are possible: either the epidemic continues to exhibit small or zero case counts until it goes extinct or a resurgence in infections, also termed a second wave, occurs. Deciphering, in real time, which of these conditions is likely is a key challenge for infectious disease epidemiology given the limited information available at low incidence [16]. Better inference of transmission under these conditions is currently considered central to designing data-informed COVID-19 intervention exit or relaxation strategies [17]. With many countries now facing potential second waves in this pandemic, estimating fluctuations in transmission during suspected epidemic troughs could be essential to achieving sustained control [10].

Here we present and develop a novel method, termed EpiFilter, for reliably estimating R_s in real time, which ameliorates the above limitations. We take an engineering inspired approach and construct an exact, recursive and deterministic (i.e. EpiFilter produces the same output for fixed input data and requires no Monte Carlo steps) inference algorithm that is quick and easy to compute both across an unfolding outbreak and in retrospect. Our method solves what is called the smoothing problem in control engineering [13]. This means we compute R_s estimates that formally integrate both forward and backward looking information. This unifies the WT method and EpiEstim, and nullifies their edge-effect issues.

Further, EpiFilter only makes a minimal Markov assumption for R_s, which allows it to avoid the strong prior window size and change-point assumptions that existing methods may apply to infer shifts in transmission [11], [9]. Using simulated and empirical data, we show that EpiFilter accurately tracks changes in R_s and provides reliable one-step-ahead incidence predictions. Moreover, we find that EpiFilter is appreciably more robust than even optimised versions of EpiEstim [14], if performing real-time inference in periods of low incidence and between epidemic waves. We illustrate the practical utility of EpiFilter on the COVID-19 incidence curve of New Zealand, which exhibits a potential second wave.

Our method, which is outlined in Fig. 1, provides a straightforward yet formally optimal (in mean squared error [18]) solution to real-time and retrospective R_s estimation. Because it couples minimal prior assumptions with maximum information extraction, it more gracefully handles periods with scarce data. Hopefully our approach will serve as a useful inference tool for investigating the risk of second waves of infection in COVID-19 and other epidemics. An implementation of EpiFilter is available at https://github.com/kpzoo/EpiFilter and its mechanics are explored and validated in the Appendix.

• Download figure
• Open in new tab

Fig 1: EpiFilter algorithm.

We summarise the construction of EpiFilter and compare its use of information to that of EpiEstim and the WT method. Left panels show that EpiEstim and the WT method consider complementary parts of the incidence curve, (blue dots), to be informative about the reproduction number R_s, while EpiFilter makes use of the entire curve. Red windows indicate the portions of the curve that inform on R_s under a window of size k as in [11], while blue windows show the portions involved in constructing the main posterior distributions in this work: p_s, r_s and q_s. Right panels outline the main assumptions (the model box) and computations (the algorithm box) necessary for realising EpiFilter, which allow us to obtain the most informative (and minimum mse) posterior distribution q_s. See the main text for the specific equations employed in our implementation [13].

METHODS

Renewal models and inference problems

We consider an infectious disease epidemic observed over some period of time 1 ≤ s ≤ t. If the incidence or number of newly infected cases at time s is I_s then is known as the incidence curve of the epidemic. For convenience, we assume that incidence is available on a daily scale so that is a vector of t daily counts (weeks or months could be used instead). A common problem in infectious disease is the inference of the transmissibility of the epidemic given this incidence curve. The renewal model [1], [19] presents a general and popular framework for investigating this problem.

This model posits that epidemic transmissibility, summarised by the effective reproduction number, R_s, generates the observed incidence as in Eq. (1). There Pois indicates Poisson observation noise, means that both distributions are equal and .

In Eq. (1) R_s describes the branching of the epidemic i.e. it captures the number of secondary cases generated per effective primary one at s, while Λ_s, which is termed the total infectiousness, defines how many effective past cases are still infectious at s. 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 known for the epidemic of interest and well approximated by its serial interval distribution and focus on inferring the complete set of R_s values i.e. given with t as the last recorded time of the incidence curve. Estimating is important, not only for learning whether the epidemic is still growing (whenever R_s > 1) or under control (i.e. if R_s < 1) [4], but also for assessing how much existing interventions can be relaxed without risking a second wave of infection [10]. This latter problem can be particularly difficult (as we will show) since, under these conditions, incidence data are scarce.

We define three key inference problems, commonly studied in engineering [18], based on how information in is recruited to construct estimates. We represent these problems in terms of the posterior distribution they induce over . Estimates are then functions of these posteriors. The first is called causal filtering, where we sequentially calculate [20]. Causal means that estimates of R_s only depend on data up to time s. Solving this problem is fundamental to real-time inference [21]. The second problem is non-causal filtering and is the complement of the first. Here we calculate , which is important for retrospective or backward-looking estimates [22].

The last problem is termed smoothing and is our main interest. To solve it we must extract all the information possible to formulate [13]. Note that p_s, r_s and q_s depend on the choice 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 the observed incidence data [18]. These models encode our assumptions about the epidemic of interest and determine how estimates tradeoff assumptions against data. We will find this latter point important for determining performance when incidence data are scarce. We next explain how the most popular approaches for estimating , the WT method and EpiEstim fit within this engineering framework.

Inference methods and low incidence

We mostly detail EpiEstim as real-time R_s estimates are the key focus of this work. EpiEstim assumes that all reproduction numbers are the same over a window back in time, of size k. Consequently, only provides information about R_s and the causal filtering distribution follows as in Eq. (2) with a shape-scale gamma prior distribution applied i.e. Gam(a, c⁻¹) [11].

Here and are sums over the window of interest.

The resulting posterior mean estimate is then and the Fisher information that contains about R_s is proportional to λ_τ(s) [23]. The observation model is 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 [14]. As a result, EpiEstim somewhat incorporates a linear moving average state space model, and assumes that p_s ≈ p_τ(s) by deeming to be uninformative. 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 and computes an estimate of R_s based on a forward-looking window [2]. It also uses the observation model of Eq. (1) and has an implicit moving average state model [3]. We abuse notation and define in this case as a window into the future of size k. The WT method assumes that and approximately solves the non-causal filtering problem. We illustrate the information windows employed by the WT method and EpiEstim and the complete causal and non-causal filtering windows in Fig. 1. The goodness of p_τ(s) and r_τ(s) as measures of p_s and r_s will depend on k and the appropriateness of the state model [14].

While EpiEstim and the WT methods are powerful tools for inferring R_s 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 p₁ = p_τ(1) = P(R₁ | I₁) and r_t = r_τ(t) = P(R_t | I_t), 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 and prolonged [23], [17]. In these periods information is scarce and the quality of estimates will depend on how well the method of choice mediates between the little available data and its inherent assumptions [15]. We illustrate this effect using EpiEstim by considering an epidemic with a sequence of n zero incidence days. This sequence is realistic for an epidemic in its tail phase, and can precede either elimination (i.e. the epidemic goes extinct) or resurgence (i.e. a second wave of infection emerges) [24]. As n increases both i_τ(s) and λ_τ(s) shrink, though i_τ(s) decays more rapidly. For every n ≥ k, i_τ(s) = 0, meaning that the shape of p_τ(s) is completely determined by the prior P(R_s).

As n increases further λ_τ(s) becomes negligible and we exactly recover our prior distribution (see Eq. (2)). This problem is acute if the window size, k, is small. Previous studies, which formally optimise k to maximise estimate reliability, found that small k is needed when inferring sharp changes in transmissibility (e.g. due to lockdowns) [14]. Moreover, as λ_τ(s) falls, the Fisher information available about R_s from shrinks, leading to an inflation in estimate uncertainty and a loss of statistical identifiability [25], [23]. An analogous effect occurs in the WT method if there is low incidence across its forward-looking window [12].

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 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 . This follows as q₁ = r₁ and q_t = p_t. 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 develop a method, termed EpiFilter, to realise these improvements in the next section.

Bayesian (forward) recursive filtering

We reformulate the causal inference problem of estimating R_s from as an optimal Markov state filtering problem. Filtering is the term used in engineering for a general class of estimation problems aimed at optimally (usually in a mean squared error (mse) sense) inferring some hidden state causally and in real time from noisy observations [18], [13]. Given 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 [20], which EpiEstim approximates. The conditional mean estimate leads to the causal minimum mse of [20].

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

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 [30].

We formulate our filter under two very mild assumptions. First, we define some closed space, , over which R_s is valid i.e. for a given resolution m, extrema R_min and R_max, and grid size then . This means R_s must take a discrete value in , the i^th element of which is . We formalise this in Eq. (4a).

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). Note that other approaches, which depend on MCMC or related sampling methods, all implicitly assume some sort of discretisation [27].

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 η is 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 [20] and epidemiology [31] 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 R_s−1. This allows variation but ensures 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 and are summarised in Fig. 1. 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 [30], [13], [22], [18], which are valid for various types of observation and state models, for our renewal model inference problem. The proof of the equations we employ can be found in these works.

Recursive filtering involves two steps: prediction and correction. The first, given in Eq. (5a), constructs a sequential prior predictive distribution, , which is informed by past data and the last state R_s−1 The second step then corrects or updates this prior into the posterior distribution, p_s, which constrains p_s using the latest observation, I_s, as in Eq. (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 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 [13]. A key difference between EpiFilter and the EpiEstim-type methods [11], [14] is that the latter approximate and with and P(R_s), 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 causal filtering problem, it is necessarily limited at the starting edge of the incidence curve, where past data are scarce 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 and hence computing q_s. To our knowledge, smoothing has not yet been explicitly considered in infectious disease epidemiology.

We specialise the general methodology from [13], [22] to obtain the recursive smoother of Eq. (6a).

We can realise Eq. (6a) exactly by taking a forward-backward algorithmic approach with p_s and p_s+1 first computed from Eq. (5). Remaining terms emerge from the state space model and by noting that q_t = p_t and iterating backwards. This approach neatly links p_s and q_s. If we assume that r_s is reasonably approximated by then we can also apply the two-filter smoothing solution of [22] to get Eq. (6b). This shows how smoothing connects EpiEstim and the WT methods, and explains why EpiFilter can be used for both real-time and retrospective inference. We summarise the main ingredients of the EpiFilter algorithm in Fig. 1.

The smoothed posterior q_s yields the conditional mean estimate , which is known to significantly improve on the mse of the filtered equivalent [20]. While filtering provides the minimum mse given causal 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 [32]. Extracting the maximum information from should engender estimates that are more robust to periods of low incidence.

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

We assume, as in [34]:. 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 [14]. Since predictions depend strongly on p_s or q_s, optimising these distributions can be crucial, for example, to forecasting second waves of infection.

RESULTS

Improved estimation at low incidence

The estimation of time-varying effective reproduction numbers when incidence is low is seen as 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 [23]. Here we compare EpiFilter, which allows exact inference but conditioned on , with EpiEstim and APEestim, which minimises the prediction error of EpiEstim by optimising its window size, at these data-poor settings [11], [14]. We use EpiEstim with a long window, as this is known to improve robustness at low incidence [24]. 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 simulate epidemics using Eq. (1)) under the serial interval distribution of Ebola virus disease available from [35]. We examine several diverse incidence trajectories (e.g. outbreaks with small sizes, large peaks or multiple modes) that all eventually decline to periods of few or no cases. We apply APEestim with optimal window , EpiEstim with a long window k and EpiFilter with state noise η (see Eq. (4b)) to estimate R_s and causally predict I_s for each epidemic. For the first two methods we compute mean predictions () as in [14] and estimates () from Eq. (2) with τ(s) delimiting the window times used. We obtain smoothed EpiFilter estimates () and causal predictions () from Eq. (6) and Eq. (7).

Our main results are in Fig. 2, which exposes how these three approaches degrade as data diminish. Each quadrant of this figure examines a different epidemic scenario. In all plots the true R_s and I_s are dashed black and dotted grey respectively, estimates are in red and predictions are blue. The central lines are conditional means and the shaded regions are their 95% confidence intervals. As our main focus is on real-time performance we do not investigate the WT method. See [11], [12] for comparisons of the WT method and EpiEstim. Our reasons for comparing to EpiEstim with a long window and APEestim follow from the Methods.

• Download figure
• Open in new tab

Figure 2: Small or dying epidemics.

We compare reproduction number estimates () and causal one-stepahead incidence predictions () from APEestim with optimised window (top panels), EpiEstim with long window k (middle panels) and EpiFilter with state noise η (bottom panels). We simulate epidemics with low daily case numbers or long tails (long sequences of zero cases) using the standard renewal model (Eq. (1)). The true R_s is in black (dashed, left panels) and its corresponding I_s in grey (dots, right panels). All mean estimates or predictions are in red with 95% confidence 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. We find that EpiFilter is more robust to small incidence (better uncertainty), whereas the other approaches can quickly decay to their prior distribution.

There we showed that since EpiEstim-type methods group data over a 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 dying or in its early growth phase, using long windows makes sense [24]. However, this reduces predictive accuracy as fluctuations in R_s are underfit.

This is the use-case for APEestim, which optimises for prediction. The top and middle panels of each quadrant of Fig. 2 illustrate the clear consequences of these tradeoffs. APEestim often selects shorter windows to better fit R_s and sequentially predict I_s. Each predicted I_s is a one-step-ahead prediction depending on the past data up to s − 1 as in Eq. (7) (we have abused notation).

While shorter windows improve estimates and predictions the reversion to the prior distribution (seen as wide confidence intervals in Fig. 2) when I_s is low is acute. Long windows resolve this latter difficulty, showing a slower rise to the prior but at the expense of getting other parts of the epidemic trajectories incorrect. Interestingly, in the bottom panels of the quadrants of Fig. 2 we see that EpiFilter has superior performance for these types of epidemics. Not only does it give good estimates and predictions throughout the trajectory of the epidemic but its rise in uncertainty as the epidemic dies is slower and more controlled. It therefore combines the advantages of APEestim and long-window EpiEstim.

Improved estimation between epidemic waves

Maintaining robust R_s estimation at low incidence is not just statistically important. Two possible outcomes may follow periods of small I_s: either the epidemic goes extinct, or a second wave of infection (also known as resurgence) occurs e.g. due to imports or unmonitored local transmission. Predicting which outcome is likely, rapidly and in real-time, is of major global concern as countries aim to relax interventions during the ongoing COVID-19 pandemic, while also minimising the risks of second waves [17]. As changes in reproduction numbers predate and signal corresponding variations in incidence, reliably identifying and inferring R_s trends in the trough preceding potential new peaks is crucial for preparedness and both timely and effective epidemic control [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 [36], [37]. Here we examine this dependence and exactly investigate this scenario, where after a low incidence period resurgence occurs. As in the previous section, we test APEestim with optimal window (top panels), EpiEstim with a long window k (middle) and EpiFilter with state noise η (bottom panels) on four diverse epidemics, which in this case all feature multiple waves or early hints of upcoming resurgence. Simulation parameters are unchanged from the last section.

Our main results are in Fig. 3 with each quadrant specifying a different epidemic trajectory. The true R_s and I_s are in dashed-black and dotted-grey respectively.

• Download figure
• Open in new tab

Figure 3: Epidemics with multiple waves.

We compare reproduction number estimates () and causal one-step-ahead incidence predictions () from APEestim with optimised window (top panels), EpiEstim with long window k (middle panels) and EpiFilter with state noise η (bottom panels). We simulate epidemics with multiple waves of infection using the standard renewal model (Eq. (1)). The true R_s is in black (dashed, left panels) and its corresponding I_s in grey (dots, right panels). All mean estimates or predictions are in red with 95% confidence 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. We find that due to the improved robustness of EpiFilter to low incidence regions, it is best able to negotiate troughs between epidemic peaks and hence infer re-emerging infectious dynamics.

All estimates ( for EpiEstim and APEestim and for Epifilter) are in red with corresponding causal predictions ( and ) in blue. Plots provide conditional means with 95% confidence intervals. For all second and additional wave epidemic examples we observe a similar tradeoff in performance between APEestim and long-window EpiEstim as in our previous analysis of Fig. 2. APEestim is better able to predict upcoming I_s terms and follow the overall trend in R_s than EpiEstim.

However, the long window of EpiEstim makes it more stable during the trough between waves, which can be advantageous for understanding transmission there.

EpiFilter once again combines the advantages of the other approaches. For every scenario in Fig. 3 it provides accurate tracking of changes in R_s with stable confidence intervals. Concurrently, its predictive performance rivals that of APEestim. EpiFilter is therefore a powerful tool when dealing with second-wave scenarios. We also find that the η = 0.1 parameter value seems to be an all-purpose heuristic, meaning that usage of EpiFilter can be somewhat simpler than APEestim and EpiEstim. The superior capability of EpiFilter in these scenarios likely results from its minimalistic assumptions (see Eq. (4)) and its increased information extraction. We next test our method on empirical incidence data.

COVID-19 in New Zealand and H1N1 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) [38] and the ongoing COVID-19 pandemic in New Zealand (up to 17 Aug 2020) [39]. The H1N1 dataset has been well-studied and so we first use this to benchmark EpiFilter. We clean this data by applying a 5-day moving average filter as recommended in [38]. Previous work [11] analysed this data with EpiEstim and found that sensible R_s estimates result when a weekly window (k = 7) is applied.

A recent study, which re-examined this epidemic with APEestim [14], shows that, while EpiEstim with k = 7 provides stable estimates for this epidemic, it is a poor causal predictor of the incidence data. Instead, an optimised window of 2 days () yields good predictions but the resulting R_s estimates are noisy. We reproduce the estimates () and predictions () from both studies in Fig. 4 and compare them against EpiFilter with η = 0.1 ( and ). The cleaned H1N1 incidence is in dotted grey, the critical R_s = 1 line in dashed black and all R_s estimates and I_s predictions (with 95% confidence intervals) are in red and blue respectively.

• Download figure
• Open in new tab

Figure 4: H1N1 influenza transmission estimates for Baltimore (1918).

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

Top and middle rows of Fig. 4 show the mentioned 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 that of EpiEstim at k = 7, yet our prediction fidelity matches 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 confidence intervals over 40 ≤ s ≤ 60.

There EpiEstim seems overconfident, and this results in a rigid overestimation of incidence. However, EpiFilter mediates its confidence to a level similar to APEestim.

We explore COVID-19 transmission patterns in New Zealand using incidence data up to 17 Aug 2020 from [39]. New Zealand presents an insightful case study because officials combined swift lockdowns with intensive testing to achieve and sustain very low incidence levels that some believed could engender local elimination of COVID-19 [40]. However, an upsurge in cases in early Aug inspired concerns about a second wave (which led to new interventions and is why we do not consider data beyond 17 Aug). Here we investigate time-varying transmission in New Zealand to see if this uptick suggests that the epidemic was resurfacing in mid Aug. 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 and that reporting delays are negligible [11]. While this makes our investigation somewhat naive, we consider this to be reasonable given the intense surveillance that New Zealand employed [41]. We also use the serial interval distribution from [42] and do not explicitly distinguish imported from local cases in our analysis. The latter could bias our study [24] but our focus is on demonstrating differences between filtering and smoothing on R_s trends instead of specifying the precise value of R_s during this period. We plot the results of our exploration in Fig. 5. More involved analyses can be performed by pre-processing the data for known delays or case ascertainment fractions as in [12], [6].

• Download figure
• Open in new tab

Figure 5: COVID-19 transmission in New Zealand.

We compute smoothed and filtered reproduction number estimates, and respectively, from the COVID-19 incidence curve for New Zealand (available at [39]). We use EpiFilter with m = 2000, η = 0.1, R_min = 0.01 and R_max = 10 with a uniform prior over the grid . The top panel shows conditional mean estimates and 95% confidence intervals for (red) and (grey). 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 panel provides smoothed one-step-ahead predictions (blue, with 95% confidence) of the actual reported cases I_s (grey). The inset gives the probability of ≤ 1 and lockdown times.

We apply EpiFilter and obtain causally filtered (, grey) and smoothed (, red) conditional mean estimates together with their 95% confidence intervals. These are in the top panel of Fig. 5 and computed from Eq. (5) and Eq. (6) respectively. The times of lockdown and release are the grey vertical dashed lines. 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 unclear trends. The smoothed estimate, by using both forward and backward-looking data largely overcomes these issues and clarifies transmission trends.

Our smoothed analysis suggests that R_s has resurged and supports the reimplementation of measures around 14 Aug. In the bottom panel of Fig. 5 we provide predictions (which are from the smoothed hence the notation ) and their 95% confidence intervals in blue against the reported incidence from [39] in dotted grey. These confirm that reasonably describes the data. In the inset we show that the P(R_s ≤ 1) also supports the resurgence hypothesis. Last, we note that transmission decreased considerably (and rise in P(R_s ≤ 1)) across the lockdown period.

DISCUSSION

Estimating time-varying trends in effective reproduction number R_s reliably and in real-time is an important and popular problem in infectious disease epidemiology

[5]. As the COVID-19 pandemic has continued to unfold, the interest in solving this problem has only elevated [8]. Initially, the focus was on characterising how R_s might respond to interventions such as lockdowns and social distancing [10]. However, as COVID-19 has progressed and countries have entered the controlled phase, concerns have shifted to trying to understand how existing interventions might be relaxed with minimum risk.

The literature on intervention exit strategies is, however, still in development, and several challenges remain to modelling the transmission behaviour of epidemics following prolonged periods of low incidence [17], [24]. While changes in R_s over this period are not the only analytic for assessing this risk, they do provide a key real-time diagnostic since upticks in R_s generally precede corresponding rises in incidence. Unfortunately, current approaches to estimating R_s under these conditions tend to be underpowered or prior-constrained [11], [23].

Here we re-examined existing methodology for inferring R_s from an engineering perspective. We observed that two of the most useful and popular inference approaches, EpiEstim [11] and the WT method [2], only capitalise on a portion of the data available, deeming either upcoming or past incidence informative (see Fig. 1). This informative portion is directly controlled by prior assumptions on the speed of possible R_s changes, which are often characterised by a window of size k. Other methods are also known to apply similarly strong change-point assumptions on R_s [23], [9]. When data are scarce these assumptions can control inference.

In control engineering a common problem, known as filtering, involves optimally (in a mse sense) estimating hidden Markov states, causally in real-time, from noisy and uncertain observations [18]. A related problem termed smoothing provides accompanying and optimal retrospective inference [13]. By reinterpreting R_s as a Markov state (Eq. (4)) observed through a noisy renewal process (Eq. (1)) and defining R_s on a predetermined grid R, we were able to construct exact filtering (Eq. (5)) and smoothing (Eq. (6)) solutions. This led to EpiFilter, which is our central contribution.

Generally, filtering and smoothing can be involved and require sophisticated sequential Monte Carlo techniques [30]. However, because our system is low dimensional and as we make only minimal assumptions about R_s, modelling it as a simple diffusion, we were able to solve these problems exactly and without complex sampling algorithms [27]. Our solutions are computationally simple and deterministic (i.e. precisely reproducible given the same data). While EpiFilter is limited by one free parameter η, our η = 0.1 heuristic seems statistically justified by its good causal one-step-ahead predictive performance (see Appendix and Fig. 2 and Fig. 3) [14].

Importantly, EpiFilter is able to look both forward and backward through the incidence data, and so maximise the information extracted [32]. This property means it combines the advantages of both EpiEstim and the WT method (see Fig. 1) and largely negates their edge-effect issues [12], offering both real-time and retrospective R_s estimations that are robust. We illustrated the advantages of EpiFilter by comparing it to EpiEstim and APEestim on many simulated examples with diverse periods of low incidence and epidemic resurgences (Fig. 2 and Fig. 3). Interestingly, we found EpiFilter able to track arbitrary changes in transmission without requiring any preset change-points or window sizes or any tuning of η.

Moreover, EpiFilter was notably better at negotiating periods of low incidence, offering a graceful degradation to its prior without sacrificing predictive accuracy. When incidence is low, it is usually beneficial to use long windows with EpiEstim [24]. This keeps R_s estimates reasonably stable. However, it often leads to poor predictions [14]. APEestim, which optimises window size for prediction fidelity, showed that in many of the simulated scenarios short windows are necessary for describing transmission patterns. Consequently, we have a tradeoff between estimate robustness and prediction accuracy.

We found that EpiFilter overcomes this tradeoff, concurrently achieving good estimates and predictions. We confirmed these advantages on empirical data from the H1N1 pandemic of 1918 (see Fig. 4). The ability to integrate the benefits of long window EpiEstim and APEestim make EpiFilter particularly useful for investigating resurgence after a period of low incidence, as it is better able to forewarn of increasing case numbers and can more speedily infer upticks in transmission (see Fig. 3). Capitalising on these properties, we performed a naive exploratory analysis of COVID-19 in New Zealand and found evidence of a potential second wave of infection after a prolonged period of control (see Fig. 5).

Balancing the assumptions inherent to a model against the data it is applied on, to produce reliable inference is a non-trivial problem that is still under active investigation in several fields [15], [36], [37]. EpiFilter, by maximising the information extracted from the incidence data and minimising its state space model assumptions, appears to strike this balance. Consequently, it performs strongly on a wide range of problems, including those involving sparse data, where other methods might struggle. Given its demonstrated advantages, straightforward formulation and theoretical underpinning, we hope that EpiFilter will be useful as a diagnostic tool for reliably warning about second waves of infection as COVID-19 unfolds.