The Poisson renewal equation model

While the effective reproductive number can be defined in multiple ways [Reference White26], we focus on the instantaneous effective reproductive number, which is considered particularly suitable for real-time monitoring [Reference Gostic1]. It is defined as the expected number of secondary infections occurring at time $ t $ that are generated by individuals infected before $ t $ , with each infective’s contribution weighted by their relative infectiousness at time $ t $ . The effective reproductive number is often denoted by $ {R}_e $ or $ {R}_t $ , as its variation over time is what makes real-time estimation relevant. However, in EpiEstim and related approaches, this temporal variation is not modelled explicitly. Instead, the reproductive number is assumed to be constant within estimation windows, usually of length $ T=7 $ days. We hence drop the index $ t $ and write simply $ R $ , while indexing the days in the estimation window by $ t=1,\dots, T $ . It is implicitly understood that the estimation exercise will be repeated for a set of rolling time windows, leading to a trajectory of estimates.

Daily or weekly incidence counts $ {X}_t,t=1,\dots, T $ are linked to $ R $ via the renewal equation [Reference Fraser27] and a Poisson assumption,

(1) $$ {X}_t \mid {X}_{t-1},\dots, {X}_0\sim \mathrm{Pois}\left({\mu}_t\right),\hskip0.48em {\mu}_t=R\times {\varLambda}_t,\hskip0.48em {\varLambda}_t=\sum_{d=1}^D{w}_d\times {X}_{t-d}. $$

Here, $ {\mu}_t $ is the expectation of the Poisson distribution, while the probabilities $ {w}_1,\dots, {w}_D $ represent a known discrete-time serial interval distribution. Serial intervals are expressed in the same time steps as incidence counts, and are assumed to have a maximum possible value of $ D $ . We note that the distributional assumption holds conditionally on past values $ {X}_{t-1} $ , $ \dots $ , $ {X}_0 $ , on which $ {X}_t $ depends through $ {\varLambda}_t $ (see Web Appendix 1 of [Reference Cori2]). The Poisson distribution is chosen as it makes Bayesian inference straightforward without the need for Monte Carlo methods [Reference Cori2].

Throughout the paper, we refer to $ {w}_1,\dots, {w}_D $ as the distribution of the serial interval (the time between the onset of symptoms in a pair of consecutive cases) rather than the generation time (the time between two consecutive infections), as the exact times of infection are typically unobserved. It is known that using the serial interval distribution, which has a larger variance than the generation time distribution, introduces some light bias in the estimation of reproductive numbers [Reference Britton and Scalia Tomba28]. However, this has become standard practice [Reference Cori2], and we do not explore this aspect in detail.

Equation (1) can be read as a Poisson regression with an identity link, no intercept and $ {\varLambda}_t $ as the sole covariate [Reference Brockhaus12, Reference Goldstein, Wakefield and Minin13]. The properties of the resulting estimates can hence be studied in the framework of generalized linear regression [Reference Fahrmeir, Kneib and Lang29]. This is most straightforward for maximum likelihood (ML) estimates; the results, however, also translate to the Bayesian approach by [Reference Cori2], which behaves similarly unless incidence is very low (Supplement B.1.2). The ML estimator takes the intuitive form

(2) $$ {\hat{R}}_{\mathrm{Po}}=\frac{\sum_{t=1}^T{X}_t}{\sum_{t=1}^T{\varLambda}_t}, $$

where the index ‘Po’ indicates the underlying distributional assumption. The estimated standard error

(3) $$ \hat{\mathrm{se}}\left({\hat{R}}_{\mathrm{Po}}\right)=\frac{{\hat{R}}_{\mathrm{Po}}}{\sqrt{\sum_{t=1}^T{X}_t}} $$

depends on $ {\hat{R}}_{\mathrm{Po}} $ and the total incidence over the estimation window, but not the variability of the time series conditioned on past values. This is a consequence of the equidispersion property of the Poisson distribution, which implies $ {\sigma}_t^2={\mu}_t $ , with $ {\sigma}_t^2 $ the conditional variance $ \mathrm{V}\mathrm{a}\mathrm{r}({X}_t\hskip2pt |{X}_{t-1},\dots, {X}_0) $ . Derivations of equations (2) and (3) are provided in Supplement B.1.1, where we also discuss the challenges due to the dynamic nature of model (1).

