Using excess deaths and testing statistics to improve estimates of COVID-19 mortalities.

Factors such as non-uniform definitions of mortality, uncertainty in disease prevalence, and biased sampling complicate the quantification of fatality during an epidemic. Regardless of the employed fatality measure, the infected population and the number of infection-caused deaths need to be consistently estimated for comparing mortality across regions. We combine historical and current mortality data, a statistical testing model, and an SIR epidemic model, to improve estimation of mortality. We find that the average excess death across the entire US is 13$\%$ higher than the number of reported COVID-19 deaths. In some areas, such as New York City, the number of weekly deaths is about eight times higher than in previous years. Other countries such as Peru, Ecuador, Mexico, and Spain exhibit excess deaths significantly higher than their reported COVID-19 deaths. Conversely, we find negligible or negative excess deaths for part and all of 2020 for Denmark, Germany, and Norway.

Introduction

The novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) first identified in Wuhan, China in December 2019 quickly spread across the globe, leading to the declaration of a pandemic on March 11, 2020 [1]. The emerging disease was termed COVID-19. As of this January 2020 writing, more than 86 million people have been infected, and more than 1.8 million deaths from COVID-19 in more than 218 countries [2] have been confirmed. About 61 million people have recovered globally.

Properly estimating the severity of any infectious disease is crucial for identifying near-future scenarios, and designing intervention strategies. This is especially true for SARS-CoV-2 given the relative ease with which it spreads, due to long incubation periods, asymptomatic carriers, and stealth transmissions [3]. Most measures of severity are derived from the number of deaths, the number of confirmed and unconfirmed infections, and the number of secondary cases generated by a single primary infection, to name a few. Measuring these quantities, determining how they evolve in a population, and how they are to be compared across groups, and over time, is challenging due to many confounding variables and uncertainties.

For example, quantifying COVID-19 deaths across jurisdictions must take into account the existence of different protocols in assigning cause of death, cataloging comorbidities [4], and lag time reporting [5]. Inconsistencies also arise in the way deaths are recorded, especially when COVID-19 is not the direct cause of death, rather a co-factor leading to complications such as pneumonia and other respiratory ailments [6]. In Italy, the clinician’s best judgment is called upon to classify the cause of death of an untested person who manifests COVID-19 symptoms. In some cases, such persons are given postmortem tests, and if results are positive, added to the statistics. Criteria vary from region to region [7]. In Germany, postmortem testing is not routinely employed, possibly explaining the large difference in mortality between the two countries. In the US, current guidelines state that if typical symptoms are observed, the patient’s death can be registered as due to COVID-19 even without a positive test [8]. Certain jurisdictions will list dates on which deaths actually occurred, others list dates on which they were reported, leading to potential lag-times. Other countries tally COVID-19 related deaths only if they occur in hospital settings, while others also include those that occur in private and/or nursing homes.

In addition to the difficulty in obtaining accurate and uniform fatality counts, estimating the prevalence of the disease is also a challenging task. Large-scale testing of a population where a fraction of individuals is infected, relies on unbiased sampling, reliable tests, and accurate recording of results. One of the main sources of systematic bias arises from the tested subpopulation: due to shortages in testing resources, or in response to public health guidelines, COVID-19 tests have more often been conducted on symptomatic persons, the elderly, front-line workers and/or those returning from hot-spots. Such non-random testing overestimates the infected fraction of the population.

Different types of tests also probe different infected subpopulations. Tests based on reverse-transcription polymerase chain reaction (RT-PCR), whereby viral genetic material is detected primarily in the upper respiratory tract and amplified, probe individuals who are actively infected. Serological tests (such as enzyme-linked immunosorbent assay, ELISA) detect antiviral antibodies and thus measure individuals who have been infected, including those who have recovered.

Finally, different types of tests exhibit significantly different “Type I” (false positive) and “Type II” (false negative) error rates. The accuracy of RT-PCR tests depends on viral load which may be too low to be detected in individuals at the early stages of the infection, and may also depend on which sampling site in the body is chosen. Within serological testing, the kinetics of antibody response are still largely unknown and it is not possible to determine if and for how long a person may be immune from reinfection. Instrumentation errors and sample contamination may also result in a considerable number of false positives and/or false negatives. These errors confound the inference of the infected fraction. Specifically, at low prevalence, Type I false positive errors can significantly bias the estimation of the IFR.

Other quantities that are useful in tracking the dynamics of a pandemic include the number of recovered individuals, tested, or untested. These quantities may not be easily inferred from data and need to be estimated from fitting mathematical models such as SIR-type ODEs [9], age-structured PDEs [10], or network/contact models [11-13].

Administration of tests and estimation of all quantities above can vary widely across jurisdictions, making it difficult to properly compare numbers across them. In this paper, we incorporate excess death data, testing statistics, and mathematical modeling to self-consistently compute and compare mortality across different jurisdictions. In particular, we will use excess mortality statistics [14-16] to infer the number of COVID-19-induced deaths across different regions. We then present a statistical testing model to estimate jurisdiction-specific infected fractions and mortalities, their uncertainty, and their dependence on testing bias and errors. Our statistical analyses and source codes are available at [17].

Methods

