key: cord-0487243-pjbebx5g
authors: Johnson, Kory D.; Beiglbock, Mathias; Eder, Manuel; Grass, Annemarie; Hermisson, Joachim; Pammer, Gudmund; Polechov'a, Jitka; Toneian, Daniel; Wolfl, Benjamin
title: Disease Momentum: Estimating the Reproduction Number in the Presence of Superspreading
date: 2020-12-16
journal: nan
DOI: nan
sha: 2c989e9816d1161339efb12b9f39ad153f74cbfd
doc_id: 487243
cord_uid: pjbebx5g

A primary quantity of interest in the study of infectious diseases is the average number of new infections that an infected person produces. This so-called reproduction number has significant implications for the disease progression. There has been increasing literature suggesting that superspreading, the significant variability in number of new infections caused by individuals, plays an important role in the spread of COVID-19. In this paper, we consider the effect that such superspreading has on the estimation of the reproduction number and subsequent estimates of future cases. Accordingly, we employ a simple extension to models currently used in the literature to estimate the reproduction number and present a case-study of the progression of COVID-19 in Austria. Our models demonstrate that the estimation uncertainty of the reproduction number increases with superspreading and that this improves the performance of prediction intervals. In a simplified model, we derive a transparent formula that connects the extent of superspreading to the width of credible intervals for the reproduction number.

A benchmark method for the estimation of the basic reproduction number R was developed by Cori et al. [2] . An influential framework that allows to quantify the phenomenon of superspreading is provided by Lloyd-Smith et al. [8] . In this report we extend the model of Cori et al. [2] to include the phenomenon of superspreading in the sense of Lloyd-Smith et al. [8] . The goal is not to derive more accurate estimates of R but rather to better quantify the uncertainty inherent in this type of estimate.

Ultimately we are interested in the estimation of R and specifically the question whether, given current case numbers, we can claim with statistical guarantees that R ≤ 1 or R > 1. Given the growing body of evidence about the existence and importance of superspreaders, we incorporate this feature into our models. We observe two important effects: first, it becomes increasingly difficult to accurately estimate the population reproduction number R in the presence of superspreading; second, models with superspreading produce prediction intervals for new cases that have improved coverage compared to those without superpreading. Both of these are demonstrated in our Austrian case-study in Section 3.

In particular, the width of a credible interval for R should decrease as a function of total number of cases used during estimation and increase with the extent of superspreading. Let S be the set of days used to estimate R in the nowcasting framework presented in Section 1.1 and assume that the (average) reproduction number does not change over time. One would then expect that credible intervals have width approximately equal to const. k s∈S I s , (1.1) for values of dispersion parameter k much smaller than 1. This is supported in models derived in Section 2.2.1. Small values of k correspond to high superspreading in our framework which follows that of Lloyd-Smith et al. [8] .

The goal of nowcasting is to get accurate estimates of the current state of the pandemic. Given that our observed infections are random observations from an underlying process, our goal is to understand the parameters of that process, particularly with respect to the reproduction number. In addition, we define a time-varying parameter we call the "momentum" of the epidemic, which is a random realization of population infectiousness at a time-point which accounts for superspreading. This is introduced formally in Section 2.1.

A benchmark for estimating R is the method developed in Cori et al. [2] implemented in the software package 'EpiEstim'. An improved extension of this framework is given in Thompson et al. [11] which accounts for variability in the generation interval (defined below). A substantial extension ('EpiNow') of the EpiEstim-package that is used to estimate R on an international level was developed by a group of researchers at the London School of Hygiene and Tropical Medicine [1] . An important overview of the challenges involved in estimating R in the current situation is given in [6] , a comparative analysis of statistical methods to estimate R is given in [10] .

If the epidemic is at an early stage, the reproduction number R and the rate of exponential growth are connected by the Euler-Lotka equation, see for instance [12, 9] for a discussion in the context of epidemiology.

As we follow the framework of Cori et al. [2] , we briefly describe their basic model. Let I 0 be the number of initial infections and I 1 , I 2 , . . . be the number of new infections on days 1, 2, . . .. By (w n ) n≥1 we denote the generation interval distribution. If D m denotes the number of people infected by a specific person on the m-th day after this person got infected, then we have for m ∈ N

