Introduction

The effective reproductive number, denoted as R_e or R_t, is the expected number of new infections caused by an infectious individual in a population where some individuals may no longer be susceptible. Estimates of R_t are used to assess how changes in policy, population immunity, and other factors have affected transmission at specific points in time [1–5]. The effective reproductive number can also be used to monitor near real-time changes in transmission [6–11]. For both purposes, estimates need to be accurate and correctly represent uncertainty, and for near real-time monitoring, they also need to be timely.

We consider 2 potential forms of bias in R_t estimates, systematic over- or underestimation and temporal inaccuracy. Misspecification of the generation interval is a large potential source of over- or underestimation, and we find that R_t estimates are most prone to this kind of bias when the true value is substantially greater or less than 1. This situation might arise at the beginning of the Coronavirus Disease 2019 (COVID-19) pandemic (when R_t is relatively high) or after particularly effective interventions (when it might be low). Over- or underestimation would have particularly strong practical consequences near the control threshold of R_t = 1, but the biases we observe are smallest in absolute terms in this range.

Another challenge is that depending on the methods used, R_t estimates may be leading or lagging indicators of the true value [4,12], even measuring transmission events that occurred days or weeks ago if the data are not properly adjusted. Temporal inaccuracy in R_t estimation is particularly concerning when trying to infer how changes in behavior have affected transmission [1–5]. Temporal inaccuracy, for instance, in the estimation of the date on which R_t falls below 1, is a focus of this Perspective. We find that it has several possible causes and can be difficult to avoid.

This Perspective focuses on the 3 main empirical methods to estimate R_t [13–16]. As an alternative to the methods reviewed here, it is possible to infer changes in transmission using a dynamical model (e.g., [3,17–19]). The accuracy and timeliness of R_t estimates obtained in this way should be assessed on a case-by-case basis, giving sensitivity to model structure and data availability.

We use synthetic data to compare the accuracy of 3 common empirical methods with the estimate R_t, first under ideal conditions, in the absence of parametric uncertainty, and with all infections observed at the moment they occur. This idealized analysis is intended to illustrate the inputs needed to estimate R_t accurately, to highlight the intrinsic differences between the methods, and to examine specific causes of bias and temporal inaccuracy 1 by 1. However, we emphasize that our idealized analyses overestimate the potential accuracy of R_t estimates obtained from real-world data, even if best practices are followed. The results show that the method of Cori and colleagues [15] is best for near real-time estimation of R_t. For retrospective analysis, the methods of Cori and colleagues or of Wallinga and Teunis may be appropriate, depending on the aims.

Later, we add realism and address practical considerations for working with imperfect data. These analyses emphasize potential errors introduced by uncertainty in the intrinsic generation interval and imperfect case observation, the need to adjust for delays in case observation and right truncation, and the need to choose an appropriate smoothing window given the sample size. Finally, we emphasize that most off-the-shelf tools leave it up to the user to account for these 5 sources of uncertainty when calculating confidence intervals. Failure to propagate uncertainty in R_t estimates can lead to overinterpretation of the results and could falsely imply that confidence or credible intervals have crossed the critical threshold.

Synthetic data

We used synthetic data to compare 3 common R_t estimation methods. Synthetic data were generated from a deterministic or stochastic SEIR model in which the transmission rate changes abruptly. Results were similar whether data were generated using a deterministic or stochastic model. For simplicity, we show deterministic outputs throughout the document, except in the section on smoothing windows, where stochasticity is a conceptual focus.

In our model, all infections are locally transmitted, but all 3 of the methods we test can incorporate cases arising from importations or zoonotic spillover [13,14,16]. Estimates of R_t are likely to be inaccurate if a large proportion of cases involve transmission outside the population. This situation could arise when transmission is low (e.g., at the beginning or end of an epidemic) or when R_t is defined for a population that is connected to others via migration.

A synthetic time series of new infections (observed at the S→E transition) was input into the R_t estimation methods of Wallinga and Teunis, Cori and colleagues, and Bettencourt and Ribeiro [13–15]. Following the published methods, we also tested the Wallinga and Teunis estimator using a synthetic time series of symptom onset events, extracted daily from the E→I transition. In the synthetic data, the generation interval followed a gamma distribution with shape 2 and rate of , which is the sum of exponentially distributed residence times in compartments E and I, each with a mean of 4 days [20]. The generation interval of COVID-19 can also be approximated by a gamma distribution [22,23]. The methods of Cori and colleagues and of Wallinga and Teunis can accommodate any positive, discrete generation interval distribution, including the discretized gamma distribution used to generate the synthetic data [13,15,21]. However, the method of Bettencourt and Ribeiro implicitly assumes the generation interval follows an exponential distribution, as in a susceptible-infectious-recovered (SIR) model, which does not include a latent period [14,20]. Thus, the method could not match the assumptions of the synthetic data.

