Quantifying superspreading for COVID-19 using Poisson mixture distributions

Introduction

For any communicable disease, the basic reproduction number, R₀, denotes the average number of secondary cases a single infected individual generates in a completely susceptible population [1, 2]. The basic reproduction number is considered to be of constant value among population members or specific population groups. However, for person-to-person transmitted infections, a complex combination of host, pathogen, and environmental factors defines the transmission potential of an infected individual, i.e. the number of other individuals a case infects during their infectious period [3, 4]. It has been shown that, for a given R₀, both the probability that an epidemic will occur and the subsequent course of the epidemic are affected by individual variation in transmission [3]. Variation in disease transmission may raise the existence of ‘superspreaders’ who infect substantially more individuals than others. When superspreading plays an important role during the epidemic, a relatively small part of infected cases will be responsible for most of the transmission, while many cases do not transmit the disease at all. Furthermore, when variation in disease transmission is present, large outbreaks can occur even if R₀ is less than one. To account for this heterogeneity, the individual number of secondary cases can be described by a random variable, whereas R₀ represents the expected value for an entire susceptible population.

The transmission potential of infected individuals can be defined as a combination of their biological infectiousness (i.e. viral shedding) and their contact behaviour [5]. It is reasonable to assume that individuals with a higher viral shedding will be more likely to transmit the infection given a contact. In addition, for a fixed level of viral shedding, infectious individuals with a higher contact rate will be more likely to generate secondary cases. Regarding SARS-CoV-2, significant individual variation in viral shedding has been reported [6] and it has been argued that small aerosols exhaled during normal speech may serve as an important transmission route [7]. Vuorinen et al. [6] investigated the possibility of SARS-CoV-2 transmission by inhalation of virus-containing aerosols, by examining a high-risk scenario where an infected individual coughs within a public indoor space. They found that there was an elevated risk of infection in case of lengthy exposure in a confined space with at least one infected individual. These results are in line with those from another study which has indicated that the virus may remain infectious as an aerosol for at least three hours [8]. Of course, not only individual characteristics such as viral shedding but also environmental characteristics such as insufficient ventilation contribute to the possibility of a superspreading event (SSE) [9].

Lloyd-Smith et al. [3] addressed heterogeneity in transmission by using the concept of an individual reproduction number as a random variable that represents the expected number of secondary cases caused by a particular infected individual. In that framework, SSEs are important realizations from the right-hand tail of the distribution of the individual reproduction number. Most studies investigating the amount of heterogeneity in disease transmission have assumed a Poisson process with rate given by the individual reproduction numbers, assumed to follow a Gamma distribution, resulting in a negative binomial offspring distribution [3, 10]. In this way, heterogeneity has often been quantified using the k parameter, with k the negative binomial dispersion parameter. This has allowed comparison between studies, where lower values of k indicate increased heterogeneity in transmission, and thus possibly a larger amount of superspreading.

Based on this framework, a substantial amount of individual variation in the transmission of SARS-CoV-2 has been described, though large differences were found between different studies. Bi et al. [11] used a negative binomial distribution to describe superspreading in the COVID-19 outbreak in Shenzhen, China, and found that about 9% of all cases were responsible for 80% of transmission. Riou and Althaus [12] estimated the negative binomial dispersion parameter k to have a median of 0.54 (90% HDI 0.014 - 6.95), with simulations suggesting that very low values of overdispersion (<0.1) are less likely. Adam et al. [13] estimated the overall mean number of secondary cases to be 0.58 (95%CI 0.45 - 0.72) with a dispersion parameter k of 0.43 (95%CI 0.29 - 0.67) in Hong Kong, indicating that 19% of cases were responsible for 80% of all local transmission. Similarly, Endo et al. [14] have used a branching process model where the number of secondary cases was assumed to follow a negative binomial distribution. Assuming R₀ to be 2.5, they estimated the dispersion parameter k to have a median of 0.1 (95% CrI 0.05 - 0.2), resulting in 80% of secondary cases being caused by about 10% of infectious cases and implying that large transmission events should be prevented in order to contain epidemic spread. Laxminarayan et al. [15] estimated the negative binomial dispersion parameter k to be 0.51 (95%CI: 0.49-0.52) using a large contact tracing dataset from two Indian states. Based on detailed contact tracing data from Hunan, China, Sun et al. [16] found that 15% of cases were responsible for 80% of transmission, and a negative binomial dispersion parameter k of 0.3. Lau et al. [17] found that superspreading was widespread across space and time, with an increasing presence towards later stages of the investigated outbreaks, highlighting the importance of maintaining social distance measures. They also found that about 2% of the most infectious cases were directly responsible for 20% of all infections.