.

We assume that a newly infected individual does not cause secondary cases on the same day, corresponding to w 0 = 0. The generation interval can be interpreted as the infectiousness profile of infected persons. The basic model of Cori et al. [2] assumes that the stochastic process (I t ) t∈Z satisfies

for a sequence of numbers R t , t ∈ Z, and where we put I m = 0 for m < 0. The latter convention could be interpreted as assuming that no cases occurred before time 0. In practice it is often assumed that the generation interval distribution is given as a Gamma distribution that has been discretized in such a way that w m = 0 for all m larger than some cut-off number ν. As a result, the sum in (1.2) will only have ν ∈ N summands and to make assertions about I t we only have to consider the case numbers I t−ν , . . . , I t−1 . As ν is a parameter that can vary between diseases, this term is kept and used throughout our model description in Section 2.1. When estimating the time-varying reproduction number, Cori et al. [2] assume that the reproductive number has stayed constant over a window of τ days. Thus (1.2) simplifies to

Note that the reproduction number in the sense of (1.3) does not denote the number of people that actually have been infected by a given individual, rather describes what one expect in an "average" evolution of the epidemic. Furthermore, while R t = R is assumed to be constant over the window of width τ , as this window moves through time the method produces estimates of R t that slowly vary over time.

The motivation for our hierarchical Bayes approach follows the framework of superspreading provided in Lloyd-Smith et al. [8] . Even if the reproduction number R t = R is constant over a small window of time, it might vary between individuals. We consider the reproduction number of a specific person with index x to be drawn randomly as

This distribution has mean R and variance R 2 /k. The degenerate case k = ∞ corresponds to the deterministic case where r x = R for all individuals. Given r x = r, this person causes Pois(r) new infections. If one integrates out the Poisson parameter r, one is left with the unconditional number of descendants which follows a negative binomial distribution with mean R and variance R + R 2 /k as in Section 2.2.1.

Apparently, (1.3) corresponds to the case where each individual has the same basic reproduction number R, i.e., k = ∞. A basic extension of (1.3) that follows the concept of random individual reproduction numbers in the sense of Lloyd-Smith et al. [8] is to assign, on day t, the individual reproduction numbers r t 1 , . . . , r t It to the I t individuals that got infected on this day. This leads to the recursion

where the individual reproduction numbers r m x are i.i.d. according to a Gamma distribution with mean R and dispersion parameter k. Note that for the degenerate case k = ∞, (1.4) recovers (1.3) . This forms the foundation of the model explained in detail in Section 2.1.

The theme of the present paper is close to that of [3] , in which heterogeneity in R between groups is explicitly modeled. While the high-level descriptions of these models sound nearly identical, those models are relevantly different than ours. In particular, [3] are interested in estimating group-specific or timevarying reproduction numbers for different geographical regions and age groups. On one hand, with sufficient group-specific data, this provides tools of a much broader scope than we present here; while on the other hand, it is assumed that within-group variability is negligibly small. Instead, we focus on aggregate data from a single geographical region but do not assume that individual variability is negligible. Rather, this is precisely the variability we are interested in modeling. Similarly, our critiques of the estimability of the reproduction number transfers to their setting as well: if within-group variability exists, group-specific reproduction numbers are more difficult to estimate than previously acknowledged.

As mentioned in the introduction, we identify an unobserved random variable which we term the "momentum" of the pandemic. This follows from a simple notational change in 1.4 according to the observation that a sum of i.i.d. Gamma random variables is also Gamma distributed with the same dispersion parameter. We rewrite 1.4 as

The terms (θ t ) t≥0 are collectively referred to as the "momentum" of the disease and will be treated as a set of nuisance parameters of the offspring distribution as our primary interest lies in estimating the hyperparameter R. Equation (2.1) describes the distribution of I t conditioned on its whole past, i.e., I s , θ s , s < t. Analogously, equation (2.2) describes θ s given its history. The difference in what we understand as the relative past originates from θ t being conceptually determined "after" I t .

