Practical considerations for measuring the effective reproductive number, Rt.

Estimation of the effective reproductive number Rt is important for detecting changes in disease transmission over time. During the Coronavirus Disease 2019 (COVID-19) pandemic, policy makers and public health officials are using Rt to assess the effectiveness of interventions and to inform policy. However, estimation of Rt from available data presents several challenges, with critical implications for the interpretation of the course of the pandemic. The purpose of this document is to summarize these challenges, illustrate them with examples from synthetic data, and, where possible, make recommendations. For near real-time estimation of Rt, we recommend the approach of Cori and colleagues, which uses data from before time t and empirical estimates of the distribution of time between infections. Methods that require data from after time t, such as Wallinga and Teunis, are conceptually and methodologically less suited for near real-time estimation, but may be appropriate for retrospective analyses of how individuals infected at different time points contributed to the spread. We advise caution when using methods derived from the approach of Bettencourt and Ribeiro, as the resulting Rt estimates may be biased if the underlying structural assumptions are not met. Two key challenges common to all approaches are accurate specification of the generation interval and reconstruction of the time series of new infections from observations occurring long after the moment of transmission. Naive approaches for dealing with observation delays, such as subtracting delays sampled from a distribution, can introduce bias. We provide suggestions for how to mitigate this and other technical challenges and highlight open problems in Rt estimation.

Introduction

The effective reproductive number, denoted as Re or R, 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 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 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 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 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 = 1, but the biases we observe are smallest in absolute terms in this range.

Another challenge is that depending on the methods used, R 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 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 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 Rt [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 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, 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 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 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. 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 Rt 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 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 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 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 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, R0 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 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.

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, 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 (ti, 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 ti (green arrows in B) or symptomatic on day ts (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.

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: 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).

Fig 2

Accuracy of R 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 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 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 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.

The method of Cori and colleagues estimates R as where I is the number of infection incidents on day t, and w 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) and current infectiousness (w) 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 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, 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 and the exponential growth rate of the epidemic, where g is the mean generation time:

Under the assumption that R evolves through time as a Gaussian process, Eq 3 facilitates efficient Bayesian R 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 estimates, especially when R 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 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],

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) infected case i, relative to the likelihood that any other case in the data infected case i: Then, the individual reproductive number of case j is defined as The case reproductive number at time t is defined as the expected value of R 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 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 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 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 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 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.

The Cori method most accurately estimates the instantaneous reproductive number in real time. It uses only past data and minimal parametric assumptions.

The method of Wallinga and Teunis estimates a slightly different quantity, the case or cohort reproductive number. The case reproductive number is conceptually less appropriate for real-time estimation but may be useful in retrospective analyses, especially those involving individual-level covariates.

In its current form, the method of Bettencourt and Ribeiro [6,14] involves restrictive assumptions that can bias R estimates, but generalized versions of the method could be accurate and computationally efficient.

Generation interval misspecification

When estimating R from observed data, misspecification of the generation interval is a large potential source of bias. Regardless of the method used, R 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 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 R0S(t) = R.

