Using excess deaths and testing statistics to determine COVID-19 mortalities

Factors such as varied 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 from January 2020 until February 2021 is 9\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\%$$\end{document}% 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 statistically insignificant or even negative excess deaths for at least most of 2020 in places such as Germany, Denmark, and Norway. Supplementary Information The online version contains supplementary material available at 10.1007/s10654-021-00748-2.


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]. As of April 2021, an estimated 130 million people have been infected with COVID-19 disease, with more than 2.8 million deaths across 219 countries [2]. About 104 million people have recovered globally. Properly ascertaining 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. Estimating 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 [4].
For example, quantifying COVID-19 deaths across jurisdictions must take into account the existence of different protocols in assigning cause of death, cataloging co-morbidities [5], and lag-time reporting. 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 at the onset of the pandemic. In the US, 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, and the elderly. Such non-random testing may lead to overestimation 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, FPR and FNR , respectively. 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 and antibody waning are still largely unknown [9,10]. Instrumentation errors and sample contamination may also result in a considerable number of false positives and/or false negatives. Specifically, at low prevalence, Type I false positive errors can significantly bias the estimation of fatality measures.
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 [11], agestructured PDEs [12], or network/contact models [13][14][15].
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 [16][17][18] to infer the number of COVID-19-induced deaths in different regions. We then present a statistical testing model to estimate jurisdictionspecific infected fractions and mortalities, their uncertainty, and their dependence on testing bias and errors. Our statistical analyses and source codes are available at [19].

Mortality measures
Different fatality rate measures have been used to quantify epidemic outbreaks [20]. One of the most common is the case fatality ratio ( CFR ) defined as the ratio between the number of confirmed "infection-caused" deaths D c in a specified time window and the number of infections N c confirmed within the same time window, CFR = D c ∕N c [21]. Depending on how deaths D c are counted and how infected individuals N c are defined, the operational CFR may vary and may exceed one if some deaths are not included in N c .
Another frequently used measure is the infection fatality ratio ( IFR ) defined as the true number of "infection-caused" deaths D = D c + D u divided by the actual number of cumulative infections to date, N c + N u . Here, D u is the number of unreported infection-caused deaths within a specified period, and N u denotes the untested or unreported infections during the same period. Thus, IFR = D∕(N c + N u ).
One major limitation 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, especially early in an outbreak when the number of cases grows faster than the number of deaths and recoveries [12]. An alternative measure that avoids case-resolution delays is the confirmed resolved mortality M = D c ∕(D c + R c ) [12], where R c is the cumulative number of confirmed recovered cases evaluated in the same specified time window over which D c is counted. One may also define the true resolved mortality via M = D∕(D + R) , 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 = R c + R u , where R c and R u are the confirmed and unreported recovered cases, The confirmed quantities are related through the total confirmed population N c = D c + R c + I c , where I c 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 IFR ≤ M . After the epidemic has long past, when the number of currently infected individuals I approaches zero, the two fatality ratios and mortality measures converge if the component quantities are defined and measured consistently, lim t→∞ CFR(t) = lim t→∞ M(t) and lim t→∞ IFR(t) = lim t→∞ M(t) [12].
The mathematical definitions of the four basic mortality measures Z = CFR, IFR, M, M defined above are given in Table 1 and fall into two categories, confirmed and total. Confirmed measures ( CFR and M) rely only on positive test counts, while total measures ( IFR and M ) rely on projections to estimate the number of infected persons in the total population N. Of the measures listed in Table 1, 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 D c , N c , and R c [22]. Estimation of IFR and the true resolved mortality M requires the additional knowledge on the unconfirmed quantities D u , N u , and R u . We now describe the possible ways to estimate these quantities, along with the associated sources of bias and uncertainty.