In the synthetic data, R₀ was set to 2.0 initially, then to 0.8 and 1.15, to simulate the adoption and later the partial lifting of public health interventions. To mimic estimation in real time, we truncated the time series at t = 150, before the end of the epidemic. Estimates from the methods of Wallinga and Teunis and Cori and colleagues were obtained using the R package EpiEstim [21]. Estimates based on the method of Bettencourt and Ribeiro were obtained by translating code from references [6,24] into the Stan language [25]. We initially assumed all infections were observed. Unless otherwise noted, the smoothing window was set to 1 day (effectively, estimates were not smoothed). To mimic the timescale of observations, we used daily time steps when generating synthetic data and performing analyses.

Comparison of common methods

The effective reproductive number at time t can be defined in 2 ways: as the instantaneous reproductive number or as the case reproductive number [15,26]. The instantaneous reproductive number measures transmission at a specific point in time, whereas the case reproductive number measures transmission by a specific cohort of individuals (Fig 1). (A cohort is a group of individuals with the same date of infection or the same date of symptom onset.) The case reproductive number is useful for retrospective analyses of how individuals infected at different time points contributed to the spread. It is a more natural choice for analyses that consider heterogeneity among individuals. For example, the case reproductive number of Wallinga and Teunis can be adapted to incorporate data on observed transmission chains [13,22,27] or to produce age-structured R_t estimates, given an age-structured contact matrix [28]. The instantaneous reproductive number is more appropriate for analyses estimating the reproductive number of the infected population on specific dates, especially when aiming to study how interventions or other extrinsic factors have affected transmission at a given point in time.

Download:

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

Fig 1. Instantaneous reproductive number as estimated by the method of Cori et al. vs. cohort reproductive number estimated by Wallinga and Teunis.

For each definition of R_t, arrows show the times at which infectors (upwards) and their infectees (downwards) appear in the data. Curves show the generation interval distribution (A, B), or serial interval distribution (C). (A) The instantaneous reproductive number quantifies the number of new infections incident at a single point in time (t_i, blue arrow), relative to the number of infections in the previous generation (green arrows) and their current infectiousness (green curve). The methods of Cori et al. and of Bettencourt and Ribeiro estimate the case reproductive number. This figure illustrates the Cori method. (B-C) The case reproductive number is defined as the average number of new infections that an individual who becomes infected on day t_i (green arrows in B) or symptomatic on day t_s (yellow arrows in C) will eventually go on to cause (blue downward arrows in B and C). The first definition applies when estimating the case reproductive number using inferred times of infection, and the second applies when using data on times of symptom onset. The method of Wallinga and Teunis estimates the case reproductive number.

https://doi.org/10.1371/journal.pcbi.1008409.g001

More formally, the instantaneous reproductive number is defined as the expected number of secondary infections occurring at time t, divided by the number of infected individuals, each scaled by their relative infectiousness at time t (an individual’s relative infectiousness is based on the generation interval and time since infection) [15,26]. The instantaneous reproductive number can be calculated exactly for a compartment model (SIR or SEIR) as follows, where β(t) is the time-varying transmission rate, S(t) the fraction of the population that is susceptible, and D the mean duration of infectiousness: (1) The methods of Cori and colleagues [15,16] and methods adapted from Bettencourt and Ribeiro [6,14,24] estimate the instantaneous reproductive number from observations. We tested their accuracy under idealized conditions, assuming perfect knowledge of the generation interval and delay distributions used to generate the synthetic data (Fig 2 and S1 Fig).

Download:

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

Fig 2. Accuracy of R_t estimation methods given ideal, synthetic data.

Solid and dashed black lines show the instantaneous and case reproductive numbers, respectively, calculated from synthetic data. Colored lines show estimates and confidence or credible intervals. To mimic an epidemic progressing in real time, the time series of infections or symptom onset events up to t = 150 was input into each estimation method (inset). Terminating the time series while R_t is falling or rising produces similar results as in S1 Fig. (A) By assuming an SIR model (rather than SEIR, the source of the synthetic data), the method of Bettencourt and Ribeiro systematically underestimates R_t when the true value is substantially higher than 1. The method is also biased as transmission shifts. (B) The Cori method accurately measures the instantaneous reproductive number. (C) The Wallinga and Teunis method estimates the cohort reproductive number, which incorporates future changes in transmission. Thus, the method produces R_t estimates that lead the instantaneous effective reproductive number and becomes unreliable for real-time estimation at the end of the observed time series without adjustment for right truncation [4,29]. In panels A and B, the colored lines show the posterior mean and the shaded region the 95% credible interval. In panel C, the colored line shows the maximum likelihood estimate and the shaded region the 95% confidence interval.

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