For increased clarity of the form of the model and the estimation methods required, we recast our model as a Bayesian Poisson regression using vector notation. This is made painfully explicit by using an arrow as in I for vectors. As our model is estimated over a set of τ days as in Cori et al. [2] , we specify the regression function over this window. To simplify notation, we use [l], for l ∈ N, to be the vector (1, 2, . . . , l). Similarly, [l, m] for l, m ∈ N is shorthand for the vector (l, l + 1, . . . , m), i.e., [l] = [1, l] . This notation will primarily be used for vector indices. Furthermore, the indices of our vectors increase in time. As such, our generation interval truncated to ν days can be condensely written as w [ν] = (w 1 , . . . , w ν ). Similarly, we have I [t−τ +1, t] = (I t−τ +1 , . . . , I t ).

As a regression model for I [t−τ +1, t] = (I t−τ +1 , . . . , I t ), equation (2.1) can be written as

This regression formulation is important as it highlights the latent variables θ that are required to fully determine the generative model. It also focuses attention on which observations are conditioned upon and which are treated as random, i.e., the τ observations to which we fit the model. This is relevant as more than τ nuisance parameters are present, namely ν + τ − 1. Observe that the earliest data point is I t−τ +1 , which itself requires a history of ν momentum values of θ to determine. This also identifies a fixed covariate matrix W which is a function of the generation interval w [ν] .

While we also think of individual reproduction numbers as changing over time due to factors such as changes in social restrictions, the assumption of constant R over a period renders this moot. For now, the dispersion parameter k is assumed constant as it is best estimated with contact tracing data instead of case count data. We set k = .072, in line with the results of Laxminarayan et al. [7] , which estimated the extent of superspreading for COVID-19 from Indian data. This is also within the range of parameter values identified in Endo et al. [4] using Chinese data.

As we have parameterized the gamma prior on θ t to have mean I t R, it is transparent below that R has an inverse-gamma distribution. Hence we use an Inv-Gamma(α, β) hyperprior where these hyperparameters are set such that R has mean 2.6 and standard deviation 2 as in Abbott et al. [1] . This yields α = 3.69 and β = 6.994. We are suppressing notation for conditioning on all observations before time t − τ + 1. Furthermore, given

A full derivation of the posterior distribution of the pair R, θ [t] given

given in Appendix A. We obtain as posterior

where the last line in the display corresponds to the chosen prior for R and the momentum values

In particular, we treat I [t−τ −ν+1, t−τ ] as constant so that the model does not further specify nuisance parameters before time t − τ − ν + 1. This prevents an infinite recursion in historical observations. Hence we need not only a prior for R, but also for

. In order to condense notation for summations in exponents, let S be the index set for the second product; i.e.,

The additional shorthand below drops "s ∈" from s ∈ S. With this notation, the posterior distribution of R given θ and I is

which is Inv-Gamma(k S I s + α, k S θ s + β). A perhaps counter-intuitive observation is that the posterior distribution of R does not depend on the generation interval w [ν] . This is the result of conditioning on θ versus integrating it out as done in Lloyd-Smith et al. [8] . In our case, this approach is infeasible as the dependence on θ is too complex. 

.

The denominator of the variance picks up an additional k term, making credible intervals wider when k is small. The dependence on θ is difficult to remove in this general setting. Section 2.2.1 considers a simpler setting in which θ can be integrated out in order to derive a transparent function for credible interval width.

To estimate this model, we alternate between a Gibbs-step to sample R and a Metropolis-Hastings step to sample θ. As E[θ s |I s , R] = I s R, we can initialize reasonable starting values for θ using various values of R such that we require little burn-in. We find total chain length to be the more important tuning parameter for valid prediction and credible intervals. In all models presented in this paper, we set ν = τ = 13 and w [ν] a discretized gamma distribution with mean 4.46 and standard deviation 2.63. Inference is conducted using the 10 6 samples that remain after a burn-in of 1,000 and thinning by 5. Model validation is presented in Appendix B.

This section describes a so-called "metaday" model, in which the momentum model can be simplified by integrating out the momentum parameters θ. This produces the standard negative binomial model discussed in Lloyd-Smith et al. [8] . Reparameterizing this model and adding suitable priors produces a conceptually and computationally simpler model. In order to apply this model to the current setting, it is necessary to model larger units than daily infections, hence the creation of "metadays". While the results are not as strong, it provides a quick approximation that still performs well. Furthermore, we can derive an explicit connection between the length of credible intervals and k.