The estimated standard error (3) can be used to derive Wald-type confidence intervals, which for a confidence level of 95% are given by $ {\hat{R}}_{\mathrm{Po}}\pm 1.96\times \hat{\mathrm{se}}\left({\hat{R}}_{\mathrm{Po}}\right) $ . Empirical coverage fractions of such intervals in simulation studies are a straightforward way to assess uncertainty quantification under varying conditions. We will address this in Section ‘Simulation study’.

Overdispersion and how to account for it

Impact on the Poisson estimator

In practice, infectious disease incidence data usually display overdispersion relative to the Poisson model. Throughout the paper, we conceive this as conditional variances exceeding conditional means, $ {\sigma}_t^2>{\mu}_t $ . Overdispersion may be due to a variety of factors, including superspreading events [Reference Lloyd-Smith20] and incomplete reporting [Reference Bracher and Held30]. The Poisson assumption will then cause our inference for $ R $ to be overconfident; see Supplementary Section A for an intuitive illustration. Well-known theory on generalized linear regression [Reference Fahrmeir, Kneib and Lang29, Sec. 5.5] implies that expression (3) underestimates the true standard error of $ {\hat{R}}_{\mathrm{Po}} $ by a factor of $ \sqrt{\phi } $ , with

(4) $$ \phi =\frac{\sum_{t=1}^T{\sigma}_t^2}{\sum_{t=1}^T{\mu}_t}. $$

The variance inflation factor $ \phi $ depends on the relationship between $ {\mu}_t $ and $ {\sigma}_t^2 $ . We note that equation (4) holds irrespective of the exact form of this relationship, even though in practice one usually makes more specific assumptions (see the next section). As we shall see in Section ‘Illustration for three outbreaks’, $ \phi $ may be quite large in practice, meaning that standard errors estimated under the Poisson assumption may be considerably too small.

Translating equation (4) to the width of the resulting confidence interval, a $ \left(1-\alpha \right) $ confidence interval from the Poisson renewal equation model has coverage probability

(5) $$ 1-2\times \Phi \left[{\Phi}^{-1}\left(\frac{\alpha }{2}\right)/\sqrt{\phi}\right], $$

where $ \Phi $ is the cumulative distribution function of the standard normal distribution. This matches the intended level $ \left(1-\alpha \right) $ only if $ \phi =1 $ , and decreases if $ \phi >1 $ . The derivation of the coverage probability (5) is provided in Supplement B.2.2.

In the following, we will outline three options to account for overdispersion, all of which have been applied in the literature. While this is typically motivated by the desire to estimate dispersion parameters [Reference Ho16, Reference Zhang, Britton and Zhou18] and general model realism [Reference Gressani31], our focus is on uncertainty quantification.

Negative binomial regression

The most common full probability model to generalize the Poisson assumption is the negative binomial distribution. We here parameterize it by its mean and a dispersion parameter, which, for example, in the R software is called the ‘size parameter’. Despite its name, a higher value of the dispersion parameter indicates lower dispersion (see below for formulas), and the Poisson distribution is the limiting case of an infinite dispersion parameter. The dispersion parameter can be handled in two ways ([Reference Cameron and Trivedi32], p.73), which we call NegBin-L and NegBin-Q, as they imply linear and quadratic mean–variance relationships. Irrespective of which is used, estimates and standard errors for $ R $ are obtained using standard regression methodology and naturally account for overdispersion.

NegBin-L

In the NegBin-L renewal equation (see, e.g., [Reference Ho16, Reference Zhang, Britton and Zhou18]),

(6) $$ {X}_t\mid {X}_{t-1},\dots, {X}_0\sim \mathrm{N}\mathrm{e}\mathrm{g}\mathrm{B}\mathrm{i}\mathrm{n}(\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}={\mu}_t,\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p}=\xi \times {\mu}_t), $$