The method of Cori and colleagues estimates R_t as (2) where I_t is the number of infection incidents on day t, and w_s is the generation interval, or the probability that s days separate the moment of infection in an index case and a daughter case [15]. Conceptually, this estimator describes the number of new infections incident on day t relative to the number (I_t−s) and current infectiousness (w_s) of individuals who became infected s days in the past and who many now be shedding virus.

The only parametric assumption required by this method is the form of the generation interval. The standard assumption is that w_s follows a discretized gamma distribution [15,21], but the estimator accepts any parametric or empirical discrete distribution with support on positive values (the same is true of the Wallinga and Teunis method [13]). Thus, when testing the method on COVID-19-like epidemic processes, we could specify the gamma-distributed generation interval of the synthetic data perfectly. When tested under idealized conditions, we found that the method of Cori and colleagues accurately estimated R_t, even tracking abrupt changes (Fig 2). This is the method we recommend for estimation of the instantaneous reproductive number.

Bettencourt and Ribeiro [14] derive an approximate relationship between the R_t and the exponential growth rate of the epidemic, where g is the mean generation time: (3)

Under the assumption that R_t evolves through time as a Gaussian process, Eq 3 facilitates efficient Bayesian R_t estimation [24]. The key disadvantage of this method is that Eq 3 is derived from an SIR model and therefore implicitly assumes that the generation interval follows an exponential distribution [20], whereas empirically, generation intervals can be heavy or light tailed, including those for COVID-19 [22,23,30]. Because Eq 3 misestimates variability in the gamma-distributed generation interval of our SEIR-type synthetic data, we find that this method produces biased R_t estimates, especially when R_t is substantially higher than 1 (Fig 2A). These biases are consistent with established theory, in which the coefficient of variation of the generation interval modulates the relationship between exponential growth rate and the reproductive number of an epidemic [20].

In its current form, we do not recommend using the method of Bettencourt and Ribeiro, given that unrealistic structural assumptions lead to bias. However, a generalized version capable of accommodating more realistic generation intervals, which implicitly involves different assumptions about the underlying epidemic process [20,31], could provide several advantages. When implemented as a Gaussian process [24], we found that the Bettencourt and Ribeiro method was computationally efficient. Also, because it penalizes large jumps in R_t across consecutive time steps, it returns smoother estimates than the method of Cori and colleagues, which is advantageous if unmodeled reporting effects, rather than bursts in transmission, are the dominant cause of variability in daily observations.

Finally, the case or cohort reproductive number is the expected number of secondary infections that an individual who becomes infected at time t will eventually cause as they progress through their infection [15,20,26] (Fig 1B and 1C). The case reproductive number can be calculated exactly at time t within the synthetic data as the convolution of the generation interval distribution w(·) and the instantaneous reproductive number described in Eq 1 [20], (4)

The method of Wallinga and Teunis [13] estimates the case reproductive number from observations. The first step is to estimate the likelihood that case j (infected at time t_j) infected case i, relative to the likelihood that any other case in the data infected case i: (5) Then, the individual reproductive number of case j is defined as (6) The case reproductive number at time t is defined as the expected value of R_j for all individuals infected at time t. An assumption common to the 3 tested methods is that all infections are observed. Below we discuss the consequences of partial observation, including that confidence or credible intervals around R_t estimates do not account for uncertainty from partial observation.

Practically speaking, there are several important differences between the case reproductive number (estimated by Wallinga and Teunis) and the instantaneous reproductive number (estimated by Cori and colleagues or Bettencourt and Ribeiro). First, the case reproductive number is shifted forward in time relative to the instantaneous reproductive number. It produces leading estimates of changes in the instantaneous reproductive number (Fig 2 and S2 Fig) because it uses data from time points after t, whereas the instantaneous reproductive number uses data from time points before t (Fig 1). Because the case reproductive number is essentially a convolution of the instantaneous reproductive number and the generation interval (Eq 4) [20], shifting the case reproductive number back in time by the mean generation interval usually provides a good approximation of the instantaneous reproductive number [2]. The case reproductive number generally changes more smoothly than the instantaneous reproductive number [26] (Fig 2B and 2C), but if a smoothing window is used, the estimates become more similar in shape and smoothness.