It is well recognized in statistical literature that the distribution underlying a data generating mechanism imposes a certain mean-to-variance relationship which in practice may be severely violated [18]. Despite this, the use of other distributions in infectious disease modelling that may just as well account for variation in disease transmission has been rather limited. Some studies suggest that SSEs follow a power-law distribution with fat tails, such as the generalized Pareto distribution [19]. Brooks-Pollock et al. [20] have used a negative binomial as well as a Poisson-lognormal distribution to model the distribution of cluster sizes for tuberculosis in the UK and the Netherlands. In this study, the Poisson-lognormal distribution provided a better fit to the UK data, indicating the importance of comparing different assumptions about the underlying distribution when variation in disease transmission is present.

To our knowledge, there are no studies that have explicitly investigated the possible bias in using the negative binomial distribution as an approximation to the underlying transmission process. We argue that it is important to compare different distributional assumptions since different distributions could portray different tail behaviour, and hence capture SSEs differently. In this work we explore the use of other Poisson mixture distributions for inference of the offspring mean and the amount of heterogeneity in disease transmission. We focus on the three-parameter generalized Gamma distribution for the individual reproduction number, because of its flexibility and the fact that it has as special cases the Gamma, Weibull, and lognormal distribution [21]. First, we carry out a simulation study to investigate the potential bias in the estimation of the offspring mean and its variance when the distribution that is fit to the data does not correspond to the actual data generating distribution. Next, we use the proposed distributions to (re-)analyze several COVID-19 datasets and investigate the impact of the considered offspring distribution on the estimation of the proportion of cases that is responsible for 80% of transmission, p_80%.

Results

Expected versus realized proportions of transmission

Based on the estimated mean and variance of the considered offspring distribution, we can obtain estimates for the proportion of cases responsible for a certain amount of transmission. There are two different approaches for obtaining these proportions (see Methods), where one is based on the distribution of the individual reproduction number [3], and the other is based on the complete offspring distribution [14]. Here we show how the different offspring distributions can affect these proportions. Figure 1a shows the expected proportion of transmission due to the 20% most infectious cases for the varying levels of heterogeneity used in the simulation study, for the different offspring distributions. For all distributions, this proportion increases with an increasing amount of heterogeneity (i.e. higher σ, which for the negative binomial results in lower k). Thus, less heterogeneity leads to a smaller proportion of transmission being attributed to the 20% most infectious cases. In case of high overdispersion, there is a substantial difference in the expected proportions between the distributions. Since for the Poisson-generalized Gamma distribution it is not possible to specify the parameters from a given mean and variance, we only estimated these proportions at the specified settings used in the simulation study, hence these are represented as dots instead of lines. Figure 1b shows the expected proportion of transmission due to a given proportion of infectious cases for the different levels of overdispersion. It can again be seen that the difference between these estimates across the different distributions increases when overdispersion increases (i.e. higher σ) and likewise there is a substantial difference in terms of the expected proportions of transmission. When taking into account the additional variation coming from the Poisson process, the same increase in the proportion of transmission due to the 20% most infectious cases is seen, again with substantial differences between the distributions when overdispersion increases (Fig. 1c and 1d). These results also hold when R > 1 (Suppl. Fig. A.2). The vertical bars in Fig. 1c represent the uncertainty in the proportion of transmission due to the discrete nature of the offspring distribution. This should be interpreted as the range of transmission that will be due to the 20% most infectious cases. When taking into account this uncertainty surrounding these point estimates, a substantial difference is seen between the Poisson-lognormal and the other distributions for higher levels of overdispersion.

• Download figure
• Open in new tab

Figure 1:

The top panel shows the expected proportion of all transmission that is (a) due to the 20% most infectious cases for different levels of overdispersion and different distributions, with the offspring mean R fixed at 0.8; and (b) due to a given proportion of infectious cases, where cases are ranked by their transmission potential, for σ = 1 (dotted), σ = 1.5 (dashed), σ = 3 (full), and the different distributions, with R fixed at 0.8. The lower panel shows the realized proportion of all transmission that is (c) due to the 20% most infectious cases, shaded vertical bars show the uncertainty surrounding the proportions at σ = {1, 1.5, 3}; and (d) due to a given proportion of infectious cases.

Application to COVID-19 data