Excess death data
An unbiased way to estimate D = D c + D u , 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.
Since COVID-19 exhibits appreciable fatality, one may reasonably expect that most "excess" deaths can be attributed to the pandemic [23][24][25][26][27]. Within each affected region, these excess deaths D e relative to "historical" deaths, are independent of testing limitations and do not suffer from highly variable definitions of COVID-induced death. Thus, D e is a more inclusive measure of virus-induced deaths than D c and can be used to estimate the total number of deaths, D e ≈ D c + D u . Moreover, using data from multiple past years, one can also estimate the uncertainty in D e . In practice, deaths are tallied daily, weekly [23,30], or in monthly aggregates [29,31] with historical records dating back J years so that for every period i there are a total of J + 1 values. We denote by d (j) (i) the total number of deaths recorded in period i from the j th previous year where 0 ≤ j ≤ J and where j = 0 indicates the current year. To quantify the total cumulative excess deaths we define d where deaths and their variances are accumulated from the first to the k th week of the pandemic. The variance in Eqs. (1) and (2) arise from the variability in the number of deaths from the same time period across J previous years.
We gathered excess death statistics from over 75 countries and all US states. Some of the data derive from open-source online repositories as listed by official statistical bureaus and health ministries [23][24][25][26][27]31]; other data are elaborated and tabulated in Ref. [29]. In some countries excess death statistics are available only for a limited number of states or jurisdictions (e.g., Turkey and Indonesia). The US death statistics that (2) Table 1 Definitions of mortality measures. Quantities with subscript "c" and "u" denote confirmed (i.e., positively tested) and unconfirmed populations. For instance, D c , R c , and N c denote the total number of confirmed dead, recovered, and infected individuals, respectively. d (j) (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 D e within all periods i = 1, … , k is thus The total number of infection-caused deaths D c + D u can be estimated using the excess deaths D e 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 Eq. (1)]

Subpopulation
Measure Z

Fatality Ratios
Resolved Mortality Excess Death Indices we use in this study is based on weekly death data between 2015-2019 [31]. For all other countries, the data collection periods are summarized in Ref. [29]. Fig. 1(a-b) shows historical death data for New York City (NYC) and Germany, while Fig. 1(c-d) plots the associated 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 second week of February 2021 (week number k = 7 ). Significant numbers of excess deaths are clearly evident for NYC, while Germany had not experienced significant excess COVID-19 deaths in the first half of 2020.
To evaluate CFR and M, only data on D c , N c , and R c are required, which are are tabulated by many jurisdictions. We can approximate the numerators of the IFR and M by using, for example, the mean excess deaths D e ≈ D c + D u defined in Eq. (2). For the denominators, estimates of the unconfirmed infected N u and unconfirmed recovered populations R u are required. In the next two sections we propose methods to estimate N u using a statistical testing model and R u using a compartmental population model.