Next is the issue of real-time observation. Estimators of the instantaneous reproductive number were partly developed for near real-time estimation and only use data from before time t (Fig 1A). Under ideal conditions without observation delays and a window size of 1 day, neither method is affected by the termination of the synthetic time series at t = 150 (Fig 2A and 2B). These methods are similarly robust if the time series ends while R_t is rapidly falling (S1A Fig) or rising (S1B Fig). (Below we discuss more realistic conditions, e.g., in which data at the end of a right-truncated time series are incomplete due to observation delays.) Unlike the instantaneous reproductive number, the case reproductive number is inherently forward-looking (Fig 1B and 1C); near the end of a right-truncated time series, it relies on data that have not yet been observed. Extensions of the Wallinga and Teunis method can be used to adjust for these missing data and to obtain accurate R_t estimates to the end of a truncated time series [4,29]. But as shown in Fig 2C, without these adjustments, the method will always underestimate R_t at the end of the time series, even in the absence of reporting delays. Mathematically, this underestimation occurs because calculating the case reproductive number involves a weighted sum across transmission events observed after time t. Time points not yet observed become missing terms in the weighted sum. Similarly, any infections that occurred before the first observed date are missing terms in the denominator of the Cori and colleagues estimator, and so the method of Cori and colleagues often overestimates R_t early in the time series.

Overall, for real-time analyses aiming to quantify the reproductive number at a particular moment in time or to infer the impact of changes in policy, behavior, or other extrinsic factors on transmission, the instantaneous reproductive number will provide more temporally accurate estimates and is most appropriate. The case reproductive number of Wallinga and Teunis considers the reproductive number of specific individuals and therefore is more appropriate for analyses aiming to incorporate individual-level covariates such as age [28] or transmission cluster membership [13] when estimating R_t.

Generation interval misspecification

When estimating R_t from observed data, misspecification of the generation interval is a large potential source of bias. Regardless of the method used, R_t estimates are sensitive not only to the mean generation time but also to the variance and form of the generation interval distribution [20].

The renewal equation is a cornerstone of demographic theory and forms the mathematical backbone of the R_t estimators described above [20]. Within the renewal equation, the generation interval mechanistically links the reproductive number R to observables such as the epidemic growth rate r or the number of new infections per day [20]. Wallinga and Lipsitch [20] describe how the exponential growth rate of the Bettencourt and Ribeiro estimator (their Eq 3.1) and continuous-time equivalents of the Cori and Wallinga and Teunis estimators (their Eqs 4.1 and 4.2) can be derived from the renewal equation model.

Originally developed in the context of population biology, the renewal equation is usually expressed as , where b(t) is the number of births at time t and n(a) is fecundity at age a, scaled by the probability of surviving to age a. When used to describe epidemic dynamics, the renewal equation model is expressed in terms of I(t), the number of infections incident at time t; S(t), the susceptible fraction; and w(·), the generation interval distribution [31–33]. Note that R₀S(t) = R_t.

(7)

The difficulty is that the “intrinsic” generation interval of the renewal equation, which is the interval needed for accurate R_t estimation, is conceptually and quantitatively different from the generation intervals observed in practice [32–34]. In the renewal equation, the generation interval describes the time distribution of all infectious contact events, whether or not the contacted individual is susceptible. (Note that a different factor in the equation, S(t), scales the probability of contacting a susceptible individual and causing a new infection.) The demographic analogy is that w(·) describes changes over time in the fecundity (infectiousness) of an index case, while S(t) determines whether the offspring of that index case are viable. In practice, we only observe generation intervals between viable pairs. Thus, unlike the intrinsic generation interval, observed intervals are sensitive to changes in the rate of susceptible depletion and can be biased estimators of the intrinsic interval and of R_t [32–34]. A related challenge is that interventions such as contact tracing and self-isolation can limit transmission late in the course of infectiousness, shortening observed generation, and serial intervals [33,35,36]. Methods for accurate estimation of the generation interval from contact tracing data involve adjusting for right truncation and accounting for population susceptibility at the times transmission pairs are observed [33,34].

The serial interval, defined as the time between symptom onset in an infector–infectee pair, is more easily observed than the generation interval and often used in its place. Although the serial and generation intervals are often conflated, failure to understand the differences between these related quantities can bias R_t estimates [22,33]. The serial interval and the generation interval have the same mean, but usually have different variances [37,38], and the serial interval can be negative (e.g., for COVID-19 [39–41]), whereas the generation interval cannot [22,33]. The intrinsic generation interval can be estimated from contact tracing data (i.e., estimates of the serial interval) [33,34].