Mortality measures

Many different fatality rate measures have been defined to quantify epidemic outbreaks [18]. One of the most common is the case fatality ratio (CFR) defined as the ratio between the number of confirmed “infection-caused” deaths Dc in a specified time window and the number of infections Nc confirmed within the same time window, CFR = Dc/Nc [19]. Depending on how deaths Dc are counted and how infected individuals Nc are defined, the operational CFR may vary. It may even exceed one, unless all deaths are tested and included in Nc.

Another frequently used measure is the infection fatality ratio (IFR) defined as the true number of “infection-caused” deaths D = Dc + Du divided by the actual number of cumulative infections to date, Nc + Nu. Here, Du is the number of unreported infection-caused deaths within a specified period, and Nu denotes the untested or unreported infections during the same period. Thus, IFR = D/(Nc + Nu).

One major issue of both CFR and IFR is that they do not account for the time delay between infection and resolution. Both measures may be quite inaccurate early in an outbreak when the number of cases grows faster than the number of deaths and recoveries [10]. An alternative measure that avoids case-resolution delays is the confirmed resolved mortality M = Dc/(Dc + Rc) [10], where Rc is the cumulative number of confirmed recovered cases evaluated in the same specified time window over which Dc is counted. One may also define the true resolved mortality via , the proportion of the actual number of deaths relative to the total number of deaths and recovered individuals during a specified time period. If we decompose R = Rc +Ru, where Rc are the confirmed and Ru, the unreported recovered cases, . The total confirmed population is defined as Nc = Dc +Rc +Ic, where Ic the number of living confirmed infecteds. Applying these definitions to any specified time period (typically from the “start” of an epidemic to the date with the most recent case numbers), we observe that CFR ≤ M and . After the epidemic has long past, when the number of currently infected individuals I approach zero, the two fatality ratios and mortality measures converge if the component quantities are defined and measured consistently, lim CFR(t) = lim M(t) and [10].

The mathematical definitions of the four basic mortality measures Z = CFR, IFR, M, defined above are given in Table I and fall into two categories, confirmed and total. Confirmed measures (CFR and M) rely only on positive test counts, while total measures (IFR and ) rely on projections to estimate the number of infected persons in the total population N. Of the measures listed in Table I, the fatality ratio CFR and confirmed resolved mortality M do not require estimates of unreported infections, recoveries, and deaths and can be directly derived from the available confirmed counts Dc, Nc, and Rc [20]. Estimation of IFR and the true resolved mortality requires the additional knowledge on the unconfirmed quantities Du, Nu, and Ru. We describe the possible ways to estimate these quantities, along with the associated sources of bias and uncertainty below.

TABLE I:

Definitions of mortality measures.

