Monitoring Health Systems by Estimating Excess Mortality

2.1 Mortality Data

For Puerto Rico, we requested individual-level mortality information with no personal identifiers from the Department of Health of Puerto Rico and obtained individual records including date, gender, age, and cause of death from January 1985 to May 2020. We also obtained daily death counts from Florida, New Jersey, and Louisiana’s Vital Statistic systems from January 2015 to December 2018, January 2007 to December 2015, and January 2003 to December 2006, respectively. Further, we obtained mortality data made public on May 29, 2020, by the CDC[12]. This dataset included weekly estimates of death counts for each state for the period of January 2017 and May 2020. Finally, we obtained individual-level mortality data from Cook County, IL, medical examiner office[13], which included demographic information for each death including age, sex, and race.

2.2 Population Estimates

Yearly population estimates for Puerto Rico were obtained from the Puerto Rico Statistical Institute. Yearly population estimates for other US states were obtained from the US Census bureau. To obtain an estimate of the population displacement after Hurricane María in Puerto Rico, we use population proportion estimates from Teralytics – a company that works with governments and private clients to assess human movement by partnering with mobile network operators. This company provided daily population proportion estimates relative to an undisclosed baseline from May 31, 2017, to April 30, 2018. This was based on all subscribers of a major, undisclosed, telecommunications company in Puerto Rico. Due to structural damage post-Hurricane Maria, proportion estimates for the four weeks after the storm were not provided. We imputed missing data with linear interpolation of observed values. We then computed population estimates by multiplying the US census population estimate on July 19, 2017, Teralytics baseline date, times the aforementioned proportions, and obtained smooth values using local regression (Supplementary Figure 1).

• Download figure
• Open in new tab

Figure 1:

Estimated event effects as percent increase over expected mortality. A) Comparison of six hurricanes. B) Estimated increase and point-wise confidence intervals during the Chikungunya epidemic in Puerto Rico. C) Estimated increase and confidence interval for the the United States from January 2017 through May 2020.

2.3 Statistical model

Extending the work developed for monitoring outbreaks [14, 15, 16, 17] we modeled death counts with the following mixed model with μ_t the expected number of deaths at time t for a typical period, 100 × f(t) the percent increase at time t due to an unusual event, ε_t a time series of auto-correlated random variables representing natural variability, and T the total number of observations. The expected counts μ_t can be further decomposed into with N_t the population at time t, α(t) a slow moving trend that accounts for changes such as the improved health outcomes we have observed during the last 20 years (Supplementary Figure 2A), s(t) a yearly periodic function representing a seasonal trend (Supplementary Figure 2B), and w(t) is a day of the week effect. We refer to f as the event effect and assume f (t) = 0 for typical periods not affected by a natural disaster or outbreak. Our model permits non-continuities in f (t) to account for sudden direct effects caused by natural disasters. Note that we can conveniently represent excess deaths at time t as μ_t × f (t) and compute excess deaths for a time period by just adding these up for the t’s in the interval of interest. We assume the vector [ε₁, …, ε _T]^T follows a multivariate normal distribution with mean 1 and variance-covariance matrix Σ determined by an auto-regressive (AR) process of order p. In practice, obtaining Maximum Likelihood Estimates (MLE) for this model is not straight-forward, in particular because the flexibility in f (t) and Σ makes the model non-identifiable in practice. To overcome this challenge, we implemented a three step approach that works well in practice, as demonstrated both by simulation and empirical validation. Further details and validation of our model and estimation procedure is provided in the Supplementary Material. The code to reproduce our analysis as well as the R package implementing the methods are available at https://github.com/rafalab/excess-mortality-paper and https://github.com/rafalab/excessmort, respectively.

• Download figure
• Open in new tab

Figure 2:

Cumulative excess mortality estimates. A) Daily excess death estimates and 95% confidence intervals for events in Puerto Rico starting on September 20, 2017, September 21, 1998, July 14, 2014, and March 1, 2020 for Hurricane María, Hurricane Georges, the Chikungunya epidemic, and the COVID-19 pandemic, respectively. B) Estimated cumulative excess mortality and reported COVID-19 cumulative deaths in the United States. C) State-specific cumulative excess death estimates from March 1, 2020 to May 16, 2020, the period associated with the COVID-19 pandemic, plotted against the excess mortality estimates from December 10, 2017 to February 24, 2018, the period associated with the 2018 flu season. The identity line is represented in red and separates states in which the increase in mortality was higher during COVID-19 from those in which the increase in mortality was higher during the 2018 flu season.

3.1 Detecting and quantifying indirect effects after hurricanes