Statistical testing model with bias and testing errors
Testing bias and different sources of uncertainty in disease testing confound the estimation of the true prevalence of a disease. If not corrected, a sampling bias due to preferential testing of symptomatic individuals, healthcare workers, and certain high-risk groups [33] may lead to an overestimation of the fraction or total number of infected individuals. Since the total number of confirmed and unconfirmed cases, N c + N u , appears in the denominator of the IFR , we develop a statistical model for its estimation in the presence of testing errors and bias in administration of tests.
Although N c + N u used to estimate the IFR includes those who have died, 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 = (N c + N u )∕N = (I + R + D)∕N for tests that include recovered and deceased individuals (e.g., antibody tests), or f = (N c + N u )∕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 N u is mapped onto the problem of identifying the true fraction f of the population that has been infected. Thus, we need to estimate f from testing results. 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, D u = 0 , and that all confirmed deaths are equivalent to excess deaths, so that D e = D c = 5 . Apart from the number of deceased, we also show the number Fig. 1 Examples of seasonal mortality and excess deaths. The evolution of weekly deaths in (a) NYC (over seven years) and (b) Germany (over six years) derived from data in Refs. [28,29]. Grey solid lines and shaded regions represent the historical numbers of deaths and corresponding confidence intervals defined in Eq. (1). Blue solid lines indicate weekly deaths, and weekly deaths that lie outside the confidence intervals are indicated by solid red lines. The red shaded regions represent statistically significant mean cumulative excess deaths D e . The reported weekly confirmed deaths d (0) c (i) (dashed black curves), reported cumulative confirmed deaths D c (k) (dashed dark red curves), weekly excess deaths d e (i) (solid grey curves), and cumulative excess deaths D e (k) (solid red curves) are plotted in units of per 100,000 in (c) and (d) for NYC and Germany, respectively. The excess deaths and the associated 95% confidence intervals given by the error bars are constructed from historical death data in (a-b) and defined in Eqs. (1) and (2). In NYC there is clearly a significant number of excess deaths that can be safely attributed to COVID-19, while in the first half of 2020 there had been no significant excess deaths in Germany. Excess death data from other jurisdictions are shown in the SI and typically show excess deaths greater than reported confirmed deaths [with Germany an exception as shown in (d)] of true infected and uninfected subpopulations and label them as true positives, false positives, and false negatives. The true number of infected individuals is N c + N u = 16 which yields the true fraction f = 16∕54 = 0.296 and an IFR = 5∕16 = 0.312.
Also shown in Fig. 2 are two examples of sampling. Biased sampling and testing are depicted by the blue contour in which 6 of the 15 individuals are alive and infected, 2 are deceased, and the remaining 7 are healthy. For simplicity, we start by assuming no testing errors. The measured infected fraction of this sample, f b = 8∕15 = 0.533 > f = 0.296 , is biased since it includes a higher proportion of infected persons than that of the entire jurisdiction. Using this biased measured fraction yields the apparent IFR = 5∕(0.533 × 54) ≈ 0.174 , which significantly underestimates the true IFR ≈ 0.312 . A less-biased sample, shown by the green contour, yields an infected fraction of 4∕14 ≈ 0.286 and an apparent IFR = 5∕(0.286 × 54) ≈ 0.324 which are much closer to the true fraction f = 0.296 and the true, "correct" IFR ≈ 0.312.
In both samples discussed above we neglected testing errors such as false positives, as indicated in Fig. 2. Tests that are unable to distinguish false positives from true negatives would yield a larger N c , resulting in an apparent infected fraction f b = 9∕15 and an even smaller apparent IFR ≈ 0.154 , as in the blue sample. By contrast, the false positive testing errors in the green sample would yield an apparent infected fraction f b = 5∕15 ≈ 0.333 and IFR ≈ 0.259 . Given a true infected fraction f, we now derive the probability of measuring the value f b under biased testing and testing errors rates FPR and FNR.
We begin by proposing a parametric form for the apparent or measured (under biased testing) infected fraction Equation (3)   Here, is the expected value of the measured and biased fraction f b , and 2 T is its variance. The parameters associated with testing are denoted T = {Q, b, FPR, FNR} , and may be time-dependent and may change from sample to sample.
Our formulae Eqs. (3), (4), and (5) provide relationships among the variables and testing parameters. In the testing context, Q, Q + , and thus f b =Q + ∕Q will be measured and used to infer f, which we have explicitly separated out in Eq. (4) above. The goal is to use f b as an estimate of f b and ultimately estimate the true infected fraction f through Eq. (3). This requires a-priori knowledge of the parameters T . For example, FPR and FNR are typically given by test specifications. To estimate b instead, one must conduct tests on a small truly random sample and compare results to those from a more widely used biased one. A summary of the main variables that we use in our statistical testing model is given in Table 2.
For a given measured value f b , s simple maximum likelihood estimate of f can be found by maximizing Eq. (4) with respect to f, with all other parameters in T specified, Additional correction terms of O(1∕ √ Q) are neglected. Note that although FNR s are typically larger than FPR s, small values of f and f b imply that f and are more sensitive to the FPR , as indicated by Eqs. (5) and (6).
If time series data for f b are available, one can evaluate the estimated 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 f would suffice to estimate corresponding mortality metrics such as the IFR . For RT-PCR testing, one generally needs to track how f b evolves in time. A rough estimate would be to use the mean of f b over the whole pandemic period to provide a lower bound of the estimated prevalence f .
The fraction f b measured under biased testing yields only the apparent IFR = D e ∕(f b N) or expected apparent IFR =D e ∕(f b N) , but Eq. (6) can then be used to estimate the true, corrected IFR ≈D e ∕(f N) . For example, under moderate bias |b| ≲ 1 and assuming FNR , FPR , f b ≲ 1 , Eq. (6) can be used to relate the expected true IFR to the expected Another commonly used representation of the IFR is is defined as the fraction of infected individuals that are confirmed [34,35]. 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 Q + and N.
While the estimate f needed to evaluate IFR 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 . Although statistical models for inferring T from data can be similarly constructed, here, we use a simple linear approximation to propagate uncertainty in the testing parameters to the squared coefficient of variation 2 f ∕f 2 of the estimated infected fraction f . The uncertainties in the mortality indices Z decomposed into the uncertainties in their respective individual components are listed in Table 3.