Table 1 shows estimates of the offspring mean R and standard deviation σ obtained by fitting the different offspring distributions to each COVID-19 dataset. In terms of AIC, in all cases the data are best described by a Poisson-lognormal distribution. Supplementary Figs. C.1, C.3, and C.5 show the fit of the different distributions to the observed offspring distribution. It can be seen that for the data from Hong Kong and India the negative binomial, Poisson-Weibull, and Poisson-generalized Gamma distributions underestimate the proportion of cases that generate only one secondary case, while this is captured well by the Poisson-lognormal distribution. For the data from Rwanda this is less evident. Goodness-of-fit plots are shown in Supplementary Figs. C.2, C.4, and C.6. For each of the three datasets (Table 1), p_80% is estimated to be substantially higher for the Poisson-lognormal than for the other distributions when based on the distribution of the individual reproduction number (Eq. (1), see Methods), and slightly higher when taking into account additional random variation from the Poisson process (Eq. (2), see Methods). Figure 2 shows these expected (left) and realized (right) proportions of transmission due to a certain proportion of cases for each distribution and each dataset. For Hong Kong, roughly 12-31% of cases are responsible for 68-100% of all transmission based on the Poisson-lognormal distribution (Fig. 2b). Based on the negative binomial distribution, roughly 14-31% of cases are responsible for 71-100% of all transmission.

View this table:

• View inline
• View popup
• Download powerpoint

Table 1:

Estimates of the offspring mean R and its standard deviation (σ) using the different mixture distributions, and their AIC value, for three COVID-19 datasets. p_80% represents the proportion of cases responsible for 80% of transmission, following Equations (1) and (2). Estimates based on the negative binomial distribution correspond to those reported in the literature for the two published datasets.

• Download figure
• Open in new tab

Figure 2:

Proportion of most infectious cases responsible for a certain proportion of transmission, based on estimates from (a-b) Hong Kong, (c-d) India, and (e-f) Rwanda. Proportions are obtained based on the distribution of the individual reproduction number (left), and based on the complete offspring distribution (right). The shaded areas in the right panels represent the uncertainty surrounding specific proportions when considering the discrete nature of the realized offspring distributions.

For India, roughly 10-29% of cases are responsible for 61-100% of all transmission based on the Poisson-lognormal distribution (Fig. 2d). Based on the negative binomial distribution, roughly 11-29% of cases are responsible for 64-100% of all transmission. For Rwanda, roughly 5-19% of cases are responsible for 44-100% of all transmission based on the Poisson-lognormal distribution (Fig. 2f). Based on the negative binomial distribution, roughly 5-19% of cases are responsible for 46-100% of transmission.

Discussion

Since most studies that aim to quantify variation in disease transmission have assumed the offspring distribution to follow a negative binomial, we investigated the impact of incorrectly assuming this distribution as an approximation to the underlying transmission process. Results from our simulation study show that when overdispersion increases, estimates of the offspring mean and especially its variance can become extremely biased when making incorrect assumptions about the underlying data generating distribution. When no variation in transmission is present, all distributions performed equally well, although there was a slightly increased bias in variance estimates when using the Poisson-Weibull or Poisson-generalized Gamma distribution. We have (re-)analyzed three COVID-19 datasets and in each case the Poisson-lognormal distribution gave the best fit to the observed data. This resulted in considerable differences in terms of the expected p_80% compared to when using a negative binomial distribution, when these proportions were based on the distribution of the individual reproduction number [3]. When accounting for the additional variation introduced by the Poisson process, the differences in these proportions of cases responsible became more subtle. For example, for the Hong Kong data the point estimate when using the Poisson-lognormal was 19.5%, compared to the previously reported 19.1% when using the negative binomial [13]. When accounting for the discrete nature of the offspring distribution, estimated ranges for these proportions were mostly overlapping for the different distributions, albeit a bit lower for the Poisson-lognormal distribution. Although most studies report p_80%, the right-sided panels in Figure 2 indicate that depending on the proportion of transmission one is interested in, there might be a more substantial difference between the distributions. This implies that different distributions have different tail properties, underlining the importance of investigating which distribution best describes the data at hand. In addition we found that the difference between the distributions increases with an increasing level of overdispersion (Figure 1).