Using our approach, we found that the death rate in Puerto Rico increased by 73% (95% CI: 58% to 87%) the days after hurricane María made landfall. The increased death rate persisted until March 2018 (Figure 1A, Supplementary Figure 3). To understand if such a large and sustained increase is common after the passage of a hurricane, we examined mortality counts in Louisiana after Katrina, New Jersey after Sandy, Florida after Irma as well as two other hurricanes in Puerto Rico, Hugo in 1989 and Georges in 1998. Although not as severe, a similar pattern to María was observed for Georges in 1998 when the death rate increased by 36% (95% CI: 24% to 48%) the days after landfall and an increased death rate persisted until December 1998 (Figure 1A, Supplementary Figure 3). The effects of Katrina on Louisiana were much more direct. On August 29, 2005, the day the levees broke, there were 834 deaths, which translates to an increase in death rate of 689% (not shown in Figure 1A). However, the increase in mortality rates for the months following this catastrophe was substantially lower than in Puerto Rico after María and Georges. None of the other hurricanes examined had direct or indirect effects nearly as high as in Puerto Rico (Figure 1A, Supplementary Figure 3).

• Download figure
• Open in new tab

Figure 3:

Comparison of estimated percent increase in mortality in Cook county, Illinois, for individuals above 40 between whites and blacks with corresponding 95% confidence intervals. A) 40 to 59 age group. B) 60 to 74 age group. C) 75 and over age group.

3.2 Quantifying the effects of epidemics

Our approach can also be used to detect and quantify the effects of epidemics or outbreaks. As an example, we fit our model to Puerto Rico mortality data from 1985 to 2020 and, apart from the hurricane seasons mentioned in the previous section, we detected an unusual increase in mortality rates from July 2014 to February 2015. This period coincides with the 2014-2015 Chikungunya outbreak [18, 19] (Figure 1B). The effects were particularly strong for individuals over 60 years and for the 0 to 4 years age group. (Supplementary Figure 4).

• Download figure
• Open in new tab

Figure 4:

Cause-specific mortality index. The solid-black line corresponds to the percent change in proportion of mortality, , attributable to a cause of interest. The shaded area represents a 95% confidence interval for . A) Bacterial infections mortality index defined as ICD10 codes between A00 and A79 for the dates around Hurricane María. B) Respiratory diseases mortality index defined as ICD10 codes between J00 and J99 for the dates around the COVID-19 pandemic in Puerto Rico.

Although mortality records for the period of the COVID-19 pandemic are not yet complete, the CDC has implemented a correction to the most recent counts and is making weekly state-level data from January 2017, to May 2020, publicly available [12]. We fit our model to these data for each state and our event effect estimates detect the pandemic for several states, as well as the 2018 flu epidemics (Supplementary Figure 6). By aggregating the results for the entire United States we confirm that the COVID-19 pandemic had a substantially larger impact than the 2018 flu season (Figure 1C).

3.3 Excess mortality estimates

Once we have fit our model, we can compute excess mortality estimates for any period as well as standard errors that take into account sampling and natural variation. We used the estimated parameters from the Puerto Rico data to compute excess deaths for the 365 days following hurricanes Georges, and Maria. We also computed these quantities for the 365 days following the start of the detected effects from the 2014 Chikungunya outbreak. Finally, we computed excess deaths in Puerto Rico during the COVID-19 pandemic, specifically from March 1, 2020, to April 15, 2020, the last day for which the death registry has complete data (typically, it takes 90 days for 99% of deaths to be registered).

For María the excess death curve grew for over six months after landfall, with a point estimate on March 20, 2018, six months after landfall, of 3,200 (95% CI: 2,800 to 3,600), and for Georges, the excess death curve grew for over three months after landfall, with a point estimate on December 21, 1998, three months after landfall, of 1,200 (95% CI: 950 to 1,500) (Figure 2A). The period of July 14, 2014, to February 20, 2015, associated with a Chikungunya outbreak was also identified to have unusually high mortality rates. However, we also noted a decrease in cumulative deaths for several months after February 20, 2015, consistent with a harvesting effect [20, 21]. A year after the start of the Chikungunya epidemic, on July 14, 2015, we observe a point estimate for excess mortality of 975 (95% CI:400 to 1,550) (Figure 2A). We note that from March 1, 2020, to April 15, 2020, associated with the COVID-19 pandemic, no effect is detected.

