Rafael Irizarry ** 2018/06/08
Late last Friday, the Puerto Rico Department of Health finally released monthly death count data for the time period following Hurricane Maria:
The news came three days after the publication of our paper describing a survey conducted to better understand what happened after the hurricane. We did not have access to these data and, after this announcement, we requested daily data. We received it on Tuesday. The data comes as a PDF file which we make public here. The code to extract the data is long and messy so we include it at the end of the post and only show the package loading part here.
1 

The wrangling produces the object official
:
head(official)
1 

Here I include an analysis showing that the newly released data confirms that the official death count of 64 is a gross underestimate and that, according to these data, the current count is closer to 2,000 than 64.
The raw data
Here is a plot of the daily data for 2017 and 2018:
1 

We clearly see the effects of the hurricane: there is a jump of between 2550 deaths per day for about 2 weeks. It also seems that the 34 following months show a higher rate than expected. But this is hard to discern due to the seasonal trend: more deaths in the winter. The plot also shows that deaths have not yet been tallied for the last couple of weeks. So we restrict the analysis and exclude data from past April 15th.
Adjusting for an aging population, seasonal effect, and depopulation
We start by estimating the population size. We know that Puerto Rico has been losing population for the last 10 years. This file includes population size estimates. We also know that many people left the island after the hurricane. Teralytics tried to estimate this using data collected from cell phones. We extracted their estimates (by hand) from this post. We then interpolated to generate an estimate of population for the days in our data:
1 

To account for differences in population sizes we compute rates.
1 

Because the population is getting older, we also need to adjust for a year effect. To be able to compare 2017 to previous years, we use only the days before September 20 to compute a yearly median rate. Because we have no post hurricane data for 2018, we also assume that 2018 has the same rate as 2017.
2015  8.176414 
2016  8.387966 
2017  8.749981 
2018  8.749981 
With these estimates in place, we can now estimate a seasonal trend using 2015 and 2016 data. We do this by fitting a periodic smooth trend using loess:
1 

Here is what the estimate population size and seasonal trend look like:
1 

Adjusted Death Rates
With all this in place, we can now adjust the raw count data and compute an observed minus expected death rate. We use loess, with a breakpoint at September 20, 2017, to estimate a smooth version of the excess death rate as a function of time.
1 

Here is what the excess death rate looks like:
1 

We can now clearly see that following the big spike right after the hurricane, the death rate continued to be higher than expected, by about 10% until past December.
This gives us the excess death rate. To compute an excess count we need to plug in a population size and a baseline rate. If we assume all of 2017 had the same population, the prehurricane yearly effect and seasonal trend as other years we get the following excess deaths by date:
1 

We can see that 64 deaths were reached during the first few days and 2,000 deaths are reached around February.
One could argue that these assumptions are not optimal because, for example, young people may have left at a higher rate past the hurricane making the expected rate a bit higher. I encourage others to change the assumptions and rerun the code to see how the final results change.
Code for wrangling PDF
Here is the code to extract the data from the PDF file:
1 