the dispersion parameter is assumed to be proportional to the mean. This implies a linear mean–variance relationship $ {\sigma}_t^2=\left(1+1/\xi \right)\times {\mu}_t $ as in the quasi-Poisson model, and can be motivated by overdispersed individual offspring distributions [Reference Zhang, Britton and Zhou18]. While some differences between NegBin-L and quasi-Poisson regression exist [Reference Hilbe33], in our applications we found similar behaviour, with point estimates almost identical to those from equation (2).

NegBin-Q

We therefore focus primarily on NegBin-Q, defined as

(7) $$ {X}_t\mid {X}_{t-1},\dots, {X}_0\sim \mathrm{N}\mathrm{e}\mathrm{g}\mathrm{B}\mathrm{i}\mathrm{n}(\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}={\mu}_t,\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{p}=\psi ), $$

and implying a quadratic mean–variance relationship $ {\sigma}_t^2={\mu}_t+{\mu}_t^2/\psi $ . This version, employed, for example, by [Reference Brockhaus12, Reference Zhang, Britton and Zhou18, Reference Gressani31, Reference Held, Höhle and Hofmann34], can be motivated by temporal variation of $ R $ over the estimation window [Reference Zhang, Britton and Zhou18]. Despite the mean structure $ {\mu}_t=R\times {\varLambda}_t $ mirroring the Poisson version, the ML estimator $ {\hat{R}}_{\mathrm{NBQ}} $ is not identical to $ {\hat{R}}_{\mathrm{Po}} $ from equation (2). No closed expression exists, but if overdispersion is sufficiently strong, the approximation

(8) $$ {\hat{R}}_{\mathrm{NBQ}}\approx \frac{1}{T}\times \sum_{t=1}^T\frac{X_t}{\varLambda_t} $$

holds; see derivation in Supplement B.3.1. Unlike in equation (2), the ratios $ {X}_1/{\varLambda}_1,\dots, {X}_T/{\varLambda}_T $ thus contribute equally to the estimate. This is because they are all equally informative under a quadratic mean–variance relationship, irrespective of the magnitude of $ {\varLambda}_t $ and $ {X}_t $ (see discussion in [Reference Ver Hoef and Boveng35]). This also explains why the approximate standard error

(9) $$ \hat{\mathrm{se}}\left({\hat{R}}_{\mathrm{NBQ}}\right)\approx {\hat{R}}_{\mathrm{NBQ}}\times \sqrt{{\left(\hat{\psi}\times T\right)}^{-1}} $$

features the number of observations $ T $ rather than the total incidence as in equation (3). It moreover depends on the estimated dispersion parameter $ \hat{\psi} $ , which reflects the conditional variability of the time series and can be computed from the model residuals (Supplement B.3). It generally holds that $ \hat{\mathrm{se}}({\hat{R}}_{\mathrm{NBQ}})\ge \hat{\mathrm{se}}({\hat{R}}_{\mathrm{Po}}) $ , that is, the estimated standard error is at least equal to its Poisson counterpart. The approximations (8) and (9) are assessed empirically in Supplementary Figure S22.

ML estimation in the negative binomial renewal equation models (6) and (7) can run into convergence issues if the data used for fitting show no or very little overdispersion. This is because the dispersion parameter $ \xi $ in the NegBin-L model (or $ \psi $ in the NegBin-Q) is usually estimated on a log-transformed scale to ensure positivity. If the highest likelihood is achieved for a very high value of $ \xi $ , the likelihood as a function of $ \mathrm{log}(\xi ) $ becomes almost flat for large $ \mathrm{log}(\xi ) $ . This causes numerical failure in the computation of point estimates and/or standard errors. Whenever this occurs, we revert to the Poisson model (1), and report the corresponding point estimates and standard errors. In such low-dispersion settings, reverting to the Poisson is numerically convenient and statistically appropriate, as the Poisson distribution can be understood as the limiting form of the negative binomial when overdispersion vanishes.

Quasi-Poisson regression has certain advantages over the negative binomial versions if the mean–variance relationship is mis-specified [Reference Wooldridge36]. In our setting, however, both equations (2) and (8) remain unbiased and consistent in this case, and only their efficiency is reduced. How strongly the two estimators differ is an empirical question we will return to in Section ‘Illustration for three outbreaks’.