Next, using the CDC data we computed excess mortality for the period of the COVID-19pandemic for which we have data, specifically for the week ending on March 15, 2020, to the week ending on May 16, 2020. When we compared these numbers to the reported numbers of COVID-19 deaths [22], we found a difference of 32,257 (95% CI: 30,211 to 34,303) deaths (Figure 2B, Table 1). Further, we found that excess deaths in the US were already above 100,000 on May 9, 2020. We also note that in the US, the number of excess deaths during the COVID-19 pandemic is much larger than the 2018 flu season, although this is not true for all states (Figure 2C, Table 1, Supplementary Figure 5). We also note that 58% of the cases come from the four northeast states: New York, New Jersey, Massachusetts, and Pennsylvania (Supplementary Figure 6). We note that Connecticut was not included in the analysis due to missing data.

3.4 Across group comparisons

At the time of writing, the Cook County Medical Examiner’s Office was making up-to-date records available online including demographic information for each individual. This office investigates any human death in Cook County that falls within the category of diseases constituting a threat to public health, among other circumstances[13]. If we consider to be the expected number of deaths recorded by the Medical Examiner’s Office, we can fit our model and interpret f (t) as the percent increase in deaths examined by this office. We fit the model to these data and observed the effects of the COVID-19 pandemic, in particular among individuals over 75. However, a persistent increase in mortality above normal levels is observed for individuals above 40 (Supplementary Figure 7). Furthermore, we found a difference in the effect between white and black individuals in March and April (Figure 3). Among those between 40 and 59 years of age, the percent increase in medical examiner death rates for Cook county was as high as 188% (95% CI: 150% to 226%) and 101% (95% CI: 68% to 134%) on April 23 for blacks and whites, respectively. Furthermore, for the 60 to 74 age group we observed a percent increase of 628% (95% CI: 560% to 696%) and 398% (95% CI: 345% to 451%) on the same date, respectively. Finally, the percent increase for individuals 75 years and older on the same date was 2,193% (95% CI: 1,996% to 2,389%) and 1,062% (95% CI: 974% to 1,150%), respectively, another statistically significant finding (Figure 3). The population size, in this database, for other racial groups was substantially smaller and were not included in the analysis.

3.5 Natural Variability and Correlated Counts

Exploration data analysis after fitting a standard Poisson model demonstrates that the daily data measurements are correlated (Supplementary Figure 8). Assuming independence, when in fact variables are correlated, leads to underestimates of the standard error of sums or averages. Cumulative excess mortality is the sum of excess deaths during a period of interest, therefore accounting for correlation allows our approach to better estimate the uncertainty of this quantity. To demonstrate this, we compared our model to a standard Poisson model and an over-dispersed Poisson model both assuming independent outcomes[14, 15, 16, 17]. As an illustrative example, we fit each of the three models to the Puerto Rico data for individuals 75 and over, and randomly selected 100 intervals of sizes L = 10, 50, and 100 days from periods with no events and computed the total number of deaths in each interval S_l = ∑_t∊L Y_t. We then used the fitted models to estimate the expected value Ê(S_l) and standard error and SÊ (S_l) and form a z-statistic Z_l = [S_l − Ê (S_l)]/SÊ(S_l). The Central Limit Theorem predicts that, if the expected value and standard errors are estimated correctly, Z_l is approximately normal with expected value 0 and standard error 1. We found that the Poisson and over-dispersed Poisson model underestimates the variance of this sum while modeling the correction improves our standard error estimates (Supplementary Figure 9).

3.6 An index based on relative proportion of specific causes of death

Limitations of our approaches for detecting and monitoring health crises in real-time are that it requires 1) complete death count data and 2) an estimate of population displacement. Note that incomplete death counts or underestimates of population displacement will result in an underestimate of the death rate. However, we if we know that a specific cause of death is associated with the event in question, we can construct an useful index that does not have these limitations. Specifically, if we redefine Y_t as the number of deaths due to a specific cause at time t and N_t as the total number of recorded deaths at time t, then the estimate of f can be interpreted as the percent increase in the proportion of deaths due to that cause, even if the data is incomplete. The underlying cause of death for the Puerto Rico mortality data was coded using the tenth version of the International Classifications of Diseases (ICD10) codes[23]. Because deaths due to bacterial infections have been reported to increase after natural disasters in locations with poor public health infrastructure[24, 25], we defined an index based on ICD10 codes between A00 and A79 as bacterial infections. We found that deaths attributable to this cause were elevated after hurricane María in Puerto Rico (Figure 4A) pointing to the possibility that this index could be used in real-time to detect challenges to the health system. During the COVID-19 pandemic, if under-reporting is occurring, one would expect causes related to the respiratory system to increase. A similar index for respiratory diseases defined as ICD10 codes between J00 and J99 computed for the COVID-19 pandemic shows no anomalies (Figure 4B).