Quantities with subscript “c” and “u” denote confirmed (i.e., positively tested) and unconfirmed populations. For instance, Dc, Rc, and Nc denote the total number of confirmed dead, recovered, and infected individuals, respectively. d((i) is the number of individuals who have died in the ith time window (e.g., day, week) of the jth previous year. The mean number of excess deaths between the periods ks and k this year is thus . Where the total number of infection-caused deaths Dc + Du appears, it can be estimated using the excess deaths over as detailed in the main text. We have also included raw death numbers/100,000 and the mean excess deaths r relative to the mean number of deaths over the same period of time from past years (see Eqs. (1)).

Subpopulation╲Measure Z	Fatality Ratios	Resolved Mortality	Excess Death Indices
Confirmed	CFR=D_cN_c	M=D_cD_c+R_c	D_eper100,000:D_c+D_u100,000
Total	IFR=D_c+D_uN_c+N_u	M=D_c+D_uD_c+D_u+R_c+R_u	relative:r=∑_i[d^(0)(i)−1J∑jJd^(j)(i)]1J∑jJ∑_id^(j)(i)

Excess deaths data

An unbiased way to estimate D = Dc + Du, the cumulative number of deaths, is to compare total deaths within a time window in the current year to those in the same time window of previous years, before the pandemic. If the epidemic is widespread and has appreciable fatality, one may reasonably expect that the excess deaths can be attributed to the pandemic [21-25]. Within each affected region, these “excess” deaths De relative to “historical” deaths, are independent of testing limitations and do not suffer from highly variable definitions of virus-induced death. Thus, within the context of the COVID-19 pandemic, De is a more inclusive measure of virus-induced deaths than Dc and can be used to estimate the total number of deaths, De ≃ Dc + Du. Moreover, using data from multiple past years, one can also estimate the uncertainty in De. In practice, deaths are typically tallied daily, weekly [21, 28], or sometimes aggregated monthly [27, 29] with historical records dating back J years so that for every period i there are a total of J + 1 death values. We denote by d((i) the total number of deaths recorded in period i from the jth previous year where 0 ≤ j ≤ J and where j = 0 indicates the current year. In this notation, , where the summation tallies deaths over several periods of interest within the pandemic. Note that we can decompose , to include the contribution from the confirmed and unconfirmed deaths during each period i, respectively. To quantify the total cumulative excess deaths we derive excess deaths per week relative to the jth previous year. Since d(0)(i) is the total number of deaths in week i of the current year, by definition . The excess deaths during week i, , averaged over J past years and the associated, unbiased variance σe(i) are given by The corresponding quantities accumulated over k weeks define the mean and variance of the cumulative excess deaths and Σe(k) where deaths are accumulated from the first to the kth week of the pandemic. The variance in Eqs. (1) and (2) arise from the variability in the baseline number of deaths from the same time period in J previous years.

We gathered excess death statistics from over 23 countries and all US states. Some of the data derive from open-source online repositories as listed by official statistical bureaus and health ministries [21–25, 29]; other data are elaborated and tabulated in Ref. [27]. In some countries excess death statistics are available only for a limited number of states or jurisdictions (e.g., Brazil). The US death statistics that we use in this study is based on weekly death data between 2015–2019 [29]. For all other countries, the data collection periods are summarized in Ref. [27]. Fig. A1(a–b) shows historical death data for NYC and Germany, while Fig. A1(c–d) plots the confirmed and excess deaths and their confidence levels computed from Eqs. (1) and (2). We assumed that the cumulative summation is performed from the start of 2020 to the current week k = K so that indicates excess deaths at the time of writing. Significant numbers of excess deaths are clearly evident for NYC, while Germany thus far has not experienced significant excess deaths.

To evaluate CFR and M, data on only Dc, Nc, and Rc are required, which are are tabulated by many jurisdictions. To estimate the numerators of IFR and , we approximate using Eq. (2). For the denominators, estimates of the unconfirmed infected Nu and unconfirmed recovered populations Ru are required. In the next two sections we propose methods to estimate Nu using a statistical testing model and Ru using compartmental population model.

Statistical testing model with bias and testing errors

The total number of confirmed and unconfirmed infected individuals Nc + Nu appears in the denominator of the IFR. To better estimate the infected population we present a statistical model for testing in the presence of bias in administration and testing errors. Although Nc + Nu used to estimate the IFR includes those who have died, depending on the type of test, it may or may not include those who have recovered. If S, I, R, D are the numbers of susceptible, currently infected, recovered, and deceased individuals, the total population is N = S + I + R + D and the infected fraction can be defined as f = (Nc + Nu)/N = (I + R + D)/N for tests that include recovered and deceased individuals (e.g., antibody tests), or f = (Nc + Nu)/N = (I + D)/N for tests that only count currently infected individuals (e.g., RT-PCR tests). If we assume that the total population N can be inferred from census estimates, the problem of identifying the number of unconfirmed infected persons Nu is mapped onto the problem of identifying the true fraction f of the population that has been infected.

Typically, f is determined by testing a representative sample and measuring the proportion of infected persons within the sample. Besides the statistics of sampling, two main sources of systematic errors arise: the non-random selection of individuals to be tested and errors intrinsic to the tests themselves. Biased sampling arises when testing policies focus on symptomatic or at-risk individuals, leading to over-representation of infected individuals.

Figure 2 shows a schematic of a hypothetical initial total population of N = 54 individuals in a specified jurisdiction. Without loss of generality we assume there are no unconfirmed deaths, Du = 0, and that all confirmed deaths are equivalent to excess deaths, so that in the jurisdiction represented by Fig. 2. Apart from the number of deceased, we also show the number of infected and uninfected subpopulations and label them as true positives, false positives, and false negatives. The true number of infected individuals is Nc + Nu = 16 which yields the true f = 16/54 = 0.27 and an IFR = 5/16 = 0.312 within the jurisdiction.

FIG. 2:

Biased and unbiased testing of a population.

A hypothetical scenario of testing a population (total N = 54 individuals) within a jurisdiction (solid black boundary). Filled red circles represent the true number of infected individuals who tested positive and the black-filled red circles indicate individuals who have died from the infection. Open red circles denote uninfected individuals who were tested positive (false positives) while filled red circles with dark gray borders are infected individuals who were tested negative (false negatives). In the jurisdiction of interest 5 have died of the infection while 16 are truly infected. The true fraction f of infected in the entire population is thus f = 16/54 and the true IFR=5/16. However, under testing (green and blue) samples, a false positive is shown to arise. If the apparent positive fraction is derived from a biased sample (blue), the estimated apparent IFR can be quite different from the true one. For a less biased (more random) testing sample (green sample), a more accurate estimate of the total number of infected individuals is when the single false positive in this sample is included, and when the false positive is excluded, and allows us to more accurately infer the IFR. Note that CFR is defined according to the tested quantities Dc/Nc which are precisely 2/9 and 2/5 for the blue and green sample, respectively, if false positives are considered. When false negatives are known and factored out CFR = 2/8 and 2/4, for the blue and green samples, respectively.

Also shown in Fig. 2 are two examples of sampling. Biased sampling and testing is depicted by the blue contour in which 6 of the 15 are alive and infected, 2 are deceased, and the remaining 7 are healthy. For simplicity, we start by assuming no testing errors. This measured infected fraction of this sample 8/15 = 0.533 > f = 0.296 is biased since it includes a higher proportion of infected persons, both alive and deceased, than that of the entire jurisdiction. Using this biased measured infected fraction of 8/15 yields IFR = 5/(0.533 · 54) ≈ 0.174, which significantly underestimates the true IFR = 0.312. A relatively unbiased sample, shown by the green contour, yields an infected fraction of 4/14 ≈ 0.286 and an apparent IFR ≈ 0.324 which are much closer to the true fraction f and IFR. In both samples discussed above we neglected testing errors such as false positives indicated in Fig. 2. Tests that are unable to distinguish false positives as negatives would yield a larger Nc, resulting in an apparent infected fraction 9/15 and an even smaller apparent IFR ≈ 0.154. By contrast, the false positive testing errors on the green sample would yield an apparent infected fraction 5/15 = 0.333 and IFR = 0.259.

Given that test administration can be biased, we propose a parametric form for the apparent or measured infected fraction to connect the apparent (biased sampling) infected fraction fb with the true underlying infection fraction. The bias parameter −∞ < b < ∞ describes how an infected or uninfected individual might be preferentially selected for testing, with b < 0 (and fb < f) indicating under-testing of infected individuals, and b > 0 (and fb > f) representing over-testing of infecteds. A truly random, unbiased sampling arises only when b = 0 where fb = f. Given Q (possibly biased) tests to date, testing errors, and ground-truth infected fraction f, we derive in the SI the likelihood of observing a positive fraction (where is the number of recorded positive tests): in which Here, μ is the expected value of the measured and biased fraction and is its variance. Note that the parameters θ = {Q, f, b, FPR, FNR} may be time-dependent and change from sample to sample. Along with the likelihood function , one can also propose a prior distribution P (θ|α) with hyperparameters α, and apply Bayesian methods to infer θ (see SI).

To evaluate IFR, we must now estimate f given and possible values for FPR, FNR, and/or b, or the hyperparameters α defining their uncertainty. The simplest maximum likelihood estimate of f can be found by maximizing with respect to f given a measured value and all other parameter values θ specified: Note that although FNRs are typically larger than FPRs, small values of f and imply that and μ are more sensitive to the FPR, as indicated by Eqs. (5) and (6).

If time series data for are available, one can evaluate the corrected testing fractions in Eq. (6) for each time interval. Assuming that serological tests can identify infected individuals long after symptom onset, the latest value of would suffice to estimate corresponding mortality metrics such as the IFR. For RT-PCR testing, one generally needs to track how evolves in time. A rough estimate would be to use the mean of over the whole pandemic period to provide a lower bound of the estimated prevalence .

The measured yields only the apparent , but Eq. (6) can then be used to evaluate the corrected which will be a better estimate of the true IFR. For example, under moderate bias |b| ≲ 1 and assuming FNR, FPR, Eq. (6) relates the apparent and corrected IFRs through .

Another commonly used representation of the IFR is . This expression is equivalent to our is defined as the fraction of infected individuals that are confirmed [30, 31]. In this alternative representation, the p factor implicitly contains the effects of biased testing. Our approach allows the true infected fraction f to be directly estimated from and N.

While the estimate depends strongly on b and FPR, and weakly on FNR, the uncertainty in f will depend on the uncertainty in the values of b, FPR, and FNR. A Bayesian framework is presented in the SI, but under a a Gaussian approximation for all distributions, the uncertainty in the testing parameters can be propagated to the squared coefficient of variation of the estimated infected fraction , as explicitly computed in the SI. Moreover, the uncertainties in the mortality indices Z decomposed into the uncertainties of their individual components are listed in Table II.

TABLE II:

Uncertainty propagation for different mortality measures.

Table of squared coefficients of variation for the different mortality indices Z derived using standard error propagation expansions [37]. We use , , , and to denote the uncertainties in the total population, confirmed cases, recoveries, and deaths, respectively. The variance of the number of excess deaths is , which feature in the IFR and . The uncertainty in the infected fraction that contributes to the uncertainty in IFR depends on uncertainties in testing bias and testing errors as shown in Eq. (A6). The term represents the covariance between Dc, Nc, and similarly for all other covariances Σe,, , , ,. Since variations in De arise from fluctuations in past-year baselines and not from current intrinsic uncertainty, we can neglect correlations between variations in De and uncertainty in Rc, Ru. In the last two rows, representing expressed in two different ways, and , respectively. Moreover, when using the SIR model to replace Du and Ru with , there is no uncertainty associated with Du and Ru in a deterministic model. Thus, covariances cannot be defined except through the uncertainty in the parameter γ = γu/μu.

Mortality Z	Uncertainties	CV^2=ΣZ2Z^2
CFR=D_cN_c	ΣD_c2, ΣN_c2, Σ_D_c,N_c	ΣD_c2Dc2+ΣN_c2Nc2−2Σ_D_c,N_cD_cN_c
IFR=D_c+D_uN_c+N_u≈D¯_efN	Σe2, ΣN2, Σ_e,N, σf2	σf2f^^2+Σe2D¯e2+ΣN2N^2−2Σ_e,ND¯_eN
M=D_cD_c+R_c	ΣD_c2, ΣR_c2, Σ_D_c,R_c	M^2(R_cD_c)^2[ΣD_c2Dc2+ΣR_c2Rc2−2Σ_R_c,D_cR_cD_c]
M=D¯_eD¯_e+R_c+R_u	Σe2, ΣR_c2, ΣR_u2, Σ_R_c,R_u	(1−M)^2Σe2D¯e2+ΣR_c2Γ^2+ΣR_u2Γ^2−2Σ_R_c,R_uΓ^2
M=D¯_eD¯_e+R_c+γ(D¯_e−D_c)	ΣR_c2, Σe2, Σ_R_c,γ, σγ2	(1−M)^2Σe2D¯e2+ΣR_c2Γ^2+(D¯_e−D_c)^2σγ2Γ^2−2(D¯_e−D_c)Σ_R_c,γΓ^2

Using compartmental models to estimate resolved mortalities

Since the number of unreported recovered individuals Ru required to calculate is not directly related to excess deaths nor to positive-tested populations, we use an SIR-type compartmental model to relate Ru to other inferable quantities [9]. Both unconfirmed recovered individuals and unconfirmed deaths are related to unconfirmed infected individuals who recover at rate γu and die at rate μu. The equations for the cumulative numbers of unconfirmed recovered individuals and unconfirmed deaths, can be directly integrated to find and . The rates γu and μu may differ from those averaged over the entire population since testing may be biased towards subpopulations with different values of γu and μu. If one assumes γu and μu are approximately constant over the period of interest, we find Ru/Du ≈ γu/μu ≡ γ. We now use , where both and Dc are given by data, to estimate and write as Thus, a simple SIR model transforms the problem of determining the number of unreported death and recovered cases in to the problem of identifying the recovery and death rates in the untested population. Alternatively, we can make use of the fact that both the IFR and resolved mortality should have comparable values and match to IFR ≈ 0.1 – 1.5% [31-33] by setting γ ≡ γu/μu ≈ 100 – 1000 (see SI for further information). Note that inaccuracies in confirming deaths may give rise to . Since by definition, infection-caused excess deaths must be greater than the confirmed deaths, we set whenever data happens to indicate to be less than Dc.

Results

Here, we present much of the available worldwide fatality data, construct the excess death statistics, and compute mortalities and compare them across jurisdictions. We show that standard mortality measures significantly underestimate the death toll of COVID-19 for most regions (see Figs. A1 and A2). We also use the data to estimate uncertainties in the mortality measures and relate them uncertainties of the underlying components and model parameters.

Excess and confirmed deaths

We find that in New York City for example, the number of confirmed COVID-19 deaths between March 10, 2020 and December 10, 2020 is 19,694 [34] and thus significantly lower than the 27,938 (95% CI 26,516–29,360) reported excess mortality cases [21]. From March 25, 2020 until December 10, 2020, Spain counts 65,673 (99% confidence interval [CI] 91,816–37,061) excess deaths [22], a number that is substantially larger than the officially reported 47,019 COVID-19 deaths [35]. The large difference between excess deaths and reported COVID-19 deaths in Spain and New York City is also observed in Lombardia, one of the most affected regions in Italy. From February 23, 2020 until April 4, 2020, Lombardia reported 8,656 reported COVID-19 deaths [35] but 13,003 (95% 12,335–13,673) excess deaths [25]. Starting April 5 2020, mortality data in Lombardia stopped being reported in a weekly format. In England/Wales, the number of excess deaths from the onset of the COVID-19 outbreak on March 1, 2020 until November 27, 2020 is 70,563 (95% CI 52,250–88,877) whereas the number of reported COVID-19 deaths in the same time interval is 66,197 [26]. In Switzerland, the number of excess deaths from March 1, 2020 until November 29, 2020 is 5,664 (95% CI 4,281–7,047) [24], slightly larger than the corresponding 4,932 reported COVID-19 deaths [35].

To illustrate the significant differences between excess deaths and reported COVID-19 deaths in various jurisdictions, we plot the excess deaths against confirmed deaths for various countries and US states as of December 10, 2020 in Fig. 3. We observe in Fig. 3(a) that the number of excess deaths in countries like Mexico, Russia, Spain, Peru, and Ecuador is significantly larger than the corresponding number of confirmed COVID-19 deaths. In particular, in Russia, Ecuador, and Spain the number of excess deaths is about three times larger than the number of reported COVID-19 deaths. As described in the Methods section, for certain countries (e.g., Brazil) excess death data is not available for all states [27]. For the majority of US states the number of excess deaths is also larger than the number of reported COVID-19 deaths, as shown in Fig. 3(b). We performed a least-square fit to calculate the proportionality factor m arising in and found m ≈ 1.132 (95% CI 1.096–1.168). That is, across all US states, the number of excess deaths is about 13% larger than the number of confirmed COVID-19 deaths.

FIG. 3:

Excess deaths versus confirmed deaths across different countries/states.

The number of excess deaths in 2020 versus confirmed deaths across different countries (a) and US states (b). The black solid lines in both panels have slope 1. In (a) the blue solid line is a guide-line with slope 3; in (b) the blue solid line is a least-squares fit of the data with slope 1.132 (95% CI 1.096–1.168; blue shaded region). All data were updated on December 10, 2020 [20, 27, 29, 36].

Estimation of mortality measures and their uncertainties

We now use excess death data and the statistical and modeling procedures to estimate mortality measures Z = IFR, CFR, M, across different jurisdictions, including all US states and more than two dozen countries.[1]. Accurate estimates of the confirmed Nc and dead Dc infected are needed to evaluate the CFR. Values for the parameters Q, FPR, FNR, and b are needed to estimate Nc+Nu = fN in the denominator of the IFR, while is needed to estimate the number of infection-caused deaths Dc + Du that appear in the numerator of the IFR and . Finally, since we evaluate the resolved mortality , through Eq. 8, estimates of , Dc, Rc, γ, and FPR, FNR (to correct for testing inaccuracies in Dc and Rc) are necessary. Whenever uncertainties are available or inferable from data, we also include them in our analyses.

Estimates of excess deaths and infected populations themselves suffer from uncertainty encoded in the variances and . These uncertainties depend on uncertainties arising from finite sampling sizes, uncertainty in bias b and uncertainty in test sensitivity and specificity, which are denoted , , and , respectively. We use Σ2 to denote population variances and σ2 to denote parameter variances; covariances with respect to any two variables X, Y are denoted as Σ. Variances in the confirmed populations are denoted , , and and also depend on uncertainties in testing parameters and . The most general approach would be to define a probability distribution or likelihood for observing some value of the mortality index in [Z, Z + dZ]. As outlined in the SI, these probabilities can depend on the mean and variances of the components of the mortalities, which in turn may depend on hyperparameters that determine these means and variances. Here, we simply assume uncertainties that are propagated to the mortality indices through variances in the model parameters and hyperparameters [37]. The squared coefficients of variation of the mortalities are found by linearizing them about the mean values of the underlying components and are listed in Table II.

To illustrate the influence of different biases b on the IFR we use from Eq. (6) in the corrected . We model RT-PCR-certified COVID-19 deaths [38] by setting the FPR = 0.05 [39] and the FNR = 0.2 [40, 41]. The observed, possibly biased, fraction of positive tests can be directly obtained from corresponding empirical data. As of November 1, 2020, the average of over all tests and across all US states is about 9.3% [42]. The corresponding number of excess deaths is , 700 [27] and the US population is about N ≈ 330 million [43]. To study the influence of variations in , in addition to , we also use a slightly larger in our analysis. In Fig. 4 we show the apparent and corrected IFRs for two values of [Fig. 4(a)] and the coefficient of variation CVIFR [Fig. 4(b)] as a function of the bias b and as made explicit in Table I. For unbiased testing [b = 0 in Fig. 4(a)], the corrected IFR in the US is 1.9% assuming and 0.8% assuming . If b > 0, there is a testing bias towards the infected population, hence, the apparent is smaller than the corrected IFR as can be seen by comparing the solid (corrected IFR) and the dashed (apparent IFR) lines in Fig. 4(a). For testing biased towards the uninfected population (b < 0), the corrected IFR may be smaller than the apparent IFR. To illustrate how uncertainty in FPR, FNR, and b affect uncertainty in IFR, we evaluate CVIFR as given in Table II.

FIG. 4:

Different mortality measures across different regions.

(a) The apparent (dashed lines) and corrected (solid lines) IFR in the US, as of November 1, 2020, estimated using excess mortality data. We set , 0.15 (black,red), FPR = 0.05, FNR = 0.2, and N = 330 million. For the corrected IFR, we use as defined in Eq. (6). Unbiased testing corresponds to setting b = 0. For b > 0 (positive testing bias), infected individuals are overrepresented in the sample population. Hence, the corrected IFR is larger than the apparent IFR. If b is sufficiently small (negative testing bias), the corrected IFR may be smaller than the apparent IFR. (b) The coefficient of variation of De (dashed line) and IFR (solid lines) with σI = 0.02, σII = 0.05, and σ = 0.2 (see Tab. II).

The first term in uncertainty given in Eq. (A6) is proportional to 1/Q and can be assumed to be negligibly small, given the large number Q of tests administered. The other terms in Eq. (A6) are evaluated by assuming σ = 0.2, σI = 0.02, and σII = 0.05 and by keeping FPR = 0.05 and FNR = 0.2. Finally, we infer Σe from empirical data, neglect correlations between De and N, and assume that the variation in N is negligible so that Σe = Σ ≈ 0. Fig. 4(b) plots CVIFR and in the US as a function of the underlying bias b. The coefficient of variation is about 1%, much smaller than CVIFR, and independent of b. For the values of b shown in Fig. 4(b), CVIFR is between 47–64% for and between 20–27% for .

Next, we compared the mortality measures Z = IFR, CFR, M, and the relative excess deaths r listed in Tab. I across numerous jurisdictions. To determine the CFR, we use the COVID-19 data of Refs. [20, 36]. For the apparent IFR, we use the representation discussed above. Although p may depend on the stage of the pandemic, typical estimates range from 4% [44] to 10% [31]. We set p = 0.1 over the lifetime of the pandemic. We can also use the apparent , however estimating the corrected IFR requires evaluating the bias b. In Fig. 5(a), we show the values of the relative excess deaths r, the CFR, the apparent IFR, the confirmed resolved mortality M, and the true resolved mortality for different (unlabeled) regions. In all cases we set p = 0.1, γ = 100. As illustrated in Fig. 5(b), some mortality measures suggest that COVID-induced fatalities are lower in certain countries compared to others, whereas other measures indicate the opposite. For example, the total resolved mortality for Brazil is larger than for Russia and Mexico, most likely due to the relatively low number of reported excess deaths as can be seen from Fig. 3 (a). On the other hand, Brazil’s values of CFR, IFR, and M are substantially smaller than those of Mexico [see Fig. 5(b)].

FIG. 5:

Mortality characteristics in different countries and states.

(a) The values of relative excess deaths r, the CFR, the with p = Nc/(Nc + Nu) = 0.1 [31], the confirmed resolved mortality M, and the true resolved mortality (using γ = 100) are plotted for various jurisdictions. (b) Different mortality measures provide ambiguous characterizations of disease severeness. (c–g) The probability density functions (PDFs) of the mortality measures shown in (a) and (b). Note that there are only very incomplete recovery data available for certain countries (e.g., US and UK). For countries without recovery data, we could not determine M and . The number of jurisdictions that we used in (a) and (c–g) are 77, 246, 73, 191, and 21 for the respective mortality measures (from left to right). All data were updated December 10, 2020 [20, 27, 29, 36].

The distributions of all measures Z and relative excess deaths r across jurisdictions are shown Fig. 5(c–g) and encode the global uncertainty of these indices. We also calculate the corresponding mean values across jurisdictions, and use the empirical cumulative distribution functions to determine confidence intervals. The mean values across all jurisdictions are (95% CI 0.0025–0.7800), (95% CI 0.0000–0.0565), (95% CI 0.0000–0.0150), (95% CI 0.0000–0.236), and (95% CI 0.000–0.193). For calculating and , we excluded countries with incomplete recovery data. The distributions plotted in Fig. 5(c–g) can be used to inform our analyses of uncertainty or heterogeneity as summarized in Tab. II. For example, the overall variance can be determined by fitting the corresponding empirical Z distribution shown in Fig. 5(c–g). Table II displays how the related can be decomposed into separate terms, each arising from the variances associated to the components in the definition of Z. For concreteness, from Fig. 5(e) we obtain which allows us to place an upper bound on using Eq. (A6), the results of Tab. II, and or on using .

Finally, to provide more insight into the correlations between different mortality measures, we plot M against CFR and against IFR in Fig. 6. For most regions, we observe similar values of M and CFR in Fig. 6(a). Althouigh we expect M → CFR and towards the end an epidemic, in some regions such as the UK, Sweden, Netherlands, and Serbia, M ≫ CFR due to unreported or incomplete reporting of recovered cases. About 50% of the regions that we show in Fig. 6(b) have an IFR that is approximately equal to . Again, for regions such as Sweden and the Netherlands, is substantially larger than IFR because of incomplete reporting of recovered cases.

FIG. 6:

Different mortality measures across different regions.

We show the values of M and CFR (a) and (using γ = 100) and with p = Nc/(Nc+Nu) = 0.1 [31] (b) for different regions. The black solid lines have slope 1. If jurisdictions do not report the number of recovered individuals, Rc = 0 and M = 1 [light red disks in (a)]. In jurisdictions for which the data indicate , we set in the denominator of which prevents it from becoming negative as long as . All data were updated on December 10, 2020 [20, 27, 29, 36].

Discussion

Relevance

In the first few weeks of the initial COVID-19 outbreak in March and April 2020 in the US, the reported death numbers captured only about two thirds of the total excess deaths [15]. This mismatch may have arisen from reporting delays, attribution of COVID-19 related deaths to other respiratory illnesses, and secondary pandemic mortality resulting from delays in necessary treatment and reduced access to health care [15]. We also observe that the number of excess deaths in the Fall months of 2020 have been significantly higher than the corresponding reported COVID-19 deaths in many US states and countries. The weekly numbers of deaths in regions with a high COVID-19 prevalence were up to 8 times higher than in previous years. Among the countries that were analyzed in this study, the five countries with the largest numbers of excess deaths since the beginning of the COVID-19 outbreak (all numbers per 100,000) are Peru (256), Ecuador (199), Mexico (151), Spain (136), and Belgium (120). The five countries with the lowest numbers of excess deaths since the beginning of the COVID-19 outbreak are Denmark (2), Norway (6), Germany (8), Austria (31), and Switzerland (33) [27][2]. If one includes the months before the outbreak, the numbers of excess deaths per 100,000 in 2020 in Germany, Denmark, and Norway are −3209, −707, and −34, respectively. In the early stages of the COVID-19 pandemic, testing capabilities were often insufficient to resolve rapidly-increasing case and death numbers. This is still the case in some parts of the world, in particular in many developing countries [45]. Standard mortality measures such as the IFR and CFR thus suffer from a time-lag problem.

Strengths and limitations

The proposed use of excess deaths in standard mortality measures may provide more accurate estimates of infection-caused deaths, while errors in the estimates of the fraction of infected individuals in a population from testing can be corrected by estimating the testing bias and testing specificity and sensitivity. One could sharpen estimates of the true COVID-19 deaths by systematically analyzing the statistics of deaths from all reported causes using a standard protocol such as ICD-10 [46]. For example, the mean traffic deaths per month in Spain between 2011–2016 is about 174 persons [47], so any pandemic-related changes to traffic volumes would have little impact considering the much larger number of COVID-19 deaths.

Different mortality measures are sensitive to different sources of uncertainty. Under the assumption that all excess deaths are caused by a given infectious disease (e.g., COVID-19), the underlying error in the determined number of excess deaths can be estimated using historical death statistics from the same jurisdiction. Uncertainties in mortality measures can also be decomposed into the uncertainties of their component quantities, including the positive-tested fraction f that depend on uncertainties in the testing parameters.

As for all epidemic forecasting and surveillance, our methodology depends on the quality of excess death and COVID-19 case data and knowledge of testing parameters. For many countries, the lack of binding international reporting guidelines, testing limitations, and possible data tampering [48] complicates the application of our framework. A striking example of variability is the large discrepancy between excess deaths De and confirmed deaths Dc across many jurisdictions which render mortalities that rely on Dc suspect. More research is necessary to disentangle the excess deaths that are directly caused by SARS-CoV-2 infections from those that result from postponed medical treatment [15], increased suicide rates [49], and other indirect factors contributing to an increase in excess mortality. Even if the numbers of excess deaths were accurately reported and known to be caused by a given disease, inferring the corresponding number of unreported cases (e.g., asymptomatic infections), which appears in the definition of the IFR and (see Tab. I), is challenging and only possible if additional models and assumptions are introduced.

Another complication may arise if the number of excess deaths is not significantly larger than the historical mean. Then, excess-death-based mortality estimates suffer from large uncertainty/variability and may be meaningless. While we have considered only the average or last values of , our framework can be straightforwardly extended and dynamically applied across successive time windows, using e.g., Bayesian or Kalman filtering approaches.

Finally, we have not resolved the excess deaths or mortalities with respect to age or other attributes such as sex, co-morbidities, occupation, etc. We expect that age-structured excess deaths better resolve a jurisdiction’s overall mortality. By expanding our testing and modeling approaches on stratified data, one can also straightforwardly infer stratified mortality measures Z, providing additional informative indices for comparison.

Conclusions

Based on the data presented in Figs. 5 and 6, we conclude that the mortality measures r, CFR, IFR, M, and may provide different characterizations of disease severity in certain jurisdictions due to testing limitations and bias, differences in reporting guidelines, reporting delays, etc. The propagation of uncertainty and coefficients of variation that we summarize in Tab. II can help quantify and compare errors arising in different mortality measures, thus informing our understanding of the actual death toll of COVID-19. Depending on the stage of an outbreak and the currently available disease monitoring data, certain mortality measures are preferable to others. If the number of recovered individuals is being monitored, the resolved mortalities M and should be preferred over CFR and IFR, since the latter suffer from errors associated with the time-lag between infection and resolution [10]. For estimating IFR and , we propose using excess death data and an epidemic model. In situations in which case numbers cannot be estimated accurately, the relative excess deaths r provides a complementary measure to monitor disease severity. Our analyses of different mortality measures reveal that

The CFR and M are defined directly from confirmed deaths Dc and suffers from variability in its reporting. Moreover, the CFR does not consider resolved cases and is expected to evolve during an epidemic. Although M includes resolved cases, its additionally required confirmed recovered cases Rc add to its variability across jurisdictions. Testing errors affect both Dc and Rc, but if the FNR and FPR are known, they can be controlled using Eq. (A3) given in the SI.

The IFR requires knowledge of the true cumulative number of disease-caused deaths as well as the true number of infected individuals (recovered or not) in a population. We show how these can be estimated from excess deaths and testing, respectively. Thus, the IFR will be sensitive to the inferred excess deaths and from the testing (particularly from the bias in the testing). Across all countries analyzed in this study, we found a mean IFR of about 0.24% (95% CI 0.0–1.5%), which is similar to the previously reported values between 0.1 and 1.5% [31-33].

In order to estimate the resolved true mortality , an additional relationship is required to estimate the unconfirmed recovered population Ru. In this paper, we propose a simple SIR-type model in order to relate Ru to measured excess and confirmed deaths through the ratio of the recovery rate to the death rate. The variability in reporting Dc across different jurisdictions generates uncertainty in and reduces its reliability when compared across jurisdictions.

The mortality measures that can most reliably be compared across jurisdictions should not depend on reported data which are subject to different protocols, errors, and manipulation/intentional omission. Thus, the per capital excess deaths and relative excess deaths r (see last column of Table I) are the measures that provide the most consistent comparisons of disease mortality across jurisdictions (provided total deaths are accurately tabulated). However, they are the least informative in terms of disease severity and individual risk, for which M and are better.

Uncertainty in all mortalities Z can be decomposed into the uncertainties in component quantities such as the excess death or testing bias. We can use global data to estimate the means and variances in Z, allowing us to put bounds on the variances of the component quantities and/or parameters.

Parts of our framework can be readily integrated into or combined with mortality surveillance platforms such as the European Mortality Monitor (EURO MOMO) project [28] and the Mortality Surveillance System of the National Center for Health Statistics [21] to assess disease burden in terms of different mortality measures and their associated uncertainty.