In order to directly relate the dispersion parameter k to the width of the credible interval, we consider the trivial generation interval in which an infected person is only infectious for a single day. Furthermore, such a simple formula comes at the cost of making several approximations as explained below. As such, results in this section should be considered heuristic as opposed to concrete.

When the generation interval w is of this form, w [1] = (1) , the model is purely Markovian and the data follow a Galton-Watson process. Recall that a Poisson(λ)-distributed random variable Y , where λ is a hyperparameter distributed according to Gamma(α, β), follows a negative binomial distribution [8] :

(2.4)

Applying (2.4) to the model from Section (2.1) yields the following distribution for the infections I t :

(2.6) Γ(I s + kI s−1 )

The structure of this likelihood suggests estimating R/(R + k) instead of R. When treating I [t−τ ] and k as fixed, Bayes' theorem yields the posterior distribution of R/(R + k):

(2.7) Given this functional form, it is natural to put a beta prior on R/(R + k). For small k, as shown in Figure 1 , this corresponds to putting an appropriate inverse-gamma prior on R, while for larger k this would correspond to a gamma prior. Therefore, to mimic the R ∼ Inf-Gamma(3.69, 6.994) prior distribution used in Section 2.1, we can use a Beta(α = 71.63,β = 3.755) prior on R/(R+k). The posterior distribution of R/(R + k) has a Beta distribution with parameters

While the hyperparameter values are not so small as to be uninformative, they are easily outweighed by the data in most settings. By a change of variables from R/(R + k) back to R, we derive the posterior distribution of R to be

(2.8)

As a final simplifying step, we compute the normal approximation of this posterior [5, Section 4.1] . To this end, the first and second derivatives of the log-posterior density are

Thus, the mode of the posterior iŝ

This yields a normal approximation of the posterior of

Consider the common setting in which β ≈ k · α. Ignoring the prior parameters momentarily, we see that t s=t−τ +1 I s and k t−1 s=t−τ I s already are of this approximate ratio. So long as the estimation window τ is not extremely small, e.g., 1 or 2, the terms in these two sums almost entirely overlap. Furthermore, while the hyperparametersα andβ are of moderate size, they also approximately satisfyβ ≈ k * α. This yields the following simplification of the variance of the normal approximation:

Hence, the approximate length of a credible interval for R behaves like const k t s=t−τ +1 I s .

It is clear that the assumption ν = 1 is highly unrealistic for COVID-19. In order to bridge this gap, we estimate the model over metadays instead of conventional days. Assuming the reproduction number R is fixed and not too large, the negative binomial model estimated using metadays is approximately equivalent to the momentum model estimated using conventional days. We define a metaday to have length equal to the mean of the generation interval, i.e.,

Given the modeling assumptions we have made for COVID-19, this means a metaday comprises approximately 4.87 conventional days. When a model is defined over metadays, setting ν = 1 is equivalent to assuming that someone is infectious over D meta days.

In order to account for non-integer-valued metadays, define D meta = D meta + D f rac , where D f rac ∈ [0, 1). For simplicity, we assume that new infections are uniformly distributed during the day. In order to not confuse subscripts indexing days and metadays, metaday times will be indicated byt instead of t. Lastly, as we are interested in using the most recent data, we care about matching the right endpoint of our time series. As such, we compute the metadays backwards from a reference day t.

Let day t be the maximal day in our data set. We define the corresponding metaday incidence,Ĩt, to bẽ

Infections for previous metadays then sum similarly over the historical data such that the metadays form a partition of days in our data set. It is easiest to represent the entire process if we allow the indices of the summation notation to be real numbers via c2 s=c1

We assume R is constant for τ days in the metaday model as in the momentum model. The corresponding parameter in the metaday model is τ meta := τ /D meta , where we account for non-integer values as before. Estimating parameters and producing forecasts requires similar modifications for non-integer values. Conceptually, however, these correspond to the same sums as before, just over metadays instead of conventional days. This can be represented concisely in the notation for real-valued summation as α =α + To compare the metaday model to the momentum model it is more sensible to compare the cumulative incidence of several, say T, days. We forecast T Dmeta metadaysĨt, wherẽ