Fig 3A illustrates the consequences of misspecifying the mean generation interval in the method of Cori and colleagues. If the mean generation interval is set too high, R_t values will typically be further from 1 than the true value—too high when R_t>1 and too low when R_t<1. If the mean is set too low, R_t values will typically be closer to 1 than the true value. These biases are relatively small when R_t is near the critical threshold of 1, but these increase as R_t takes substantially higher or lower values (Fig 3). Therefore, biases from misspecification of the generation interval may be greatest early in an epidemic, when the sensitivity of high R_t values to bias may be compounded by limited data and highly uncertain generation interval estimates.

Download:

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

Fig 3.

Biases from misspecification of the generation interval mean (A) or variance (B). Demonstrated using the method of Cori et al. GI, generation interval.

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

The consequences of misspecifying the form and variance of the generation interval distribution are illustrated in Figs 2A and 3B. As explained above, biases in the Bettencourt and Ribeiro estimator arise from misspecification of the form of the generation interval, even if the mean is correctly specified. Similarly, the accuracy of the Cori and colleagues and Wallinga and Teunis estimators in Fig 2B and 2C depends on accurate specification of the generation interval; in practice, if the mean, variance, or form of a pathogen’s true generation interval is uncertain, R_t estimates obtained using these methods can be biased.

EpiEstim [21] allows users to account for uncertainty in the mean and standard deviation of the generation interval by resampling over a range of plausible values [15,21]. Similarly, Bayesian methods such as EpiNow2 [42] and applications of the Bettencourt and Ribeiro method [24] allow users to specify the prior variance of the mean and standard deviation. Uncertainty around an incorrect value can widen the resulting 95% interval but will not shift the assumed central value toward the truth and will not correct bias in the central R_t estimates.

Joint estimation of both R_t and the serial interval is possible, depending on data quality and magnitude of R_t [13,43,44], and the EpiEstim [16,21] package provides an off-the-shelf option for joint estimation. However, these off-the-shelf methods should be used with caution, as they estimate the observed serial interval, not the intrinsic generation interval, and do not account for changes over time in behavior or susceptibility.

Adjusting for delays

Estimating R_t requires data on the daily number of new infections (i.e., transmission events). Due to lags in the development of detectable viral loads, symptom onset, seeking care, and reporting, these numbers are not readily available. All observations reflect transmission events from some time in the past. In other words, if d is the delay from infection to observation, then observations at time t inform R_t−d, not R_t (Fig 4). Obtaining temporally accurate R_t estimates thus requires assumptions about lags from infection to observation. If the distribution of delays can be estimated, then R_t can be estimated in 2 steps: first by inferring the incidence time series from observations and then by inputting the inferred time series into an R_t estimation method. Alternatively, the unobserved time series could be inferred simultaneously with R_t or treated as a latent state. Such methods are now available in a development version of the R package EpiNow2 [42, 45].

Download:

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

Fig 4. R_t is a measure of transmission at time t.

Observations after time t must be adjusted. ICU, intensive care unit.

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

Simple but mathematically incorrect methods for the inference of unobserved times of infection have been applied to COVID-19: convolution and temporal shifts. The errors introduced by these methods may be tolerable if delays to observation are relatively short and not highly variable and if R_t is not rapidly changing. But when dealing with longer or more variable observation delays, or when aiming to infer the timing of changes in R_t accurately, these methods may not be sufficient.

One method infers each individual’s time of infection by subtracting a sample from the delay distribution from each observation time. This is mathematically equivalent to convolving the observation time series with the reversed (backward) delay distribution, but convolution does not accurately infer the underlying time series of infections from observations [46–48]. The forward delay distribution has the effect of spreading out infections incident on a particular day across many days of observation. This blurring into the future is biologically realistic and reflects individual variation in disease progression and care seeking. To recover the original time series of infections from observations requires deblurring. Instead, as illustrated in S3 Fig, backward convolution unrealistically spreads them out further. An unintended consequence of added blurring from backward convolution is that it can help smooth over weekend effects and other observation noise. But a crucial pitfall is that this blurring also smooths over true variation in R_t: peaks, valleys, and changes in slope of the latent time series of infection events. Convolution and other approaches that blur or oversmooth can therefore prevent or delay detection of changes in R_t and can impede accurate inference of the timing of these changes (Fig 5C).

Download:

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