Using compartmental models to estimate resolved mortalities
Since the number of unreported or unconfirmed recovered individuals R u required to calculate the total resolved mortality M is not directly related to excess deaths nor to positive-tested populations, we use an SIR-type compartmental model to relate R u to other inferable quantities [11]. Both unconfirmed recovered individuals and unconfirmed deaths D u are related to unconfirmed infected individuals I u who recover at rate u and die at rate u . The equations for the cumulative numbers of unconfirmed recovered individuals and unconfirmed deaths, Recorded positive tests Infected fraction under biased testing False negative rate can be directly integrated to find R u (t) = ∫ t 0 u (t � )I u (t � )dt � and D u (t) = ∫ t 0 u (t � )I u (t � )dt � . 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 R u ∕D u ≈ u ∕ u ≡ . We now use D u ≈ D e − D c to estimate R u ≈ (D e − D c ) and write M as An estimate for the expected value M can be obtained by substituting the excess deaths D e in Eq. (8) with D e from historical data.
Thus, a simple SIR model transforms the problem of determining the number of unreported death and recovered cases in M 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 M should have comparable values and match M to IFR ≈ 0.1 − 1.5 % [35-37] by setting ≡ u ∕ u ≈ 100 − 1000 (see SI for further information). Note that inaccuracies in confirming deaths may give rise to D c >D e . Since by definition, infection-caused excess deaths must be greater than the confirmed deaths, we set D e − D c = 0 whenever data happens to indicate D e to be less than D c .

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. 1, A1, and A2). We also use the data to estimate uncertainties in the mortality measures and relate them to uncertainties of the underlying components and model parameters.
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 March 30, 2021 in Fig. 3. We observe in Fig. 3a that excess deaths in countries like Egypt, Mexico, Peru, Russia, and South Africa are significantly higher than the corresponding number of confirmed COVID-19 deaths. In particular, they were about three times higher than the number of reported COVID-19  [22,29,31,40] in Egypt, Peru, Russia, and South Africa. 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. 3b. We performed a least-square fit to calculate the proportionality factor m arising in D e = mD c and found m ≈ 1.087 (95% CI: 1.052-1.121). That is, across all US states, the number of excess deaths is about 9% larger than the number of confirmed COVID-19 deaths in the period from January 2020 until February 2021.

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, M across different jurisdictions, including all US states and more than 75 countries. 1 Accurate estimates of the confirmed infecteds N c and the confirmed dead D c are needed to evaluate the CFR . Values for the parameters Q, FPR , FNR , and b are needed to estimate N c + N u = f N in the denominator of the IFR , while D e is needed to estimate the number of infection-caused deaths D c + D u that appear in the numerator of the IFR and M . Finally, since we evaluate the resolved mortality M , through Eq. (8), estimates of D e , D c , R c , , and FPR , FNR (to correct for testing inaccuracies in D c and R c ) 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 Σ 2 e and 2 f . The uncertainty in f depends on the uncertainties arising from finite sampling sizes, uncertainty in bias b and uncertainty in test sensitivity and specificity, which are denoted 2 b , 2 I , and 2 II , respectively. We use Σ 2 to denote variances in populations and 2 to denote variances in intrinsic parameters; covariances with respect to any two variables X, Y are denoted as Σ X,Y . Variances in the confirmed popula- and also depend on uncertainties in testing parameters 2 I and 2 II . The most general approach would be to define a probability distribution or likelihood for observing a value of the mortality index Z ∈ [z, z + d z] . As outlined in the SI, these probabilities depend on the probability densities of the components of the mortalities, which in turn may depend on hyperparameters that define these probability densities. Here, we simply assume uncertainties that are propagated to the mortality indices through variances in the model parameters and hyperparameters [41]. 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 3.
To illustrate the influence of different biases b on the IFR we use f from Eq. (6) in the corrected IFR ≈ D e ∕(f N) . We model RT-PCR-certified COVID-19 deaths [42] by setting the FPR = 0.05 [43] and the FNR = 0.2 [44,45]. The observed, possibly biased, fraction of positive tests f b can be directly obtained from corresponding empirical data.
As of February 6, 2021, the average of f b over all tests and across all US states is about 9.3% [46]. The corresponding number of excess deaths is D e = 536,617 [29] and the US population is about N ≈ 330 million [47]. To study the influence of variations in f b , in addition to f b = 0.089 , we also use a slightly larger f b = 0.15 in our analysis. In Fig. 4 we use the D e from the US and show the expected apparent and corrected IFR s for two values of f b [Fig. 4a] and the coefficient of variation CV IFR [Fig. 4b] as a function of the bias b. For unbiased testing [ b = 0 in Fig. 4a], the corrected IFR in the US is 3.6% assuming f b = 0.089 and 1.4% assuming f b = 0.15 . If b > 0 , a higher proportion of the infected population is tested, hence, the expected apparent IFR =D e ∕(f b N) is smaller than the true IFR , as can be seen by comparing the solid (corrected IFR ) and the dashed (apparent IFR ) lines in Fig. 4a. 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 CV IFR as given in Table 3.
The first term in uncertainty 2 f ∕f 2 given in Eq. (A6) in the SI 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 b = 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 using Eq. (2), neglect correlations between D e and N, and assume that the variation in N is negligible so that Σ e,N = Σ N ≈ 0 . Fig. 4b plots CV IFR and CV D e in the US as a function of the underlying bias b. The coefficient of variation CV D e is about 1%, much smaller than CV IFR , and independent of b. For the values of b shown in Fig. 4b, CV IFR is between 51-68% for f b = 0.089 and between 20-27% for Next, we compared the mortality measures IFR , CFR , M, M , and r listed in Table 1 across numerous jurisdictions.
To determine the CFR , we use the COVID-19 data of Refs. [22,40]. For the expected apparent IFR , we use the representation IFR = pD e ∕N c discussed above. Although p may depend on the stage of the pandemic, typical estimates range from 4% [48] to 10% [35]. We set p = 0.1 over the lifetime of the pandemic. We can also directly use IFR =D e ∕(f N) , however estimating the corrected IFR requires evaluating the bias b. In Fig. 5a, we show the running values (up until March 30, 2021) of the relative excess deaths r , the CFR , the apparent IFR , the confirmed resolved mortality M, and the true resolved mortality M across all regions. In all cases we set p = 0.1, = 100 . As illustrated in Fig. 5b, 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 expected total resolved mortality M 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. 3a. On the other hand, Brazil's values of CFR , IFR , and M are substantially smaller than those of Mexico [see Fig. 5b].
Here, the CI measures the variation across jurisdictions. To calculate ⟨M⟩ and ⟨M⟩ , we excluded countries with incomplete recovery data. The distributions plotted in Fig. 5c-g can be used to inform our analyses of uncertainty or heterogeneity as summarized in Table 3. For example, the overall variance Σ 2 Z can be determined by fitting the corresponding empirical Z distribution shown in Fig. 5c-g. Table 3 displays how the related CV 2 Z can be decomposed into separate terms, each arising from the variances associated to the components in the definition of Z. For concreteness, from Fig. 5e we obtain CV 2 IFR = Σ 2 IFR ∕⟨IFR⟩ 2 ≈ 8.05 which allows us to place an upper bound on 2 b using Eq. (A6), the results of Table 3, and or on 2 . Finally, to provide more insight into the correlations between different mortality measures, we plot M against CFR and M against IFR in Fig. 6. For most regions, we observe similar values of M and CFR in Fig. 6a. Although we expect M → CFR and M → IFR towards the end of an epidemic, in some jurisdictions such as the UK, the Netherlands, and Sweden, M ≫ CFR due to unreported or incomplete reporting of recovered cases. About 30% of the regions that we show in Fig. 6b have an IFR that is approximately equal to M . Again, for regions such as Sweden and the Table 3 Uncertainty propagation for different mortality measures. Table of squared coefficients of variation CV 2 = Σ 2 Z ∕Z 2 for the different mortality indices Z derived using standard error propagation expansions [41]. Here, we have used the expected values IFR and M and decomposed the CVs through variances about the mean values of all parameters such as D e . We use Σ 2 to denote the uncertainties in the total population, confirmed cases, recoveries, and deaths, respectively. The variance of the number of excess deaths is Σ 2 e , which feature in the IFR and M . The uncertainty in the infected fraction 2 f that contributes to the uncertainty in IFR depends on uncertainties in testing bias and testing errors as shown in Eq. (A6). The term Σ D c ,N c represents the covariance between D c , N c , and similarly for all other covariances Σ e,N , Σ D c ,R c , Σ R c ,R u , Σ R c , . Since variations in D e arise from fluctuations in past-year baselines and not from current intrinsic uncertainty, we can neglect correlations between variations in D e and uncertainty in R c , R u . The last two rows represent M expressed in two different ways, Γ ≡D e + R c + R u and D e + R c + (D e − D c ) , respectively. Moreover, when using the SIR model to replace D u and R u with D e − D c ≥ 0 , there is no uncertainty associated with D u and R u in a deterministic model. Thus, covariances cannot be defined except through the uncertainty in the parameter Netherlands, M is substantially larger than IFR because of incomplete reporting of recovered cases.