Comparing model fits and testing for the presence of conditional overdispersion

To assess whether and to what degree the negative binomial models improve the fit relative to the Poisson, we compare the AIC values and employ likelihood ratio (LR) tests. The quasi-Poisson model cannot be included in this comparison, as it does not have a proper likelihood function. In the LR test, we test the null hypothesis of no conditional overdispersion. Specifically, we test $ {H}_0:1/\xi =0 $ (NegBin-L) and $ {H}_0:1/\psi =0 $ (NegBin-Q) against one-sided alternatives. To account for the fact that the null hypotheses are on the border of the parameter space we employ the correction suggested by Molenberghs and Verbeke [Reference Molenberghs and Verbeke37].

Setup of the simulation study

To support our claim that Poisson assumptions lead to overconfidence in the estimation of $ R $ we conduct a simulation study. We combine different distributional assumptions for the data-generating process and the model used for $ R $ estimation, and compute empirical coverage fractions of Wald confidence intervals. These are given by the fraction of cases in which confidence intervals with a given nominal level contain the true reproductive number $ R $ .

In the main manuscript, we focus on time series generated from the NegBin-L renewal equation (6). We generate 1,000 trajectories from each of four scenarios, which are defined by a combination of two parameters: two values of the reproductive number $ R\in \{1.5,\hskip2pt 2.5\} $ and two levels of overdispersion $ \xi \in \{2,\hskip2pt 0.2\} $ . For all scenarios, we used the same shifted gamma serial interval distribution with a mean of 7.5 days and a standard deviation of 2.1 days, cut off at $ D=14 $ days. All trajectories are initialized by the same sequence of values $ \Big(2 $ , $ 7 $ , $ 6 $ , $ 4 $ , $ 3 $ , $ 7 $ , $ 7 $ , $ 6 $ , $ 6 $ , $ 6 $ , $ 4 $ , $ 9 $ , $ 3 $ , $ 7\Big) $ . In Supplementary Section C.2.1, we include additional simulation scenarios, where the trajectories are initialized by higher values around 100 in order to produce different trajectories with overall higher incidence. After this deterministic initialization, we generate random trajectories of length 28 days. The first 14 days serve to avoid dependence on the specific initialization sequence (‘burn-in period’), while the last 14 days are used to estimate the value of $ R $ . We emphasize that the true $ R $ is a constant throughout the 28 days, such that in terms of the mean structure the assumptions shared by all considered models are fulfilled and the theoretical arguments from the previous sections apply. While in reality, smooth variation of $ R $ would be expected, our simulation study is concerned with the behaviour of EpiEstim-like methods if all assumptions other than the Poisson distribution are fulfilled. It therefore inherits the constraint of time-constant $ R $ values. An adapted simulation where the reproductive number varies smoothly over time is available in Supplementary Section C.2.3.

The serial interval is motivated by the estimate [Reference Crowcroft38] provide for respiratory syncytial virus (RSV; see also [Reference Vink, Bootsma and Wallinga39]). In Supplementary Section C.2.2, we provide further simulation scenarios with serial intervals that are motivated by influenza (mean 3.7 days, sd 2.1 days [Reference Morgan40]) and measles (mean 13.7 days, sd 1.5 days [Reference Fine41]).

For each simulated trajectory, $ R $ values were estimated under the four distributional assumptions discussed in Section ‘Overdispersion and how to account for it’. The Poisson and quasi-Poisson models were fitted using the stats::glm function in R [42], while gamlss::gamlss [Reference Stasinopoulos and Rigby43] was employed for the negative binomial models. To assess the impact of small sample biases, we applied two different estimation window lengths (7 and 14 days) across all scenarios. As described in Section ‘Negative binomial regression’, in cases of convergence issues due to low degrees of observed overdispersion, we reverted from the negative binomial models to the Poisson version.

To make our assessment comprehensive, we ran corresponding simulation studies using the NegBin-Q and Poisson versions of the model as the data-generating processes. In addition, and similarly to [Reference Challen and Danon44], we study a setting where none of the considered models matches the data-generating process. To this end, we simulated time series from an overdispersed branching process model. Details on all of these settings are provided in Supplementary Section C.