Fig 5. Pitfalls of simple methods to adjust for delays to observation when estimating R_t.

Infections back calculated from (A) observed cases or (B) observed deaths either by shifting the observed curve back in time by the mean observation delay (shift), by subtracting a random sample from the delay distribution from each individual time of observation (convolve), or by deconvolution (deconvolve). Only the deconvolved time series is adjusted for right truncation. Deconvolution most accurately recovers peaks or valleys in the true infection curve. Shifting is less accurate, and convolution is least accurate. Errors from back-calculation increase with the variance of the delay distribution (B vs. A). (C) Posterior mean and credible interval of R_t estimates from the Cori et al. method. Inaccuracies in the inferred incidence curves affect R_t estimates, especially when R_t is changing (here, R_t was estimated using shifted values from panels A and B). Finally, we note that shifting the observed curves back in time without adjustment for right truncation leads to a gap between the last date in the inferred time series of infection and the last date in the observed data, as shown by the dashed lines and horizontal arrows in panels A–C.

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

The second simple-but-incorrect method to adjust for lags is to shift either the raw inputs (the observed time series) or outputs of R_t estimation on the time axis. R_t estimates obtained by applying the methods of Cori and colleagues or Bettencourt and Ribeiro to unadjusted data will lag the true instantaneous R_t by roughly the mean delay from infection to observation. Unadjusted estimates from Wallinga and Teunis are less lagged, because the lead intrinsic to this estimator, relative to the instantaneous reproductive number, partially offsets lags to observation (Fig 2C). For example, the net lag of Wallinga and Teunis applied to unadjusted times of symptom onset is roughly the incubation period (lag) minus the generation time (amount by which case R leads instantaneous R). For pathogens with similar generation and incubation intervals, the net lag would be small, providing a reasonable approximation to the instantaneous R_t.

Unlike backward convolution, temporal shifting does not further blur the observed time series. Thus, if the mean delay is known accurately, this method is preferable to subtracting samples from the delay distribution (Fig 5A and 5B). However, shifting the input time series does not undo the blurring effect of the original delay, which, like backward convolution, can impede accurate inference of changes in R_t. Shifting inputs or R_t estimates by a fixed amount also fails to account for realistic uncertainty in the true mean delay, which will not be known exactly and might change over time.

More reliable methods to reconstruct the incidence time series are now under development. Given a known delay distribution, a potential solution is to infer the unlagged signal using maximum likelihood deconvolution. This method was applied to AIDS cases, which feature long delays from infection to observation [48], and in the reconstruction of incidence from mortality times series for the 2009 H1N1 (swine flu) pandemic [46]. It is now being applied to COVID-19 [8,49]. Fig 5 shows an example of deconvolution applied following the methods in [46]. The method of [48] is implemented in the backprojNP function within the surveillance R package [50,51]. In principle, deconvolution can more accurately estimate the latent time series than temporal shifting or backward convolution, but the method is sensitive to misspecification of the mean, variance, or form of the delay distribution and the stringency of the stopping condition of the deconvolution algorithm. It can also be difficult to quantify uncertainty in the deconvolved time series [41] and to implement deconvolution while adjusting for right truncation.

A potential alternative to deconvolution is R_t estimation models that include forward delays to observation in the inference process or that treat the time series of infections as latent states. Such methods are in development within the R package EpiNow2 [42]. An additional advantage of inferring the time series of infections jointly with R_t is the seamless integration of various sources of uncertainty, e.g., in R_t and reporting. By comparison, the 2-step approach of first transforming the observed time series and then calculating R_t requires users to propagate uncertainty from the back-calculation step into the R_t estimation step. A final advantage of models that include forward delays to observation is that they could facilitate inference from multiple populations or data streams simultaneously [52,53,64]. For example, by assuming that cases, hospitalizations, and deaths all arise from a common infection process, these methods might be able to infer the incident time series of infections more accurately and precisely, potentially while also estimating delays and changes in ascertainment for specific data sources (e.g., outpatient cases).

Deconvolution or R_t estimation methods that include a forward observation process are particularly useful when delays to observation are relatively long and variable and in analyses that require accurate inference of the timing and speed of changes in R_t. If delays to observation are relatively short, or if R_t is not substantially changing, then deconvolution may not be necessary. For example, when working with synthetic case data in which the mean delays to observation are short and known accurately, the underlying infection curve (Fig 5A) and underlying R_t values (Fig 5C) can be recovered with reasonable accuracy simply by shifting the observed time series. But longer and more variable delays to observation worsen inference of the underlying incidence curve (Fig 5B). In turn, this makes it more difficult to infer the speed and timing of abrupt changes in R_t and to relate those changes to policies, behaviors, or other extrinsic epidemic drivers at specific points in time. For example, simply shifting a time series of observed deaths by the mean delay does not accurately recover the underlying curves of infections or R_t (Fig 5B and 5C).