Our T cumulative metaday forecast is then given via

In the case study of the next section, we forecast the total weekly cases, i.e., T = 7.

This section primarily focuses on understanding the evolution of the reproduction number in Austria between April and November. As the momentum model effectively needs τ + ν observations to be fit, this is approximately as early as estimates can be provided for Austria. Our goals are three-fold: to demonstrate the increase in estimated variability in R due to superspreading, to provide valid prediction intervals for new cases, and to compare to similar models without superspreading. Some results will be shown for other countries to help establish the validity of our method, but the focus is on Austrian data.

An important component of estimating the reproduction number on a given date is to account for the delay distribution between infections and observed cases as discussed in Gostic et al. [6] . If a delay of length d occurs between infection and observation, then an infection observed at time t actually occurred on day t − d. In this case, we have a "true infection history" that is distinct from the reported case numbers. Abbott et al. [1] estimate and sample true potential infection histories given observed case numbers. As our primary goal is to understand the uncertainty in estimating R as opposed to providing best in class predictions of R for a given date, we ignore this complication. This allows us to take as model input the historical 7-day moving average of reported cases and to compare methods with simple, transparent input. As a result, however, the x-axis of "date" in the following results should be interpreted as "prediction date", not an estimate of the true infections or reproduction number on this date.

Data on the progression of COVID-19 in Austria is shown in Figure 2 . This graph includes curves for the raw infection data as reported by the European Center for Disease Prevention and Control (Raw), the 7-day moving average of Raw (Raw (MA)), each sampled infection history (Sampled Inf.), and the daily median of the sampled infection histories (Sampled Inf. (M)). Observe that the boundary of the "band" created by the sampled infection histories is not smooth, as it is created from 1,000 separate faded lines. Note that using sampled infection histories effectively shifts the time series backward in time.

In order for the infection histories to approximately match the reported case numbers, we have aligned them in time.

As mentioned in Section 2.1, we sample one million total samples of R and the momentum vector θ. To forecast future cases, we use an individual sample of parameters and run the momentum model for a specified period of time. Our graphs show results for the average number of new cases over the following week. There is no additional smoothing of the raw data or predictions. As our input is the 7-day moving average, our prediction is the 7-day-ahead forecast of this moving average. Forecasts and intervals are shown every day between April 1, 2020 and October 31, 2020. In all of the following graphs, we plot predictions and intervals from three models: the momentum model with k = .072, the metaday model of Section 2.2, and the EpiEstim model of Cori et al. [2] . We label the EpiEstim model "Epi*", as the estimates are produced directly via equation (3.1) below instead of using the EpiEstim R package. As in Cori et al. [2] , we fix a generation interval, as opposed to taking samples of a generation interval estimated from a separate data source as in Thompson et al. [11] . As a result, we are not comparing to the best in class model within the EpiEstim/EpiNow framework, but with a model of corresponding complexity to the momentum model. Other improvements to the modeling framework could then be built on top of the momentum model as they have been for the model of Cori et al. [2] .

To estimate the model of Cori et al. [2] , we estimate the parameters of the Cori et al. [2] posterior distribution directly from the infection data:

where a and b are the shape and rate parameter of the gamma prior distribution on R. The the historical 7-day moving average, we estimate this posterior distribution, draw one million samples for R, and run the corresponding data generating process (1.2) for the required number of days. Figure 3 shows the difference between models with and without superspreading on Austrian data. All of the prediction curves track the observed cases similarly as forecasts are only produced for the coming week. The difference in their prediction intervals, however, is significant. Most notably, the intervals for the momentum model with k = .072 are much wider than those of Epi*. Similarly, the metaday variant of this model produces intervals which are wider still. If one looks closely, one can see that the Epi* model lags behind the true data when cases are both increasing and decreasing. On the other hand, the momentum and metaday models "overshoot" the peaks and troughs in the time series. This is likely due to the model estimating greater momentum during these periods.

