Excess deaths and Hurricane Mar\'ia

We clarify the distinction between direct and indirect effects of disasters such as Hurricane Mar\'ia and use data from the Puerto Rico Vital Statistics System to estimate monthly excess deaths in the immediate aftermath of the hurricane which struck the island in September of 2017. We use a Bayesian linear regression model fitted to monthly data for 2010--16 to predict monthly death tallies for all months in 2017, finding large deviations of actual numbers above predicted ones in September and October of 2017 but much weaker evidence of excess mortality in November and December of 2017. These deviations translate into 910 excess deaths with a 95 percent uncertainty interval of 440 to 1,390. We also find little evidence of big pre-hurricane mortality spikes in 2017, suggesting that such large spikes do not just happen randomly and, therefore, the post-hurricane mortality spike can reasonably be attributed to the hurricane.

ogy for estimating excess deaths to months before the hurricane actually happened. It is impossible for the hurricane to cause excess deaths before it even existed so we hope and expect that the differences between observed and predicted values for the pre-hurricane months will cluster around 0, aside from some random fluctuation. This procedure is analogous to testing whether improvements in health that followed patients' consumption of a new drug were truly caused by the drug by determining whether similar improvements followed the consumption of a placebo in a control group. We find that no other month between February and August of 2017 displays differences between observed and predicted values at all close in magnitude to those for September and October, although January does display a rather large difference. These null results for the pre-hurricane period of February through August and, to a lesser extent, for January are reassuring because they suggest that huge prediction errors do not tend to just happen randomly and for no reason.
Thus, they give us confidence that our excess-deaths findings for September and October of 2017 truly reflect reality rather than being mere artefacts of random fluctuations.
There exists a fairly large and growing literature on excess deaths attributed to wars (e.g. Roberts et al. (2004); Burnham et al. (2006); Coghlan et al. (2006); Degomme and Guha-Sapir (2010), strong claims for the accuracy of such estimates (Wise, 2017) and even proposals to incorporate estimates of excess war deaths into the Sustainable Development Goals (SDG)framework (Alda and McEvoy, 2017). Pre-war death rates have provided the typical counterfactual baseline for these estimates although one prominent estimate used, unconvincingly (see Human Security Report Project (2011)), the death rate for the whole of Sub-Saharan Africa (Coghlan et al., 2006) as a counterfactual baseline. These estimates are plagued by exceptionally wide error bands, sometimes approaching plus or minus 100 percent (Roberts et al., 2004;Burnham et al., 2006;Coghlan et al., 2006), and a methodological defect that conflates violent deaths caused directly by war with non-violent deaths indirectly attributable to war (Spagat andvan Weezel, 2017, 2018b,a). High-profile studies have been further discredited due to unrepresentative samples (Coghlan et al., 2006;Human Security Report Project, 2011) and fabricated data (Burnham et al., 2006;Spagat, 2010). Yet recent surveys of estimates of excess war deaths take an uncritical (Alda and McEvoy, 2017), even laudatory (Wise, 2017), view while ignoring the uncertainty and underplaying the controversy that surrounds these estimates. Some strong claims in the literature, e.g., that non-violent war deaths dwarf violent ones (Wise, 2017) or that excessdeath point estimates should be incorporated into SDG's (Alda and McEvoy, 2017), are seriously undermined by a realistic appraisal of the uncertainty surrounding excess wardeath estimates.
For several reasons, the case of Hurricane Mara provides a good learning opportunity for the war-deaths literature in addition to the intrinsic importance of understanding the event itself. First, in contrast to many war settings for which death rates must be estimated through surveys, we have solid, if not perfect, month-by-month data on deaths in Puerto Rico both before and after the hurricane. Yet we can also observe a quasi experiment in which a methodology typical of that used in the war-deaths literature was implemented in Puerto Rico by Kishore et al. (2018) who used a survey to estimate both pre and post hurricane death rates before reliable post hurricane data became available. The results of Kishore et al. (2018) was a central estimate well above what we know now to be a reasonable number and an extremely wide uncertainty interval, the latter of which is a hallmark of the war death literature. Second, the counterfactual thinking appropriate for assessing excess deaths in Hurricane Mara is relatively simple compared to war-based counterfactuals. This is true, despite the fact that the treatment of migration after the hurricane has a potentially large effect on excess death estimates several months after the event, as Santos-Burgoa et al. (2018) point out. Nevertheless, realistic war counterfactuals are likely to also include migration in addition to further complicating factors such as prewar droughts, coups or economic sanctions. Thus, we might expect estimates of excess deaths resulting from Hurricane Mara to provide an upper bound on the quality attainable for estimates of excess deaths resulting from wars.
The paper proceeds as follows. First, we briefly clarify the distinction between deaths caused directly by the hurricane and deaths caused indirectly by the hurricane. Next we present our model and estimates. We then compare our methodology with those of some of the other excess deaths estimates. Finally, we argue for the application of our methods for measuring excess deaths in future violent events.
It is important to distinguish between deaths caused directly by the hurricane and deaths caused indirectly by the hurricane. If the hurricane blew down a tree which then struck a man and killed him that would be a direct death. Such deaths are relatively easy to identify and tallying them up is a perfectly valid exercise. Yet it is unlikely that direct deaths will capture the full lethal effect of the hurricane. It is plausible, for example, that just after the hurricane X out of every 100 heart attack victims died while just before the hurricane only Y out of 100 died where X > Y . If so then we can speak of an "excess" of heart attack victims due to the hurricane. And, of course, we can extend this excess deaths concept to encompass deaths from all immediate causes, not just heart attack deaths. When post-hurricane death rates systematically exceed pre-hurricane rates then we can view the difference as hurricane-caused excess. All the above-quoted studies do this.
We should recognize that our estimates of excess deaths are in the aggregate and not at the individual level. Excess death models, including our own, do not predict whether any particular death is an 'excess' death but rather whether, in the case of Puerto Rico, the death count increased above previous expectation after the hurricane. Consider, for example, a particular heart attack victim who died before reaching a hospital amidst a shortage of ambulance drivers during the hurricanes aftermath. Would this man have survived if the hurricane had never happened? There may be no clear and convincing answer to this question. This man may have died anyway, even in the best of circumstances. Or maybe he could have been saved if an ambulance had arrived unusually quickly and contained an unusually experienced paramedic. There can be many relevant factors governing the fates of particular people and even analysts with good inside information, such as the GWU team, will struggle to identify particular individuals whose fate flipped from survival to death because of, and only because of, the hurricane. Thus, a big advantage to the statistical approach to measuring excess deaths is that it operates in terms of probabilities rather than certainties as is appropriate for an environment within which yes-no answers are elusive. In contrast, the intentions of many people in Puerto Rico to create a list of all the people killed by the hurricane (see Sutter and Shoichet (2018)  However, the PRVSS has stopped releasing data and is not responding to emails according to Santos-Lozada (personal communication). Meanwhile, the PRIS released the most recent figures in the face of a possibly debilitating reorganization (Guglielmi, 2018;Wade, 2018).
We have no original insight into the behind-the-scenes statistical drama in Puerto Rico which is important but which has no practical effect on any of the excess deaths estimates made using the PRVSS data. Nevertheless, it is worth noting that the PRVSS numbers used by Santos-Lozada and Howard (2018b,a) in their analysis are extremely close, but do not exactly match, the numbers released by PRIS, even for early months. We did not repeat our analysis using the earlier PRVSS numbers because the PRVSS never released an October 2017 number: Santos-Lozada and Howard (2018b,a) impute the October tally in their study. These data discrepancies must be resolved before there can be a fully definitive account of hurricane-related deaths in Puerto Rico.  (table 1). The data also shows that the number of deaths in January 2017 -well before the hurricane struck the island -is also significantly above average.
To gain better insight into the potential impact of the hurricane on mortality we standardize the time-series data to account for historical variation. Specifically, for each observation x t we subtract that month's mean µ m during the baseline period (2010-16) and divide by that month's standard deviation to obtain σ m :

Empirical framework
Bayesian methods are applied to analyse the data and estimate the number of excess deaths.
The model can be denoted in terms of a normal likelihood function as: where y t is the actual reported number of deaths at time t and x t indicates the month of the year. α is a constant which, in this case, represents the estimated number of deaths in Notes: Standard deviation in parentheses.
Source: Santos-Lozada and Howard (2018a,b); Varela (2018) the baseline month (January) while θ m is the estimated difference in deaths relative to the baseline (i.e. January) for each individual month (Feb-Dec). In other words, we regress the number of deaths (y 1 , ..., y n ) on a vector of month indicators (x 1 , ..., x n ) omitting the baseline month so that the model is not overdetermined. Due to the time-series nature of the data the errors (σ 2 y ) are adjusted for serial correlation using an AR(1) correction; modeling the error as t = ρ t−1 + η t .
The model's key feature is that the predicted (or fitted) number of deaths y i depends on the month of the year through the estimated parameter distribution θ m . As such it accounts for systematic differences between months due to, for instance, seasonality. The model is fitted to the data covering the years 2010-16 and the estimated posterior distribution is subsequently used to generate predictions for the expected number of deaths for each month in 2017, including the pre-hurricane months which serve to check the predictions against the data. The predicted values for September-December serve as a counterfactual for a hypothetical world in which the hurricane never happened. For the pre-hurricane months (January-August) the predicted values serve as a placebo test; the expectation is that the predictions are close to the historical ranges subject to some random variation.
In short, the excess mortality estimates are based on deviations in post-hurricane monthly tallies from the observed historical patterns for each month.
One constraint the analysis faces is that there is, unfortunately, only a limited amount of information to construct a baseline. The available data only goes back to 2010 and, in any case, the further we go back in time the less useful the data potentially are for predicting the future. Given the relatively small sample size (N = 84) it is important to account for uncertainty. An advantage of the Bayesian framework is that it efficiently incorporates all the available information, including in small samples (Zellner, 1988). Importantly, this framework treats the data as fixed and the parameters random. Thus, under the samplingbased Bayesian method used in this study (Hamiltonian Monte Carlo), the quality of the inference is controlled not by asymptotic properties but by the number of samples taken (Kruschke, 2014).
Ideally in a small sample an informative prior is specified as an initial estimate (Dunson, 2001;Kadane, 2015;McNeish, 2016). However, given that all the available data is needed to estimate the parameters, this means that no information is left to formulate informative priors . Therefore, we use a diffuse prior distribution that relies on the default settings in R's 'brms' package (Bürkner, 2017) -an interface for Stan (Carpenter et al., 2017) -which follow a Student t distribution (t (3, 0, 1)). This entails that the estimates will be close to what would be obtained using frequentist methods and, additionally, that they are sensitive to the idiosyncrasies of the data. The parameters are estimated using 4 Markov chains with 2,000 iterations each, the first 1,000 of which serve as warm-up. The model converged based on visual inspection of the trace plots and theR statistics of 1.
'Hurricane season' is a binary indicator taking value 1 for the months June-November and 0 otherwise; 'Dry season' is a binary indicator taking a value 1 for the months December-March and 0 otherwise.
LOO − IC, leave-one-out information criterion; RM SE, Root mean squared error; N , number of data points.
σ 2 adjustment is for AR(1) autocorrelation as described in text.
We make the following observations. First, the average predicted value for September 2017 shows a huge difference, suggesting that there were many hurricane-related deaths. There To construct the main excess deaths estimates the predicted distribution is subtracted from the reported mortality number for September and October 2017.  To test the strength of the estimates a placebo test is conducted on the pre-hurricane months (January through August) in 2017. The essential idea is that there was no major event (treatment) during these months so we would expect that the death tallies for January-August 2017 to fluctuate within the predicted ranges during this period. If this assumption is violated -i.e. if we observe some large predictions errors during the prehurricane period -, then we should discount our belief in the claim that the post-hurricane prediction errors signify excess deaths attributable to the hurricane. In other words, if huge random fluctuation in the monthly death tallies are routine then it would seem plausible that the post-hurricane prediction errors are just part of this routine of substantial fluctuations rather than being caused by Hurricane Mara. However, the results show that the 2017 monthly tallies for February through August all fall well within the 95% UI of the predicted posterior distribution as displayed by figure 4b. One exception to this is the January tally which fails this placebo test by exceeding the upper bound of its 95% UI by about 70 deaths. This is of course just one test out of eight and not a resounding failure; it could be that the January 2017 value is a random occurrence. Perhaps something did occur during January 2017 which nudged the death tally for that month beyond the range of predicted deaths based on historical mortality patterns. Pending further investigation this estimate should be interpreted as a small warning sign that the excess death estimate could potentially be on the high side.
Focusing on a larger period, covering September through December 2017, leads to an excess deaths estimate of 920 (95%UI: [-40; 1870]). This central estimate is close to the September-October value but the bottom of the UI now dips below 0. As already noted, the November and December prediction errors look to be within the range of normal fluctuations. Applying the second type of placebo test shows that the average prediction error  At this stage the scale of the death toll attributable to Hurricane Mara during September and October of 2017 is fairly clear. There is a broad consensus among multiple studies based on data from the Puerto Rico Vital Statistics System that the excess death toll from the hurricane is probably somewhere around 1,000, with most studies, although not our own, putting the estimate somewhat above this figure. Even the initially resistant government of Puerto Rico has broadly recognized the reality discussed above with their updated figure of 1,427, which, surprisingly is on the high end of the independent vitalstatistics-based work. Unfortunately, the government has not presented details that would allow us to assess this figure although these details may be forthcoming when the GWU team releases further findings.
The only dissonant note comes from the survey-based Kishore et al. (2018) study which we think is best viewed as a quasi experiment that compares the survey approach for measuring excess, commonly used in war environments, to a vital-statistics-based approach in an environment where these statistics are highly reliable if not perfect. The outcome is sobering for the excess war-deaths literature: a central estimate that is far too high and an very wide uncertainty interval reaching beyond plus or minus 83 percent. Also discouraging are the t shirts, hats and placards displaying the number 4,645 at demonstrations in Puerto Rico which can easily be found on the internet.
We caution against projects attempting to make a list of all the particular individuals whose deaths, including non-violent ones, are attributable to the hurricane. Causes of death tend to be multiple and varied so while the hurricane will have played a role in many deaths it will often not be the single unique cause of death. In other words, the excess death concept is statistically, not individually, based.
We hope that our paper will have a lasting effect on methodologies used to estimate excess deaths following violent events in the future. We emphasize the utility of two main points from our methodology. First, there is the Bayesian regression approach that efficiently integrates all relevant data. Second, there are the placebo tests that help us distinguish between random fluctuations in death tallies and systematic movements in these tallies caused by a violent event such as a hurricane. Hopefully, the use of these techniques can reduce the scale of future controversies.