Over the last year since hurricanes Irma and Maria struck Puerto Rico, a number of analyses have attempted to estimate excess mortality attributable to the storms, which few doubted was higher than the 64 official deaths the government had certified.1–8 These have generated a great deal of discussion concerning how such estimates are made and, by extension, their validity. One line of criticism, as the director of the Federal Emergency Management Agency recently stated, is that it is hard to know what is accurate because “the numbers are all over the place.”9 In this commentary, we discuss differences in each analysis in terms of methods employed and assumptions made. We show that, despite these differences, all analyses are generally consistent in their estimation of the mortality impact of these storms.
Because we will never observe a world in which the hurricanes did not occur, estimating excess mortality including direct and indirect deaths requires a counterfactual estimate of what mortality would have been had the disaster not occurred, then comparing this to the actual number of observed deaths over the same period. In standard vital registration-based analyses, differences in point estimates produced by these can be attributed to differences in (1) the time period over which excess deaths are evaluated; (2) the length of the counterfactual period; (3) adjustments made for changes to the population at risk of mortality over time, including, importantly in the case of the hurricanes that struck Puerto Rico, population displacement; (4) adjustments made for long and short term (e.g., seasonal) variation in mortality; and (5) the specification of the functional form of the relationship between factors identified in (3) and (4).
The Table presents all prior publicly available estimates of excess mortality after hurricane Maria, as well as replications of the registration-based estimates using published methodologies and the most recent vital registration data. The replicated estimates eliminate potential differences owing to the completeness or quality of mortality data used. For the one survey-based estimate in the Table, we present results from a re-analysis of the data, correcting for bias in estimation of the mortality rate.6 The GW Institutional Review Board exempted the project as it posed minimal risk. The data and analytic files for the replication and re-analyses can be found at https://prstudy.publichealth.gwu.edu/datasets. The top panel of the Table shows estimates that do not take population displacement after the storms into account, the bottom those that do. The most obvious source of variation is the period over which the number of excess deaths is evaluated. We have provided a separate column for each period.
ESTIMATES NOT ACCOUNTING FOR POPULATION DISPLACEMENT AFTER THE STORM
The first three estimates presented in the Table1,4,5 were produced in the weeks immediately after hurricane Maria with aggregate vital registration data either publicly available or provided by the Puerto Rico Vital Statistics Registry. All compare daily numbers of deaths in a specified period after the storm to the average number of deaths in a comparable period before the storm. In the Centro de Periodismo Investigativo (CPI) analysis, 985 excess deaths were estimated to have occurred in the 40 days after Maria, and 1065 over the months of September and October 2017, compared with the same periods in 2016.4 Using this same methodology, but with the most recently updated registration data, we estimate slightly higher numbers for both periods.
Limiting the counterfactual period to the prior year has the advantage that long-term trends in mortality and the population at risk have limited influence on estimated differences. The disadvantage, however, is that random variation in mortality from year to year has a potentially greater impact. In the New York Times (NYT) analysis, the counterfactual estimate is based on average daily deaths between September 20 and October 31 from 2 years prior, 2015 and 2016.1 This mitigates somewhat the impact of both long-term trends and random annual variation, although these may still be large. They estimate 1052 excess deaths in this period. Replicating their analysis with the updated registration data, we again estimate a slightly higher number.
Rivera and Rolke5 estimated 822 excess deaths from September 20th to October 31st, calculating average daily deaths in the 42 days after the storm. Their counterfactual period is restricted, however, to the 19 days of September before the storm. The disadvantage to this strategy is two-fold. First, it assumes that mortality rates stay constant from month to month. Puerto Rico generally experiences an increase in deaths from September to October, and we could expect this to have happened in 2017 as well. Second, it magnifies the potential for random variation in daily death rates to influence the counterfactual estimate, potentially a reason their excess mortality estimate is the lower bound relative to others over the same period.
Santos-Lozada and Howard2 take a longer-term approach to the counterfactual. They begin by estimating 95% confidence intervals for average monthly deaths in September and October over the period 2010–2016. They then compare observed mortality in September and October 2017 to the upper bound of the counterfactual interval for each month. This yields estimates of 518 excess deaths in September and 567 in October (Santos-Lozada 1 in the Table). Updating their analysis, the authors estimate 459 excess deaths in September, 564 in October, and 116 in November (Santos-Lozada 2 in the Table).3 Using an identical methodology (without imputing October 2014 deaths as they do), our replicated estimates are extremely close to those reported in their first analysis. Although this approach has the advantage of controlling for long-term variation in seasonality relative to other analyses discussed so far, it increases the likelihood that long-term trends in both the level and variability of mortality as well as changes in the structure of the population at risk over the period will influence the counterfactual estimates.
None of the registration-based analyses discussed thus far attempted to adjust their counterfactual estimates for differences in the distribution of the population at risk. The more recent George Washington University (GW) analysis, commissioned by the Government of Puerto Rico, having the advantage of complete, individual vital registration records, does. Counterfactual mortality is modeled using data from September 2017 to February 2018, adjusting for long-term changes in mortality, seasonality, and the distribution of the population by age, sex, and socioeconomic status as well as interactions between these.8 One set of models makes no adjustment for population displacement after the storm (labeled the “census” scenario), whereas another does (labeled the “displacement scenario”). Despite these differences, under the census scenario, excess mortality estimated from the preferred model is quite close to those of virtually all others.
The main conclusion thus far is that, despite substantial differences in methodology, these estimates are remarkably consistent. Omitting the Rivera and Rolke5 estimate based on the highly constrained counterfactual period, estimates of excess mortality from September to October 2017 ranged between 1023 and 1077 deaths. Secondarily, analyses limiting the counterfactual period to the prior 1 or 2 years tended to produce slightly higher estimates than those which employed longer counterfactual periods, likely because of differences in the treatment of seasonality.
ESTIMATES ACCOUNTING FOR POPULATION DISPLACEMENT AFTER THE STORM
Thus far, we have only considered analyses that do not take population displacement into account. Puerto Rico, however, experienced massive out-migration after September 2017, substantially lowering the population at risk of mortality.10–12 Two analyses based on vital-registration data have adjusted the estimated population at risk to account for population displacement.7,8 A third, based on a household survey of mortality after the storm, does so implicitly, estimating the population at risk directly.6 The GW estimate is derived using the same methodology described above, additionally adjusting for cumulative displacement in each month. Excess mortality from September to October is estimated to be 1271 (95% CI: 1154, 1383), from September to December 2098 (95% CI: 1872, 2315), and from September to February 2975 (95% CI: 2658, 3290).8
The more recent registration-based estimate by Acosta and Irizarry7 takes a similar, but in some ways quite different approach. Although they do not model changes to the population at risk, they do adjust for population size and displacement. Daily mortality for the period January 2015 to June 2018 is modeled using functional forms for long-term and short-term seasonal trends close to those used by GW. In contrast to that analysis, however, they explicitly parameterize variation not associated with these trends to estimate excess death rates directly. Despite these important methodologic differences, they estimate 3433 (95% CI: 3189, 3676) excess deaths from 20 September 2017 to 15 April 2018, consistent with the GW estimates and a persistent but diminishing effect of the storms after February.
The final analysis discussed here takes a completely different approach to estimating excess mortality. In a widely reported article, Kishore et al.,6 using household survey data, estimated 4645 excess deaths from 20 September through 31 December 2017. They also report an estimate of 5045 excess deaths, an 8.6% increase over their primary estimate, after post-hoc adjustment for potential survivor bias attributable to single-individual households.6 These estimates are much larger than those presented by GW and Acosta and Irizarry, which cover longer periods. This is primarily because of bias introduced in their calculation of the post-Maria mortality rate. The reported estimates are derived by aggregating the total number of deaths observed in the sample and dividing this by their estimate of population size. When this is done, every death, instead of every household rate, is given equal weight. To estimate an unbiased rate using their sampling methodology, the count of household deaths needs to be conditioned on (or divided by) household size. Our re-analysis of their data, unweighted for sampling methodology and omitting adjustment for single-individual households as their primary estimate did, finds that conditioning on household size as an offset produces an estimate of 2285 excess deaths (95% CI: −816, 5387) using their population estimate (which is almost identical to that in the GW analysis) and assumed baseline mortality rate. Weighting for the stratified sampling methodology, which is appropriate, we estimate 2022 excess deaths (95% CI: −3796, 7839). Adjusting this estimate upward for potential single-individual household bias by the same proportion reported in their analysis, we would arrive at an estimate of 2196 excess deaths. This unbiased estimation method produces point estimates of excess mortality that are nearly identical to the GW estimate for the period September to December, although they are not statistically different from zero. We further show in our re-analysis code that one can arrive at a similar unbiased estimate (albeit with a much larger confidence interval) by randomly selecting one individual per household for analysis. Survey analyses such as this are invaluable in emergency or disaster contexts where, unlike Puerto Rico, no complete vital registration system exists. Large sampling variation associated with estimates from them, and sample size requirements for minimizing it must be considered in sample design, however.
Far from being “all over the place”, this summary has shown that, considering differences in the time period excess deaths are calculated over, all estimates produced thus far are remarkably consistent with each other. Minor differences in the registration-based estimates may be attributable to differences in counterfactual periods used, adjustments for changes to the population at risk of mortality by such factors as age and sex, adjustments for seasonality and long-term mortality trends, specifications of functional form of these factors and estimation strategies. Even in the case of the Kishore et al. analysis, a completely different methodology produces almost identical estimates to those produced with registration data once bias in the published estimate is corrected.
The major difference in these estimates has been shown to be owing to whether a particular analysis accounted for population displacement after the storm. This is highlighted by the GW study that produces estimates consistent with others under both the census and displacement scenarios. Because migration away from the island did occur, we believe estimates accounting for displacement are the most accurate and that the emerging scientific consensus produced by this work is that between September 2017 and late winter to early spring 2018, somewhere between 3000 and 3500 deaths in Puerto Rico may be attributed to the effects of hurricanes Maria and Irma.
ABOUT THE AUTHORS
JOHN SANDBERG is an Associate Professor in GW’s Department of Global Health, specializing in demographic methods. His current research concerns the association between social network processes and health outcomes. CARLOS SANTOS-BURGOA is a Professor in GW’s Department of Global Health. He was Dean of the School of Public Health of Mexico at the National Institute of Public Health, Director General at Mexico´s Ministry of Health and Senior Advisor and Acting Department Director at the Pan American Health Organization. AMIRA ROESS is an Assistant Professor in GW’s Department of Global Health, trained as an epidemiologist. Her current research concerns preparedness and response and the ecology and evolution of emerging diseases, including environmental, social, behavioral, and economic drivers of emerging diseases. ANN GOLDMAN-HAWES is an epidemiologist and economist whose research includes the design and implementation of protocols for cost analysis and cost effectiveness studies of population health intervention initiatives. CYNTHIA M. PEREZ is a Professor at the University of Puerto Rico School of Public Health, Department of Biostatistics and Epidemiology. Her current work involves the evaluation of the association between oral microbiota and HPV infection and the longitudinal assessment of the bidirectional relationship between periodontitis and glucose homeostasis. ALEJANDRA GARCIA-MEZA is a Research Associate in GW’s Department of Global Health. She currently also works as a consultant for the World Bank and the Pan American Health Organization. LYNN GOLDMAN is the Dean of the Milken Institute School of Public Health, a former Assistant Administrator for Toxic Substances in the Environmental Protection Agency, a member of the National Academy of Medicine, and serves on the National Academy of Medicine Council, as well as the governing board of the National Academy of Sciences.