Our analyses indicated that the negative binomial distribution often underestimates the proportion of cases that generate only one secondary case, thereby possibly overestimating the importance of superspreading events. This overestimation was observed when the proportion of cases responsible was obtained based on the distribution for the underlying individual reproduction number. When accounting for the Poisson process, super-spreading was found to be only slightly more important when using a negative binomial distribution to describe the data, compared to a Poisson-lognormal distribution. A negative binomial distribution enables easy comparison between different studies through its dispersion parameter k [3]. However, this should not be a reason to only use negative binomial offspring distributions. The results from different studies can also be compared by their estimated p_80%, which is often reported as well. It should be noted that there are different approaches for obtaining these proportions, hence care should be taken when comparing these results between studies. Lloyd-Smith et al. [3] assume SSEs to be realizations from the right-hand tail of the distribution of the individual reproduction number, hence their approach is based on this continuous distribution. In contrast, Endo et al. [14] have based these estimates on the complete offspring distribution, taking into account additional variation arising from the discrete Poisson process. In this way, the second approach accounts for more heterogeneity. However, if the effective contact process in reality is not a Poisson process, this approach may result in biased estimates.

In this study we have considered the three-parameter Poisson-generalized Gamma distribution, which has as special cases the Poisson-lognormal, Poisson-Weibull, and negative binomial distribution. Although the Poisson-generalized Gamma distribution has the advantage of being very flexible due to the additional parameter, the disadvantage is that because of this added complexity the estimation is computationally more extensive, especially for large datasets. Furthermore, parameter estimation can be difficult because different parameter sets can give rise to the same density function. In general, parameter estimation becomes more difficult when the amount of overdispersion is high and incorrect assumptions about the underlying data generating distribution are made. For that reason we were not able to fit the Poisson-Weibull and Poisson-lognormal distributions in some scenarios of our simulation study. This occurred when the data were highly overdispersed (k < 0.1), which is less likely to be encountered in practice [12].

Inference of the amount of heterogeneity in transmission is paramount for identifying a disease’s potential of superspreading. Correctly quantifying this heterogeneity is important because it affects estimates of other epidemiological parameters, modulates the degree of unpredictability of an epidemic, and needs to be taken into account when modeling disease control and planning control strategies [22]. When there is evidence of substantial superspreading, control measures should focus on limiting the potential for SSEs to occur by restricting large events and avoiding crowding in other public spaces. Typically, when control measures are taken, one aims to prevent transmission from those cases expected to have a high individual reproduction number, without knowing whether they will actually realize these secondary cases. Control measures thus act on the individual reproduction number and the expected transmissions, whereas their effect will be observed at the level of realized transmissions. Also, because of the increased speed at which the epidemic spreads when SSEs are present, heterogeneity could lower the doubling times [23].

Detailed contact tracing data are needed to obtain empirical offspring distributions, but these are often not available. Therefore this work should be extended such that the considered distributions can be used to infer the offspring mean R and its overdispersion from final size data [10], which are often more readily available. Rock et al. [24] have mentioned the distinction between ‘super-spreaders’ and ‘super-shedders’, who are both responsible for an above average number of secondary cases but for different reasons. In a meta-regression analysis, Chen et al. [25] investigated the relationship between the dispersion parameter k and respiratory viral load (rVL). They found that heterogeneity in rVL facilitates variation in individual infectiousness and hence may be associated with overdispersion in the number of secondary cases. Future work should thus aim to disentangle heterogeneity coming from variation in contact rates versus heterogeneity coming from variation in viral shedding. Furthermore, assuming a homogeneous Poisson process for effective contacts is likely a simplification of the real contact process, hence the use of other distributions to describe the contact process should also be considered. In a recent study, Wong & Collins [19] have suggested that the offspring distribution for SARS-CoV-2 is fat-tailed, which is consistent with a generalized Pareto distribution. We did not find evidence that a discrete Pareto offspring distribution would better describe the data used in this study (see Supplementary Text). Overall, the results of the present study suggest that, whenever possible, several distributions should be compared in terms of their fit to the observed data before making conclusions on the amount of heterogeneity by simply assuming a negative binomial offspring distribution.

Methods

Data availability

The data from Hong Kong and India used in this study are publicly available. The empirical offspring distribution from Rwanda and the code to generate data as used in the simulation study are available on GitHub (https://github.com/cecilekremer/PoiMixtSS).

Code availability

All relevant R code used in this study is available on GitHub (https://github.com/cecilekremer/PoiMixtSS).

Author contributions

NH, CF, CK, and AT contributed to the conceptualization of the study. CK, AT, NH, SB, HM, and SV conducted the literature search. CK, AT, SB, HM, and SV contributed to the analysis code; CK and AT performed the analysis. NH, CF, CK, AT, and CA contributed to the interpretation of results. CK and AT drafted the manuscript. All co-authors critically reviewed and revised the manuscript. CK and AT contributed equally to this work.

Competing interests

All authors have declared no competing interests.