Another advantage of working with observations nearer the time of infection, such as times of symptom onset among newly symptomatic individuals, is that they provide more information about recent transmission events and therefore allow R_t to be estimated in closer to real time (Fig 5C) [47]. Of course, this advantage could be offset by sampling biases and reporting delays. Users will need to balance data quality with the length of the observation delay when selecting inputs.

Further investigation is needed to determine the best methods for inferring infections from observations if the underlying delay distribution is uncertain. If the delay distribution is severely misspecified, all 3 approaches (deconvolution, shifting by the mean delay, or convolution) will incorrectly infer the timing of changes in incidence. In this case, methods such as deconvolution or shifting by the mean delay might more accurately estimate the magnitude of changes in R_t but at the cost of spurious precision in the inferred timing of those changes. Ideally, the delay distribution could be inferred jointly with the underlying times of infection or estimated as the sum of the incubation period distribution and the distribution of delays from symptom onset to observation (e.g., from line list data).

Adjusting for right truncation

Near real-time estimation requires not only inferring times of infection from the observed data but also adjusting for missing observations of recent infections. The absence of recent infections is known as “right truncation.” Without adjustment for right truncation, the number of recent infections will appear artificially low because they have not yet been reported [4,29, 54–58]. Thus, adjusting for right truncation is particularly important in analyses with the goal of near real-time R_t estimation.

Fig 5 illustrates the consequences of failure to adjust for right truncation when inferring times of infection from observations. Subtracting the mean observation delay m from times of observation (“shift” method in Fig 5A and 5B) leaves a gap of m days between the last date in the inferred infection time series and the last date in the observed data. This hampers recent R_t estimation (Fig 5C). Inferring the underlying times of infection by subtracting samples from the delay distribution (“convolve” method in Fig 5A and 5B) dramatically underestimates the number of infections occurring in the last few days of the time series.

The simplest approach is to drop estimates on the last few dates or to flag them as unreliable [8]. But many methods are available to adjust for right truncation, which can improve real-time analysis [42,54,59]. Based on past reporting delays, these methods infer the total number of infections, observed and not-yet-observed, at the end of the time series.

In short, accurate near real-time R_t estimation requires both inferring the infection time series from recent observations and adjusting for right truncation. Errors in either step could amplify errors in the other. Joint inference approaches for near real-time R_t estimation, which simultaneously infer times of infection and adjust for right truncation, are now in development [42].

Accounting for incomplete observation

The effect of incomplete case observation on R_t estimation depends on the observation process. If the fraction of infections observed is constant over time, R_t point estimates will remain accurate and unbiased despite incomplete observation [13,15,44,60,61]. Data obtained from carefully designed surveillance programs might meet these criteria. But even in this best-case scenario, because the estimation methods reviewed here assume all infections are observed, confidence or credible intervals obtained using these methods will not include uncertainty from incomplete observation. Without these statistical adjustments, practitioners and policy makers should beware false precision in reported R_t estimates.

Sampling biases will also bias R_t estimates [60]. COVID-19 test availability, testing criteria, interest in testing, and even the fraction of deaths reported [62] have all changed over time. If changes in ascertainment are well understood and measurable, they can be incorporated into or used to adjust R_t estimates. For example, several R_t estimation dashboards currently adjust for testing effort [6,8,11]. Observed hospitalizations and deaths may be less sensitive to changes in test availability and testing effort, but may be more sensitive to other factors, such as hospital bed availability. Hospitalizations and deaths will also vary in their representativeness of mean transmission rates, depending on which age groups are being infected. No matter what data are used, 1 potential solution is to flag R_t estimates as potentially biased in the few weeks following known changes in data collection or reporting. At a minimum, practitioners and policy makers should understand how the data underlying R_t estimates were generated and whether they were collected under a standardized testing protocol.

Smoothing windows

Because R_t estimators incorrectly assume all infections are observed, day-of-week reporting effects and stochasticity in the number of observations per day can cause spurious variability in R_t estimates, especially if the number of observations per day is low [15]. The Cori method incorporates a sliding window to smooth R_t estimates, but the size of the smoothing window can affect both the temporal and quantitative accuracies of estimates. Larger windows effectively increase the sample size by drawing information from multiple time points but blur what may be biologically meaningful changes in R_t. Some smoothing approaches can also cause R_t estimates to lead or lag the true value (Fig 6).