Discussion
In this paper, we review running COVID mortality metrics starting from the onset of the pandemic through December 2020, or later dates where data were available. 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 [17]. 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 [17].
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 ten countries with the largest numbers of excess deaths since the beginning of the COVID-19 outbreak (all numbers per 100,000) are Peru (358) [29]. From this point of view, as of March, 2021, these countries have not experienced statistically significant higher mortality due to COVID-19.
The proposed use of excess deaths in standard mortality measures may provide a more meaningful estimate 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, specificity, and sensitivity.
Underlying our use of excess deaths D e for evaluating the disease IFR and mortality M is the assumption that behavioral changes during the pandemic (social distancing, maskwearing, lockdowns) had no appreciable effect on death. For example, the mean traffic deaths per month in Spain between 2011-2016 is about 174 persons [49], so any pandemic-related changes to traffic volumes would have little impact considering the much larger number of COVID-19 deaths. However, other deaths due to, e.g., increased suicides and deferred medical treatment may contribute significantly to D e . 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 [50]. 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 [17], increased suicide rates [51], 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 M (see Table 1), is challenging and only possible if additional models and assumptions are introduced.
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., , the underlying error in the determined number of excess deaths can be estimated using historical death statistics from the same jurisdiction. Uncertainties  Table 3) in mortality measures can also be decomposed into the uncertainties of their component quantities, including that of the positive fraction f that in turn depend on uncertainties in the testing parameters. While we have considered only the average or last values of f b , our framework can be straightforwardly extended and dynamically applied across successive time windows, using e.g., Bayesian or Kalman filtering approaches.
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 [52] complicates the application of our framework. A striking example of variability is the large discrepancy between mean excess deaths D e and confirmed deaths D c across many jurisdictions which render mortalities that rely on D c suspect.
Finally, we have not partitioned the excess deaths or mortalities into subpopulations in age or other attributes such as sex, co-morbidities, occupation, etc. 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 M defined in Table 1 may provide different characterizations of disease severity across jurisdictions due to differences in testing bias and reporting protocols. The propagation of uncertainty and coefficients of variation that we summarize in Table 3 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 . 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 M . The number of jurisdictions that we used in (a) and (c-g) are 136, 247, 127, 191, and 75 for the respective mortality measures (from left to right). All data from January 2020 onwards were included and last updated on March 30, 2021 [22,29,31,40] mortality measures are preferable to others. If the number of recovered individuals is being monitored, the resolved mortalities M and M should be preferred over CFR and IFR , since the latter suffer from errors associated with the time-lag between infection and resolution [12]. For estimating IFR and M , 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 D c 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, it also requires confirmed recovered cases R c , adding to its variability across jurisdictions. Testing errors affect both D c and R c , 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.44% (95% CI: 0.0-2.0%), which is similar to the previously reported values between 0.1 and 1.5% [35][36][37]. • In order to estimate the true resolved mortality M , an additional relationship is required to estimate the unconfirmed recovered population R u . In this paper, we propose a simple SIR-type model in order to relate R u to measured excess and confirmed deaths through the ratio of the recovery rate to the death rate. The variability in reporting D c across different jurisdictions generates uncertainty in M 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 1) should 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 M 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 [30] and the Mortality Surveillance System of the National Center for Health Statistics [23] to assess disease burden in terms of different mortality measures and their associated uncertainty.

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s10654-021-00748-2. Data availability All datasets used in this study are available from Refs. [23-27, 29, 31]. The CDC excess death numbers and CIs are based on an overdispersed Poisson generalized linear model [31]. For the excess death data of Ref. [29], we determined CIs using the linear regression model class in R. The source codes used in our analyses are publicly available at [19].
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.