To assess accuracy, Table 1 shows, for each method, the proportion of true weekly new cases that fall within the prediction intervals over the prediction period. Coverage is shown for the 50% and 90% prediction intervals for the raw infection data. For all countries, the momentum models with superspreading provide coverage which is far closer to the nominal level for both the 50% and 90% intervals. Clearly coverage is still not exact, and all models perform worse on the Czech data. It is still notable that the momentum models provide approximate coverage in these cases even with the inherent messiness of the COVID-19 case data. Figure 4 shows the same graphs but for the Czech Republic and Croatia. The disease progression in the Czech Republic is similar to that of Austria over the shown period. Croatia is a common Austrian tourist destination and the disease progression is markedly different there than in Austria. The estimated coverage probabilities of the prediction intervals are also shown in Table 1 . The story remains the same as before: coverage is far better for the momentum model with superpreading than without. Similarly, Epi* appears shifted relative to the observed cases, particularly for the Czech data. Here we see that the momentum model performs better than the metaday model, particularly around peaks in the time series.

As the reproduction number is unobserved, we are unable to compare our predictions within a supervised setting as we compared our model forecasts. Given the previous discussion though, we see that the additional variability provided by the momentum model is needed to provide prediction intervals with approximate coverage. Figure 5 shows the 90% credible intervals for R derived by the momentum and Epi* models. The figure clearly demonstrates that the intervals for R are drastically different: with superspreading, intervals for R are roughly 2-3 times as wide as those without. This could have potentially large implications for policy making as we know that relatively small changes in the size of R can lead to large differences in the number of new cases if the disease is allowed to progress unchecked.

At the beginning of our estimation period, at the time when restrictions were being relaxed in Austria, it quickly becomes infeasible to claim that the reproduction number is below 1; i.e., the credible intervals estimated around May/June include the value 1. While there have been some fluctuations over the summer, the majority of the period was worsening though remained indistinguishable from R = 1. This fall, however, has seen long periods with reproduction numbers significantly greater than 1, even with our comparatively wide credible intervals. Intervals are, in general, asymmetric, and skewed toward higher values. Figure 6 shows the corresponding graphs but for the Czech Republic and Croatia. The story is again primarily the same, in which the momentum model with superspreading produces much wider credible intervals. One obvious feature of the Croatian data, however, is a steep decline and subsequent steep increase in June. This corresponds to a large increase and plateau in cases as seen in Figure 4 . The Epi* model estimates that R increases to well over 2 within a short period of time before decreasing again to previous levels. Alternatively, in the same period, the momentum model provides a noticeably lower median estimate but with an incredibly wide interval. Further exploration of the feature is warranted, though it is reasonable that such a large deviation over a small window of time should produce significantly more uncertainty in the value of the underlying parameter, particularly when the model is estimated under the assumption that R is constant over τ = 13 days. Within the momentum model, such short-term deviations can be captured by an increase or decrease in disease momentum instead of just an increase in R. On the other hand, this feature appears to show a flaw within the metaday model, as both the estimated R and interval estimate have extreme spikes. This is likely due to the shortterm nature of the case increase and the metaday model only using roughly 3 metadays for estimation. 

In this paper, we provide a simple extension of the Cori et al. [2] model to account for superspreading. This "momentum" model, incorporates unobserved random variables which drive the process of new infections. Even if case numbers and R are relatively small, the presence of superspreaders can increase the momentum of the disease beyond what would be expected if all individuals have the same infectiousness. The momentum model produces wider credible intervals and wider posterior predictive intervals. We find that these wider intervals significantly improve the coverage of the prediction intervals. Lastly, we derived a simple equation to relate the width of credible intervals to the degree of superspreading via a simplification of the momentum model.

Our distributional assumptions are as follows: We want to calculate the joint distribution:

In the last step we used the conditional independence properties for I t and θ t−1 , respectively. Repeating this process to separate I [t−τ +2,t] and θ [t−τ +1,t−1] from the rest yields:

Now, focusing on the last term, we have

