Implementation

We recall and describe the implementation of methods to estimate the serial interval distribution and reproduction numbers in epidemics. We also propose tools to explore the sensitivity of estimates to required assumptions.

Defining a generation time distribution

The generation time is the time lag between infection in a primary case and a secondary case. The generation time distribution should be obtained from the time lag between all infectee/infector pairs [8]. As it cannot be observed directly, it is often substituted with the serial interval distribution that measures time between symptoms onset. In our software package, the ‘generation.time’ function is used to represent a discretized generation time distribution. Discretization is carried out on the grid [0,0.5), [0.5, 1.5), [1.5, 2.5), etc.… where the unit is a user chosen time interval (hour, day, week…). Several descriptions are supported: “empirical” requiring the full specification of the distribution, or parametric distributions taken among “gamma”, “lognormal” or “weibull”. In the latter case, the mean and standard deviation must be provided in the desired time units.

A function (‘est.GT’) is also provided to estimate the serial interval distribution from a sample of observed time intervals between symptom onsets in primary cases and secondary cases by maximum likelihood.

Estimation of initial reproduction numbers

Reproduction numbers may be estimated at different times during an epidemic. In the following, we recall methods for estimating the “initial” reproduction number, i.e. at the beginning of an outbreak, and for estimating the “time-dependent” reproduction number at any time during an outbreak, as well as the required hypotheses for the methods. Proposed extensions and options implemented in the software are also presented.

Attack rate (AR)

In the classical SIR model of disease transmission, the attack rate (AR : the percentage of the population eventually infected) is linked to the basic reproduction number [9], by $R_0=-\frac{log\left(\frac{1-\mathit{\text{AR}}}{S_0}\right)}{\mathit{AR}-\left(1-S_0\right)}$ where S ₀ is the initial percentage of susceptible population. The required assumptions are homogeneous mixing, closed population, and no intervention during the outbreak.

Exponential growth (EG)

As summarized by Wallinga & Lipsitch [4], the exponential growth rate during the early phase of an outbreak can be linked to the initial reproduction ratio. The exponential growth rate, denoted by r, is defined by the per capita change in number of new cases per unit of time. As incidence data are integer valued, Poisson regression is indicated to estimate this parameter [6, 10], rather than linear regression of the logged incidence. The reproduction number is computed as $R=\frac{1}{M\left(-r\right)}$ where M is the moment generating function of the (discretized) generation time distribution. It is necessary to choose a period in the epidemic curve over which growth is exponential. We propose to use the deviance based R-squared statistic to guide this choice. No assumption is made on mixing in the population.

Maximum likelihood estimation (ML)

This model, proposed by White & Pagano [11], relies on the assumption that the number of secondary cases caused by an index case is Poisson distributed with expected value R. Given observation of (N _0, N ₁, …, N _T ) incident cases over consecutive time units, and a generation time distribution w, R is estimated by maximizing the log-likelihood $\mathit{\text{LL}}\left(R\right)={\displaystyle {\sum }_{t=1}^T\text{log}\left(\frac{e^{-{\mu }_t}{{\mu }_t}^{N_t}}{N_t!}\right)}$ where μ _t = R ∑ _{i = 1} ^t N _{t − i} w _i . Here again, the likelihood must be calculated on a period of exponential growth, and the deviance R-squared measure may be used to select the best period. No assumption is made on mixing in the population.

The approach assumes that the epidemic curve is analysed from the first case on. If this is not the case, the initial reproduction number will be overestimated, as secondary cases will be assigned to too few index cases: we implemented a correction as described in Additional file 1: Supplementary material S1. It is also possible to account for importation of cases during the course of the epidemic.

Sequential bayesian method (SB)

This method, although introduced as “real-time bayesian” by its authors, more exactly allows sequential estimation of the initial reproduction number. It relies on an approximation to the SIR model, whereby incidence at time t + 1, N(t + 1) is approximately Poisson distributed with mean N(t)e ^{(γ(R − 1))}[12], where $\frac{1}{\gamma }$ the average duration of the infectious period. The proposed algorithm, described in a Bayesian framework, starts with a non-informative prior on the distribution of the reproduction number R. The distribution is updated as new data is observed, using $P\left(R|N_0,\dots ,N_{t+1}\right)=\frac{P\left(N_{t+1}|R,N_0,\dots ,N_t\right)P\left(R|N_0,\dots ,N_t\right)}{P\left(N_0,\dots ,N_{t+1}\right)}$. In other words, the prior distribution for R used on each new day is the posterior distribution from the previous day. At each time, the mode of the posterior may be computed along with the highest probability density interval. As before, the method requires that the epidemic is in a period of exponential growth, i.e. it does not account for susceptible depletion; it implicitly uses an exponential distribution for the generation time; and assumes random mixing in the population.

Estimation of time dependent reproduction numbers (TD)

The time-dependant method, proposed by Wallinga & Teunis [13], computes reproduction numbers by averaging over all transmission networks compatible with observations. The probability p _ij that case i with onset at time ti was infected by case j with onset at time tj is calculated as ${\mathit{p}}_{ij}=\frac{N_{i}w\left(t_i-t_j\right)}{{\displaystyle {\sum }_{i\ne k}{N}_{i}w\left(t_i-t_k\right)}}$. The effective reproduction number for case j is therefore R _j = ∑ _i p _ij , and is averaged as $R_t=\frac{1}{N_t}{\displaystyle {\sum }_{\left\{t_j=t\right\}}{R}_j}$ over all cases with the same date of onset. The confidence interval for Rt can be obtained by simulation. Correction for real time estimation, where not yet observed secondary cases are taken into account is possible [14]. It is possible to account for importation cases during the course of the epidemic.