Download:

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

Fig 6. Accuracy of R_t estimates given smoothing window width and location of t within the smoothing window.

Estimates were obtained using synthetic data drawn from the S→E transition of a stochastic SEIR model (inset) as an input to the method of Cori et al. Colored estimates show the posterior mean and 95% credible interval. Black line shows the exact instantaneous R_t calculated from synthetic data.

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

Lags can be particularly severe when using the conventions suggested by Cori and colleagues, in which R_t is reported on the last date in a given window, rather than on the middle date. Although this convention returns R_t estimates to the last date in the time series, which is convenient for real-time estimation, R_t estimates reported at the end of a window are based entirely on data from the past and therefore lag the instantaneous R_t (Fig 6B). Instead, we recommend reporting R_t at the midpoint of the smoothing window (Fig 6A), which produces estimates that are more accurately oriented in time. An apparent disadvantage is that a centered window precludes R_t estimation on the last w/2 time units, where w is the width of the window. However, we argue that temporal accuracy is preferable and that if daily counts are low, failure to produce estimates on the last few days in the time series realistically acknowledges uncertainty about recent trends. Thus, for Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) and other pathogens with short timescales of infection, near real-time R_t estimation requires large enough daily counts to permit a small window (e.g., a few days).

Although the sliding window increases statistical power to infer R_t, it does not by itself accurately calculate confidence intervals. Thus, underfitting and overfitting are possible. The risk of overfitting in the Cori method is determined by the length of the time window that is chosen. In other words, there is a trade-off in the window length between picking up noise with very short windows and over-smoothing with very long ones. To avoid this, one can choose the window size based on short-term predictive accuracy, for example, using leave-future-out validation to minimize the 1-step-ahead log score [63]. Proper scoring rules such as the ranked probability score can be used in the same way, and a time-varying window size can be chosen adaptively [42].

Conclusions

We tested the accuracy of several methods for R_t estimation in near real time and recommend the methods of Cori and colleagues [15], which are currently implemented in the R package EpiEstim [21]. The Cori method estimates the instantaneous rather than the case reproductive number and is conceptually appropriate for near real-time estimation. The method uses minimal parametric assumptions about the underlying epidemic process and can accurately estimate abrupt changes in the instantaneous reproductive number using ideal, synthetic data.

Most epidemiological data are not ideal, and statistical adjustments are needed to obtain accurate and timely R_t estimates. First, considerable preprocessing is needed to infer the underlying time series of infections (i.e., transmission events) from delayed observations and to adjust for right truncation. Best practices for this inference are still under investigation, especially if the delay distribution is uncertain. The smoothing window must also be chosen carefully, potentially adaptively, and daily counts must be sufficiently high for changes in R_t to be resolved on short timescales. To avoid biases in R_t estimates, the generation interval distribution must be estimated and specified accurately. Finally, to avoid false precision in R_t, uncertainty arising from delays to observation, from adjustment for right truncation, and from imperfect observation must be propagated. The functions provided in the EpiEstim package quantify uncertainty arising from the R_t estimation model but currently not from uncertainty arising from imperfect observation or delays.

Work is ongoing to determine how best to infer infections from observations and to account for all relevant forms of uncertainty when estimating R_t. Some useful extensions of the methods provided in EpiEstim have already been implemented in the R package EpiNow2 [42,45], and further updates to this package are planned as new best practices become established.

But even the most powerful inferential methods, extant and proposed, will fail to estimate R_t accurately if changes in sampling are not known and accounted for. If testing shifts from more to less infected subpopulations or if test availability shifts over time, the resulting changes in case numbers will be ascribed to changes in R_t. Thus, structured surveillance also belongs at the foundation of accurate R_t estimation. This is an urgent problem for near real-time estimation of R_t for COVID-19, as case counts in many regions derive from clinical testing outside any formal surveillance program. Deaths, which are more reliably sampled, are lagged by 2 to 3 weeks and still subject to biases in underreporting. The establishment of sentinel populations (e.g., outpatient visits with recent symptom onset) for R_t estimation could thus help rapidly identify the effectiveness of different interventions and recent trends in transmission.

Code availability

All code for analysis and figure generation is available at https://github.com/cobeylab/Rt_estimation.

Supporting information

Acknowledgments

We are grateful to Michael Höhle for helpful comments. This work was completed in part with resources provided by the University of Chicago’s Research Computing Center.