In the last equation, we used the fact that the individual θ s are conditionally independent given the vector I. At this point, the terms P (θ s | I [t−τ −ν+1,t−τ ] , R) become problematic. Knowledge of the terms I m for m > s certainly should shed some insight on the value of θ s ; however, it is not clear how this can be feasibly handled. It is not possible to prevent the occurrence of such terms due to the hierarchical nature of this model: the distribution of I s requires previous θ values, which in return demand the inclusion of previous I values ad infinitum. This problem could be avoided by modeling all data from the start of the pandemic, at which point we could confidently set all values of I and θ corresponding to times prior to the onset of the pandemic to 0. This, however, would require treating the value of R as fixed for the entire pandemic, rendering our approach irrelevant as this assumption is clearly false.

As a solution, we propose putting a prior distribution on these problematic θ s such that P (θ s | I [t−τ −ν+1,t−τ ] , R) ∼ Γ(I s k, k/R), essentially disregarding the additional information provided by future observations. Using a different prior, such as setting P (θ s | I [t−τ −ν+1,t−τ ] , R) = δ RIs -which has the appeal of creating terms such as those in Cori et al. [2] -is statistically unsound, as we would draw different θ s from different types of distributions.

All this taken together yields: Using an inverse-gamma prior on R and using the densities of the other terms as discussed before evaluates to the same likelihood as in the main text.

Here we summarize estimation results for simulated data in order to more precisely show the effect of superspreading in a setting in which true parameters are known. The coverage and length of intervals are shown in Figure B .7. All simulations were run using an initial sequence of cases that had constant value of 50. Simulations were repeated 20 times in order to asses coverage probabilities. Of greatest initial import is verifying that the 90% credible intervals for R indeed cover the true value with approximately nominal probability. The case R = 1 is of primary importance, as it represents the bright-line between the epidemic growing or shrinking. That we have nearly exact coverage in this setting is indication that our credible intervals do not achieve coverage merely by being extremely wide. Furthermore, the intervals for θ also cover the true values with the specified probability in the R = 1 and R = 1.5 cases. With our initial sequence of cases and R = .7, the epidemic sometimes dies out, which can be missed by the model. As such, the R = .7 cases has worse coverage over the momentum parameters θ.

After establishing coverage, our motivation for modeling superspreading is verified by looking at the lengths of the credible intervals: for k small, our intervals need to be extremely wide. In fact, the interval for k = .1 is approximately 2.5 times longer than the interval for k = 10 for both R = .7 and R = 1. For R = 1.5, the estimation problem becomes relatively easy as case numbers grow substantially. This leads to very small credible intervals.

As the explicit conditional distribution of the momentum parameters θ is intractable, we present a summary of the samples observed through the MCMC simulation in Figure B .8. This includes all 25 momentum parameters required when τ = ν = 13 as well as R. As R = 1.5 in this setting, one can observe that the scale increase for θ s as s increases. It is clear that the parameters vary widely through MCMC estimation, even though the are initialized at the marginal MLE:θ s = I sR . Multiple chains are run, each with a separate initial value forR. When k is small, variability in θ is large, requiring both tuning of the proposal distribution and long chains to be simulated in order to overcome high auto-correlation in the MCMC draws of θ. 

Estimating the time-varying reproduction number of sars-cov-2 using national and subnational case counts

A new framework and software to estimate time-varying reproduction numbers during epidemics

Modeling the heterogeneity in covid-19's reproductive number and its impact on predictive scenarios

Estimating the overdispersion in covid-19 transmission using outbreak sizes outside china

Bayesian Data Analysis

Epidemiology and transmission dynamics of covid-19 in two indian states

Superspreading and the effect of individual variation on disease emergence

Estimating epidemic exponential growth rate and basic reproduction number

A comparative analysis of statistical methods to estimate the reproduction number in emerging epidemics with implications for the current covid-19 pandemic

Improved inference of time-varying reproduction numbers during infectious disease outbreaks

How generation intervals shape the relationship between growth rates and reproductive numbers

This appendix derives the posterior distribution of R and θ [t−τ −ν+1,t−1] given the relevant observable past, i.e., I [t−τ −ν+1,t] . We briefly restate some basic properties and definitions of our model.Let w i denote the expected proportion of future infections caused by an infected person which occur on day i after infection. Let ν denote the length of infectiousness, i.e., w ν+k = 0 for all k > 0. Lastly, τ denotes the number of days over which we assume R is constant.