Data

Provisional death counts for 2019–2020 were regularly published on the CDC website (https://www.cdc.gov/nchs/covid19/covid-19-mortality-data-files.htm) on Wednesdays at 5 p.m. Reported deaths were categorised by Morbidity and Mortality Weekly Report (MMWR) week of publication and by US jurisdiction where the death occurred. For the current study, 22 snapshots were collected with publication dates between September 2020 and the first week of February 2021. One snapshot published the week of 16 September 2020 (MMWR week 38) was omitted for technical reasons, which was not critical for the study. The time period containing the most recent week with non-zero deaths covered MMWR week 34 (week ending date: 22 August 2020) in the earliest collected snapshot through week 53 of 2020 (week ending date: 2 January 2021) in the last four snapshots. The death counts for the week of publication and for the preceding week as well as for the weeks of 2021 were likely to be missed in any given snapshot because of zero reported counts; all non-zero death counts less than 10 were masked by the CDC for privacy reasons. The reporting jurisdictions included the 50 states with New York state separated into two jurisdictions: New York City and the rest of the state. Additionally, the District of Columbia and Puerto Rico were among the total 53 jurisdictions.

Historical records of weekly deaths from 2014 to 2018 were retrieved from the same source as the provisional counts for 2019–2020. The structure of the dataset was analogous except that it was not subject to any changes in the future. The jurisdictional counts of reported COVID-19 deaths were assessed via the daily trends published by the CDC (https://covid.cdc.gov/covid-data-tracker/#trends_dailytrendscases (accessed 26 February 2021)).

Reporting delay: parametric estimation (independent and partial pool model)

The reporting delay distribution describes the distribution of time periods between the occurrence of an event and its reporting to the system. The probability distribution function of the reporting delay is usually modelled using one of three unimodal distributions with positive support (f _i(○; θ),  i = 1, 2, 3): the gamma, Weibull or lognormal distributions. The set θ consists of two parameters: the mean and the standard deviation (s.d.) of the reporting delay distribution.

To estimate the reporting delay, the death count $d_w^{( {s, j} ) }$ reported on week w by jurisdiction j in any of the earlier snapshots s = 1, …, (S − 1) was compared with the death count $d_w^{( {S, j} ) }$ reported in the latest snapshot S [Reference Akhmetzhanov11, Reference Tsuzuki12]. Poisson likelihood was used to infer the unknown parameters θ _j:

(1)$$L^{(j)} \left(\theta_j\semicolon \; \left\{d_w^{(s, j)} \right\} \right) = \mathop \prod \limits_{s = 1}^{S-1} \mathop \prod \limits_{\matrix{\small_{w \lt T(s)} \cr d_w^{(s, j)} \ne 0}} {\rm Poisson}_{\rm pmf} \left(d_w^{(S, j)} \semicolon \;E \left[d_w^{(S, j)} \right] \right),$$

(2)$$E \big [ {d_w^{( {S, j} ) } } \big ] = d_w^{( {s, j} ) } \cdot \displaystyle{{F( {T( S ) -w + 0.5\semicolon \;\theta_j} ) -F( {0.5\semicolon \;\theta_j} ) } \over {F( {T( s ) -w + 0.5\semicolon \;\theta_j} ) -F( {0.5\semicolon \;\theta_j} ) }}, \;$$

for any w, j and s = 1, …, (S − 1). The second equation accounts for the continuity factor [Reference White and Held6]. Here, T(s) denotes the publishing times of the snapshots s, and Poisson_pmf(d;E[d]) is the probability mass function:

$${\rm Poisso}{\rm n}_{{\rm pmf}}( {d\semicolon \;E[ d ] } ) = \displaystyle{{{( {E[ d ] } ) }^d\cdot\exp\! ( {-E[ d ] } ) } \over {d!}}.$$

To estimate variation in reporting delays across jurisdictions, two different approaches were employed. In the first approach, the reporting delay for each jurisdiction was estimated independently, such that each likelihood L ^(j) (j = 1, …, 53) was maximised with respect to θ _j. In the second approach, a partial pool model was used to infer the common shared mean reporting delay and its s.d. [Reference Flaxman22–Reference Gelman24]. In the latter context, the reporting delays for various jurisdictions were closely related to each other, sharing a common mean value μ. Any deviations from the shared value of the mean were modelled using a Student's t-distribution:

(3)$$m_j\;\sim \;{\rm Student}\hbox{\underscore}{t}( {\nu , \;\mu , \;\sigma_\mu } ) , \;$$

where the other two parameters were the degree of freedom ν and the standard error of the mean $\sigma _\mu$. The Student's t-distribution (3) was chosen over the normal distribution because it is less sensitive to outliers. Like the nonparametric estimation described below, the first approach showed promise when used to nowcast the number of deaths that have yet to be reported. The second approach was used to identify the common mean μ, and the corresponding P-values (percentiles of the Student's t-distribution) for detecting outliers (i.e. jurisdictions significantly deviating from others in their reporting delays).

A negative binomial distribution could have been used instead of the Poisson distribution in equation (1). However, simulations showed that the value of the overdispersion parameter in the negative binomial distribution approached an arbitrarily large value, implying equivalence of the negative binomial and Poisson likelihood functions. A similar conclusion was reached in another relevant study [Reference McGough15].

Mixture model

Although the reporting delay distribution was chosen from one of three unimodal distributions, this selection induced a constraint to the modelling framework by imposing a structural prior. Another approach to account for all three distributions within a single model is to consider mixtures of distributions. This strategy provides a greater degree of flexibility because each distribution contributes to the total likelihood proportionally to the relative weights π _i ($\sum\nolimits_i {\pi _i} = 1$), subject to the data fit. By contrast with the common practice in formulating mixture models, where each component distribution has its own set of parameters (e.g. each of the three distributions would have their own means and s.d.s), we assumed that all distributions shared the same set of parameters (mean and s.d.). This ensures a higher convergence probability of implemented Markov Chain Monte-Carlo simulations [Reference Keller and Kamary25].

Alternatively, the best-fit distribution could be selected based on information criteria (e.g. the widely applicable information criterion (WAIC) or ‘leave-one-out’ information criterion (LOOIC) [Reference Gelman24, Reference Watanabe26]). However, integrating out unobserved (latent) variables from the model, such as the death counts masked by CDC, can be challenging [Reference Watanabe27]. The mixture model implements all three component distributions based on relative weights π _i. In this case, there was no need to integrate out latent variables or to manually calculate likelihoods for each data point, as is required using other methods [28].

Following these assumptions, the total likelihood for the mixture model was defined as follows:

$$L_{\rm \Sigma }( {\{ {\theta_j} \} \semicolon \;\{ {w, \;s} \} } ) = \mathop \prod \limits_{\,j = 1}^J \mathop \sum \limits_{i = 1}^3 \pi _i^{( j ) } L_i^{( j ) } ( {\theta_j\semicolon \;\{ {w, \;s} \} } ) , \;$$

where $\pi _i^{( j ) }$ are relative weights ($\sum\nolimits_i {\pi _i^{( j ) } } = 1$) and the subscript i indicates one of the three distributions (i = 1, 2, 3). The component likelihoods $L_i^{( j ) }$ are given by equation (1) and the expected deaths $E[ {d_w^{( {S, j} ) } } ] _i$ (an internal argument of the likelihoods) are given by (2) respective to each distribution i:

(4)$$E[ {d_w^{( {S, j} ) } } ] _i = d_w^{( {s, j} ) } \cdot\displaystyle{{F_i( {T( S ) -w + 0.5\semicolon \;\theta_j} ) -F_i( {0.5\semicolon \;\theta_j} ) } \over {F_i( {T( s ) -w + 0.5\semicolon \;\theta_j} ) -F_i( {0.5\semicolon \;\theta_j} ) }}, \;$$

where F _i denotes the cumulative distribution function of the component distribution i. The posterior probability for each component distribution $q_i^{( j ) }$ could be then determined using the equation:

$$q_i^{( j ) } = \displaystyle{{\pi _i^{( j ) } L_i^{( j ) } ( {\theta_j\semicolon \;\{ {w, \;s} \} } ) } \over {\mathop \sum \nolimits_{i = 1}^3 \pi _i^{( j ) } L_i^{( j ) } ( {\theta_j\semicolon \;\{ {w, \;s} \} } ) }}.$$

Reporting delay: nonparametric estimation

For nonparametric estimation of the reporting delay, the reverse-time discrete hazard was defined as previously described [Reference Lawless7, Reference Höhle and der Heiden8]: g _j(d) = Pr(delay = d | delay ≤ d) = f _j(d)/F _j(d). Here, the variable d was introduced such that a zero value (d = 0) corresponds to the death count reported within the first 2 weeks (equivalently, within the first 10 days because all snapshots were published on Wednesdays rather than on the last day of the week). Other values (d = 1, 2, …, D) correspond to reporting delays of weeks, respectively. The upper bound D denotes the maximum delay, implying that F _j(delay ≥ D) = 1. Finally, g _j(0) = 1 was imposed, and other hazards g _j(d > 0) were found by fitting the probability distribution functions $f_j( d ) = g_j( d ) \prod\nolimits_{i = d + 1}^D {( {1-g_j( i ) } ) }$ to the data. Equations (1), (4) were used, accounting for the only difference in defining the cumulative distribution functions:

(5)$$F_j( {T( s ) -w\semicolon \;\{ {g_j( d ) } \} } ) = \left\{{\matrix{ {\mathop \prod \limits_{i = T( s ) -w-1}^D ( {1-g_j( i ) } ) \;\;\;{\rm if}\;\;\;T( s ) -w-1 \le D, \;} \cr {{\rm 1\;\;\;otherwise, \;}} \cr } } \right.$$

where the parameter D was set to 20 weeks in the simulations.

Nowcasting procedure

To predict the number of deaths not yet reported by the surveillance system, a prospective nowcasting framework was applied [Reference Akhmetzhanov11, Reference Cowling29]. The number of yet unreported deaths on a given week was sampled from a negative binomial distribution that followed the failure-counting interpretation [Reference Champredon30]. The first parameter of the negative binomial distribution (the number of ‘failures’) was the number of already reported deaths during that week, whereas the second parameter (the probability of ‘success’) was the cumulative distribution function of the reporting delay counted from week of death, w, to the publication date of the latest snapshot.

Expected excess mortality

The expected weekly number of deaths was estimated using a Poisson linear regression model [Reference Weinberger21] involving a seasonal component but neglecting to adjust for severe influenza and associated pneumonia. The posterior median and the 95% upper bound were set as two thresholds. The range of differences between the nowcasted number of deaths and each of these thresholds was then reported as excess deaths as in previous studies [5, Reference Kawashima20]. All negative differences were assigned to zero. The reader is referred to the Appendix for additional mathematical details of the statistical framework.

Technical details

To infer individual mean reporting delays and perform nowcasting, nonparametric estimation of the reporting delay distribution was used. A parametric estimation was implemented only for verification purposes (Fig. 1). A partial pool model was used to calculate the common mean of the reporting delay shared across jurisdictions. Because of excessive computational time requirements, only a lognormal distribution was implemented in the partial pool model. The choice of the lognormal distribution was guided by its dominant selection while fitting the mixture models for various jurisdictions (Appendix Fig. 1).

Fig. 1. (a) Mean reporting delay by jurisdiction using different estimation approaches (legend). Error bars indicate the 95% CI for individually estimated reporting delays using a parametric model. Dashed line indicates a common mean delay inferred from the partial pool model. The entire shaded area indicates the 95% CI for the common mean delay, whereas the dark shaded area covers the interquartile range of the posterior. (b) Relationship between the fraction of deaths reported within the first 10 days and the mean reporting delay by jurisdiction obtained from non-parametric estimation of the reporting delay distribution. Dashed line indicates an estimate of 61% cited in technical notes of CDC [5]. (c) Correlation between number of reported COVID-19 deaths per 100 000 from September to December 2020 and the mean reporting delay by jurisdiction. Solid line is obtained from a linear regression model. Shaded area indicates 95% CI.

Statistical inference was conducted within the Bayesian framework realised in CmdStan (version 2.26, https://mc-stan.org). Pre- and post-processing of the data and results were performed in the Python environment (version 3.8). The code snippets are available at http://github.com/aakhmetz/Excess-mortality-in-US-2020.