The difficulty is that the “intrinsic” generation interval of the renewal equation, which is the interval needed for accurate R 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 [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 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 values will typically be further from 1 than the true value—too high when R>1 and too low when R<1. If the mean is set too low, R values will typically be closer to 1 than the true value. These biases are relatively small when R is near the critical threshold of 1, but these increase as R 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 values to bias may be compounded by limited data and highly uncertain generation interval estimates.

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.

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 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 estimates.

Joint estimation of both R and the serial interval is possible, depending on data quality and magnitude of R [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.

The intrinsic generation interval is required to correctly define the relationship between R and incident infections.

The intrinsic generation interval is rarely observable, and care must be taken to estimate it from proxies such as the serial interval.

Adjusting for delays

Estimating R 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, not R (Fig 4). Obtaining temporally accurate R estimates thus requires assumptions about lags from infection to observation. If the distribution of delays can be estimated, then R 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 estimation method. Alternatively, the unobserved time series could be inferred simultaneously with R or treated as a latent state. Such methods are now available in a development version of the R package EpiNow2 [42, 45].

Fig 4

R is a measure of transmission at time t.

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

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 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 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: 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 and can impede accurate inference of the timing of these changes (Fig 5C).

Fig 5

Pitfalls of simple methods to adjust for delays to observation when estimating R.

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 estimates from the Cori et al. method. Inaccuracies in the inferred incidence curves affect R estimates, especially when R is changing (here, R 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.

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 estimation on the time axis. R estimates obtained by applying the methods of Cori and colleagues or Bettencourt and Ribeiro to unadjusted data will lag the true instantaneous R 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.

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. Shifting inputs or R 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 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 is the seamless integration of various sources of uncertainty, e.g., in R and reporting. By comparison, the 2-step approach of first transforming the observed time series and then calculating R requires users to propagate uncertainty from the back-calculation step into the R 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 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. If delays to observation are relatively short, or if R 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 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 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 (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 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 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).

Estimating the instantaneous reproductive number requires data on the number of new infections (i.e., transmission events) over time. These inputs must be inferred from observations using assumptions about delays between infection and observation.

Inferring the unlagged time series of infections using deconvolution, or within an R estimation model that includes forward delays, can improve accuracy.

A less accurate but simpler approach is to shift the observed time series by the mean delay to observation. If the delay to observation is not highly variable and if the mean delay is known exactly, the error of this approach may be tolerable. A key disadvantage is that shifting by a fixed amount of time does not account for uncertainty or individual variation in delay times.

Sampling from the delay distribution to impute individual times of infection from times of observation accounts for uncertainty but blurs peaks and valleys in the underlying incidence curve, which, in turn, compromises the ability to rapidly detect changes in R.

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 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 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 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 estimation, which simultaneously infer times of infection and adjust for right truncation, are now in development [42].

Due to reporting delays, infections at the end of a growing time series will be undercounted. To avoid systematic unerestimation of R on the most recent dates, adjust for the right truncation using 1 of many available methods or truncate the time series to the last date with complete reporting.

Accounting for incomplete observation

The effect of incomplete case observation on R estimation depends on the observation process. If the fraction of infections observed is constant over time, R 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 estimates.

Sampling biases will also bias R 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 estimates. For example, several R 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 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 estimates were generated and whether they were collected under a standardized testing protocol.

R point estimates will remain accurate given the imperfect observation of cases if the fraction of cases observed is time-independent and representative of a defined population. But even in this best-case scenario, confidence or credible intervals will not accurately measure uncertainty from imperfect observation.

Changes over time in the type or fraction of infections observed can bias R estimates. Structured surveillance with fixed testing protocols can reduce or eliminate this problem.

Smoothing windows

Because R 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 estimates, especially if the number of observations per day is low [15]. The Cori method incorporates a sliding window to smooth R 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. Some smoothing approaches can also cause R estimates to lead or lag the true value (Fig 6).

Fig 6

Accuracy of R 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 calculated from synthetic data.

Lags can be particularly severe when using the conventions suggested by Cori and colleagues, in which R is reported on the last date in a given window, rather than on the middle date. Although this convention returns R estimates to the last date in the time series, which is convenient for real-time estimation, R estimates reported at the end of a window are based entirely on data from the past and therefore lag the instantaneous R (Fig 6B). Instead, we recommend reporting R 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 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 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, 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].

If R appears to vary abruptly due to underreporting, a wide smoothing window can help resolve R. However, wider windows can also lead to lagged or inaccurate R estimates.

If a wide smoothing window is needed, report R for t corresponding to the middle of the window.

To avoid overfitting, choose a smoothing window based on short-term predictive accuracy [63] or use an adaptive window [42].

Conclusions

We tested the accuracy of several methods for R 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 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 to be resolved on short timescales. To avoid biases in R estimates, the generation interval distribution must be estimated and specified accurately. Finally, to avoid false precision in R, 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 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. 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 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. Thus, structured surveillance also belongs at the foundation of accurate R estimation. This is an urgent problem for near real-time estimation of R 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 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.

Real-time accuracy when R is rising or falling.

(A–C) Alternate version of Fig 2 in which the time series ends on the day R first hits its minimum value after falling abruptly (time 67, yellow point) or 8 days after the changepoint (time 75). (D–F) The time series ends on the day R stops rising (time 97, yellow point) or 8 days later (time 105). Estimates of the instantaneous reproductive number (A, B, D, E) remain accurate to the end of the time series, and estimates do not change as new observations become available in the 8 days following the changepoint. As in the main text, estimates of the unadjusted case reproductive number (C, F) depend on data from not-yet-observed time points. These estimates become more accurate as new observations are added to the end of the time series (orange vs. blue). Methods to infer the number of not-yet-observed infections can help make estimates of the case reproductive number more accurate in real time [4, 29]. All panels show fits to the time series of new infections and assume all infections are observed instantaneously. Solid black line shows the instantaneous reproductive number, and dashed black line shows the case reproductive number. Colored lines and confidence region show posterior mean and 95% credible interval (A, B, D, E) or maximum likelihood estimate and 95% confidence interval (C, F).

(TIF)

Click here for additional data file.

Smoothed estimates of Cori et al. and Wallinga and Teunis.

Both were estimated using a 7-day smoothing window on a synthetic time series of new infections, observed without delay. The estimates of Cori et al. and Wallinga and Teunis are similar in shape when smoothed, but the estimate of Wallinga and Teunis (the case reproductive number) leads that of Cori et al. (the instantaneous reproductive number) by roughly 8 days, or the mean generation interval. Solid colored lines and confidence regions show the posterior mean and 95% credible interval (Cori et al.) or maximum likelihood estimate and 95% confidence interval (Wallinga and Teunis). Dotted and dashed lines show the exact instantaneous reproductive number and case reproductive number, respectively.

(TIF)

Click here for additional data file.

Why is deconvolution needed to recover latent times of infection?

(A) Consider 1,000 individuals all infected at time 100 (vertical line shows the mean). (B) Now consider the times at which these individuals are observed. Logically, tobservation = tinfected+u, where u is a random variable describing the delay between infection and observation. Mathematically, this is a convolution of the infection time and the delay distribution. Because u has non-zero variance, observation times are not only shifted into the future but also blurred across many dates. This blurring is biologically realistic; due to the variability in disease progression and care seeking, individuals with the same date of infection will not necessarily be observed at the same time. (C) Using the observations in panel B, we aim to recover the latent times of infection shown in panel A. Doing so would require not only shifting into the past but also removing the variance introduced by the observation process, which can be achieved by deconvolution. Instead, as demonstrated here, a common strategy is to subtract u from the times of observation, effectively repeating the convolution shown in panel B, but this time moving backward in time rather than forward. This is not the correct inverse operation. It fails to remove variance introduced by the observation process (the forward convolution) and adds new biologically unrealistic variance, further blurring the inferred times of infection. (D) Shifting the times of observation by the mean delay E[u] is also incorrect, as it does not remove the variance from the forward convolution in panel B. But if the mean delay time is known exactly, this approach is preferable to panel C, as it avoids adding even more variance. Ultimately, deconvolution methods would be needed to recover panel A from the observations in panel B while properly accounting for uncertainty.

(TIF)

Click here for additional data file.