Estimating reproduction numbers

The ‘estimate.R’ function applies the methods described above to a given epidemic curve. Several methods may be used at the same time by listing them in the “methods” argument. In the session code presented in the appendix, a generation time distribution typical of influenza is defined, using a Gamma distribution with mean 2.6 days and standard deviation 1 day [17]. An example dataset (Germany.1918) is loaded from the package.

Initial inspection of the incidence data shows that the exponential growth period takes place during the first 30 days of the epidemic curve. Sensitivity analyses may help refine the choice of an optimal time window for exponential growth (see below). Here, we applied all methods (except the attack rate) on the first 32 days (epidemic peak) of the epidemic curve and reported estimates in Table 1. Surprisingly, although the analysis uses the same data, the estimates range in a relatively large interval, with up to 15% variation (from 1.2 to 1.4). Moreover, confidence or credible intervals do not always overlap (Figure 1A). The fit of each model to the data is however quite similar in all cases, except for the SB method which fits very poorly (Figure 1B).

Method	Default	Optimal
	R	R
	[ 95% CI ]	[ 95% CI ]
EG	1.34	1.56
(optimal time window: 7:22)	[ 1.33 ; 1.36 ]	[ 1.50 ; 1.62 ]
ML	1.21	1.54
(optimal time window: 11:22)	[ 1.16 ; 1.27 ]	[ 1.42 ; 1.66 ]
SB	1.20	1.38
	[ 1.11 ; 1.28 ]	[ 1.25 ; 1.51 ]
TD	1.40	1.40
	[ 1.09 ; 1.73 ]	[ 1.09 ; 1.73 ]

See text for details regarding the methods. All estimates were obtained using the first 32 days of data (default column) or the best fitting time window (“optimal” column). For the SB method, the optimal reported estimate was obtained on day 22, as this date best fits the end of the exponential growth period. For the TD method, daily estimates were averaged over the 32 first days.

Sensitivity analysis

The EG and ML methods require the user to select the time period over which growth is exponential. By default this is taken as the time period from first case to the date of maximum incidence. However, a better choice is possible using the deviance R-squared statistic over a range of possible time periods. The largest R-squared value corresponds to the period over which the model of analysis fitted the data best: we select this period to provide estimates. To look for this time period, the function ‘sensitivity.analysis’ systematically computes the deviance R-squared statistic over a range of time periods chosen by the user. A plot can be obtained that displays the largest R-squared value over time periods of increasing length (see Figure 2A), and the corresponding estimates can be displayed according to the chosen time window (Figure 2B). Here, this analysis shows that the portion of the epidemic curve that best fitted exponential growth in the EG method was of length 15, and more precisely in this case between time units 7 and 22. The estimate of the reproduction number was 1.56 [ 1.50 ; 1.62 ]. For a large choice of time windows, the estimates of the reproduction number remained within the 95%CI of the best fit, suggesting that the estimate was robust to change in the period of exponential growth. The estimates obtained using the best fitting time window for methods are reported in Table 1. Interestingly, when the “best fitting” time period was used for the EG and ML method, the variability between estimates was reduced.

A second issue worth considering is how estimates change according to the choice of the generation time distribution. A function was developed that systematically computes the reproduction ratio over a range of user chosen generation time distribution. In our example, we varied the generation time distribution for the EG method (see Figure 3): as expected, the estimates increased with the mean generation time [18]. Using the same epidemic curve, the reported reproduction ratio ranged between 1.3 and 2 when the mean generation time goes between 1.5 and 5 days.

Comparison of methods

We conducted a simulation study to investigate how estimates of the reproduction number changed with the method of estimation, the role of over-dispersion in the secondary cases distribution and of epidemic curve aggregation at increasingly large time steps. Epidemics were simulated using a branching process, with no restriction on the number of susceptibles to allow exponential growth. For each case, the latent and infectious period were sampled in distributions typical of influenza (gamma distributions with mean+/−sd 1.6+/−0.3 days for latency and 1+/−1 days for infectious period [19]) yielding a generation time interval with mean 2.6 days. The epidemic was started with one incident case at time t=0, then for each incident case, the number of secondary cases was sampled from a negative binomial distribution with mean β I and variance k β I, where I was the individual’s duration of infectious period, β the effective contact rate, and k the overdispersion parameter. The actual time of infection of secondary cases was sampled uniformly during the infectious period of the index case. In this model, the reproduction number is R = β E(I), so that it was possible to calibrate β to obtain the desired value of R. We ran simulations with 3 values for R (1.5, 2 and 3) and overdispersion parameter k to 1 (no overdispersion) and 4 (large overdispersion).

All epidemics were simulated over a period corresponding to approximately the 6 first generations of cases (24 days). For each combination of R and k, 4000 epidemics with more than 20 cases were simulated. Epidemic data were then aggregated daily, by 3 days periods and by 6 days periods. Estimation of the reproduction ratio was made, and the bias and mean squared error were calculated.