It's not the heat, it's the vulnerability: attribution of the 2016 spike in heat-associated deaths in Maricopa County, Arizona

Extreme heat poses public health risks in many cities across the globe, especially in regions where heat is persistent and severe (e.g. Hartz et al 2013, Harlan et al 2014). Quantifying the public health burden of heat in places that experience this hazard most intensely, as well as understanding the causal drivers of that burden, may be instructive as cities prepare for a warmer future. Here, we explore factors that may have contributed to an unprecedented number of heat-associated deaths that occurred in one of the hottest metropolitan areas of the United States in the summer of 2016—one of the warmest seasons on record.

Many studies point toward a future with higher average temperatures as well as more intense and frequent extreme heat events as a result of greenhouse gas forcing and urbanization (e.g. Georgescu et al 2014, Maloney et al 2014, Russo et al 2014). These changes have also been projected to increase heat-associated morbidity and mortality (e.g. Åström et al 2013, Kingsley et al 2016, Petkova et al 2016). However, it is unclear how reliably historical models of the temperature-health association at short time scales (days) produce estimates of variability at longer time scales (e.g. seasons, years). Available evidence suggests that long-term variability may be difficult to model (Guo et al 2012). There is high uncertainty regarding the contribution of climate change as a driver of the future public health burden of heat in comparison with other factors such as population growth, demographic transition, and adaptation strategies (Hondula et al 2015, Petkova et al 2016, Gosling et al 2017). The interplay between these factors has rarely been examined retrospectively, but the limited evidence available suggests that social factors play a substantial role in influencing long-term trends in heat-associated mortality (e.g. Bobb et al 2014, Ng et al 2016).

One of the planet's warmest urbanized regions is the southwestern United States. Understanding the impact of heat in this region, where conditions are already severe throughout a prolonged warm season, may provide insights into how other cities may be able to cope with projected warming (Centers for Disease Control and Prevention (CDC) Environmental Public Health Tracking Program 2018). Maricopa County, Arizona, which spans most of the rapidly growing Phoenix Metropolitan Area, is home to more than four million people (US Census Bureau 2017). In 2016, the month of June was the warmest to date in the US National Oceanographic and Atmospheric Administration Climate Division within which Maricopa County is located, the month of July was the second-warmest, and the six-month warm season spanning May through October was the sixth warmest (NCEI 2018). In the same year, the Maricopa County Department of Public Health (MCDPH) reported a record high number of heat-associated deaths (MCDPH 2017). The 155 confirmed heat-associated deaths in 2016 exceeded the previous record by 48 cases and was more than 82% above the standardized prior ten-year average.

The objective of our analysis is to quantify the extent to which the abnormally high heat-associated death count of 2016 was driven by meteorology using heat-associated death records in concert with weather observations. In doing so, we hope to shed more general insights into the factors that reduce or exacerbate the impacts of heat and how preparedness and response strategies can improve. This knowledge is especially important to build in light of projections of large increases in heat-associated mortality in many cities resulting from global and urban warming.

Maricopa County, Arizona, is located in the Sonoran Desert in the southwestern United States and has a hot, arid climate. Summer daily maximum temperatures routinely exceed 38 °C (100 °F) and an average year brings 18 days of daily maximum temperatures of 43 °C (110 °F) or higher at Phoenix Sky Harbor Airport (CDC Environmental Public Health Tracking Program 2018). Rapid urbanization over the past several decades has led to large increases in daily minimum temperatures, a distinctive signature of the urban heat island effect (Chow et al 2012, Georgescu et al 2013). The county accounts for the majority of the geographical extent and population of the Phoenix-Mesa-Chandler Metropolitan Statistical Area, the 12th largest by population in the United States (US Census Bureau 2017).

Several regional agencies have made significant efforts to assess and reduce the public health burden of extreme heat. Among them, MCDPH operates an advanced and distinctive heat mortality surveillance program. Health authorities began more closely monitoring and classifying heat-associated deaths in 2006 after a highly publicized heat wave in 2005 resulted in greater local awareness of heat-associated deaths (Yip et al 2008). Since 2006, MCDPH has tracked and investigated heat-associated deaths, collecting data on location, time, demographics, and other circumstantial evidence for cases in which heat is suspected to have been the immediate cause of death or a contributing factor. The health department obtains this information from two sources: (1) the Maricopa County Office of the Medical Examiner, which flags deaths that it suspects to be heat-associated, and provides circumstances of death; and (2) the Arizona Department of Health Services, Office of Vital Registration, which contains a database of confirmed heat-associated cases. MCDPH extracts confirmed heat-associated cases by querying the Vital Registration database using key phrases (e.g. 'heat exposure', 'exhaustion'). The program also looks for certain International Classification of Disease-10 (ICD-10) codes that might indicate a heat-associated death (e.g. X30 (Exposure to excessive natural heat), T67.X (Effects of heat and light), P81.0 (Environmental hyperthermia of newborn)). MCDPH's efforts have culminated in a record of heat-associated deaths that is more inclusive than queries based on ICD codes alone (e.g, Berko et al 2014) but more precise than the construct of excess all-cause deaths that is often used for temperature–mortality studies (e.g. Anderson and Bell 2011, Gasparrini et al 2015). More details of the MCDPH heat surveillance programs are available online (https://maricopa.gov/1858/Heat-Surveillance).

MCDPH heat surveillance data from 2006 to 2016 were used for this analysis. This time period encompassed all years for which MCDPH has implemented their enhanced surveillance system for heat-associated deaths. The time series of heat-associated deaths spans the full calendar year, although the vast majority of cases (∼90%) are reported in the months May–September. Daily heat-associated counts were standardized to account for longer-term changes in population size and composition (e.g. age structure) as well as overall population health status. We directly standardized the data by adjusting the daily heat-associated death counts based on the total deaths from all causes reported in the county each year. We surmised that total deaths would be proportional to population size and composition as well as overall population health. All-cause annual mortality data were sourced from the Arizona Department of Health Services Office of Vital Statistics (http://pub.azdhs.gov/health-stats/report/ahs/ahs2015/pdf/5e1.pdf). All data were standardized to 2015, the most recent year for which total mortality counts are available.

Daily meteorological data were obtained from the National Centers for Environmental Information for Phoenix Sky Harbor Airport, an official first-order weather station located near downtown Phoenix. We obtained daily maximum, mean, and minimum dry bulb temperature, and examined the daily distribution of temperatures in 2016 compared to the ten years prior, with an emphasis on the five-month warm season (May–September) that accounts for the majority of heat-associated deaths. We did not incorporate humidity into the exposure models given prior research demonstrating little difference between temperature–mortality models that include or exclude humidity in the study setting (e.g. Harlan et al 2014). We also incorporated information from an experimental National Weather Service (NWS) product called HeatRisk (https://wrh.noaa.gov/wrh/heatrisk/). HeatRisk provides a daily classification of observed or forecast meteorological conditions based on the relative severity of hot weather compared to climatology and is being used as a decision support tool in NWS Weather Forecast Offices across the American Southwest at present.

Our methodological approach was to hindcast the 2016 heat-associated death counts based on temperature–mortality associations derived from the period 2006–2015 and meteorological observations from 2016. Because there is no consensus regarding the optimal statistical model or exposure variable for associating weather and heat-associated deaths in this region (e.g. Hondula et al 2014, Petitti et al 2016) or others (e.g. Barnett et al 2010), we used a number of different modeling approaches and exposure variables to create hindcasts for 2016. To minimize the introduction of bias into the hindcasting approach, the functional form of all models was determined ahead of detailed review of observations or model output.

We created three different types of models: (1) simple temperature–mortality; (2) complex temperature–mortality; (3) HeatRisk. The models estimate the daily association between the exposure variable and the risk of heat-associated deaths based on the observed time series of heat-associated deaths for the entire population of Maricopa County and the observed daily weather. For the simple and complex temperature–mortality models, we tested maximum, mean, and minimum temperature as the exposure variable. The simple temperature–mortality models took the following form:

$$\begin{eqnarray*}&&\mathrm{Log}(\text{scaled}\_\text{heatdeaths})={s}({\rm{temp}},\,{k}=4),\end{eqnarray*}$$

where scaled_heatdeaths is the standardized daily heat-associated death count from the MCDPH surveillance program, s is a natural penalized thin plane smoothing spline with k knots, and temp is the exposure variable (maximum, mean, or minimum temperature). We did not control for seasonality in these models because we anticipated that any seasonality in heat-associated death counts would primarily arise from meteorology and not confounding factors (following Petitti et al 2016). The complex temperature–mortality models took the following form:

$$\begin{eqnarray*}&&\begin{array}{l}\mathrm{Log}(\text{scaled}\_\text{heatdeaths})={s}(\text{complex}\_\text{temp},\,{\rm{by}}\\ \,=\,{\rm{as}}.{\rm{factor}}({\rm{month}}),\,{k}=4,)\end{array}\end{eqnarray*}$$

where complex_temp variable was applied to take into account the lagged effect of temperature on human mortality (e.g. Gasparrini et al 2015). The complex_temp variable was defined as n-day moving average of maximum, mean, or minimum temperature of the current and previous days. Alternative variants of the moving average from two-day (lag 0–1) to five-day (lag 0–4) were used in the complex model, and the moving average providing the best fit (based on the generalized cross validation score) was used as the final model. The HeatRisk model took the following form:

$$\begin{eqnarray*}&&\mathrm{Log}(\text{scaled}\_\text{heatdeaths})={\rm{as}}.{\rm{factor}}({\rm{HeatRisk}}),\end{eqnarray*}$$

where HeatRisk variable represents the NWS's HeatRisk classification (five categories: No Elevated Risk, Low Risk, Moderate Risk, High Risk, Very High Risk).

We generated hindcasts of 2016 heat-associated mortality using the observed daily temperatures from 2016 as input to each model. Daily hindcasts of heat-associated deaths were generated and summed by year. Model performance was evaluated by examining the deviance explained in the daily historical (2006–2015) heat-associated death records. Subsequently we reviewed the extent to which the model captured annual variability in heat-associated deaths counts, with an emphasis on how the hindcast for 2016 compared to other years. Confidence intervals (C.I.) for annual sums were calculated as the square root of the annual sum of the squared radii of the daily C.I. All analysis was conducted using R version 3.4.2, accessed through RStudio version 1.1.383. Time series modeling relied on the mgcv package (Wood 2017).

Climatological assessment of 2016

Warm season daily mean temperatures in Maricopa County in 2016 followed a progression similar to the prior ten-year average, peaking in June (daily average mean temperature = 34.9 °C) and July (daily average mean temperature = 36.2 °C). The highest daily mean temperatures were observed during two prolonged periods in which the average daily mean temperature was above 37 °C, one occurring 19–28 June and another that spanned most of the latter half of July. The highest daily mean temperature of the summer was recorded on 20 June (39.4 °C). Despite the high frequency of very warm days in 2016, compared to other years in the period 2006–2015, there were relatively few days at the highest end of the temperature distribution (figure 1). For example, there were nearly 15 fewer days with daily mean temperatures at or above 32 °C than the average of the previous ten years. More importantly from a health perspective, there were only eight days with daily mean temperatures above 38 °C compared to an average of 10.6 over the previous decade, and zero days with daily mean temperatures of 40 °C or above. An even more pronounced absence of days at the high end of the temperature distribution was observed for minimum temperatures, whereas the pattern of maximum temperatures was more closely aligned with the previous ten years (not shown).

Zoom In Zoom Out Reset image size

Temperature–mortality models

The mean standardized annual heat-associated death count in Maricopa County over the period 2006–2015 was 85.6, with a standard deviation of 24.6. In 2016 the number of reported cases (155) was more than 85.5%, or 2.8 standard deviations, above the historical average (figure 2).

Zoom In Zoom Out Reset image size

Simple temperature–mortality models all showed a clear and strong relationship with historical daily heat-associated deaths. Figure 3 illustrates the modeled relationship between daily heat-associated deaths and daily mean temperature, which was a better predictor of the relationship over the period 2006–2015 than minimum or maximum temperature. Differences between the models were, however, very small. The daily mean temperature model with zero lag accounted for 50.7% of deviance in daily standardized heat-associated deaths; minimum temperature accounted for 49.3% and maximum temperature accounted for 46.9%.

Zoom In Zoom Out Reset image size

There was a much weaker association of the relationship between modeled and observed heat-associated deaths when the model predictions and observations were aggregated by year (figure 2). The modeled annual heat-associated death counts fell in a narrow range (74.1–95.6) compared to the standardized annual observations (43.9–121.6). Model uncertainty was low relative to the inter-annual variability in the observations; the 95% C.I. for the annual heat-associated deaths counts (not shown on figures) was never wider than ±1.0 deaths per year. The aggregated counts from the zero lag mean temperature model explained approximately 19.7% of the variance in the aggregated annual observations. Using this model, we would have anticipated 79.7 heat-associated deaths in 2016 (95% C.I. = 78.8, 80.6) based on the observed temperatures. Predictions from models based on minimum and maximum temperature were 77.6 (76.9, 78.4) and 84.1 (83.0, 85.1) deaths, respectively.

Regardless of the exposure variable chosen, the prediction for 2016 and the upper bound of its 95% C.I. using a simple temperature–mortality model falls below the historical average number of heat-associated deaths. In fact, based on mean temperature, the prediction would have ranked as the 3rd lowest over the 11 years for which model predictions are available. The model consistently under predicted mortality across the three months of the warm season with the most heat-associated deaths, with June, July, and August associated with under-predictions of 42%, 51%, and 50%, respectively. Similarly, daily errors in model predictions were inversely correlated with temperatures, with larger errors occurring on the hottest days. The mean absolute error for daily predictions on days with daily mean temperatures above the 95th percentile in 2016 was 1.47 deaths, whereas on days between the 75th and 95th percentile the mean absolute error was 0.97 deaths.

Complex temperature–mortality models provided a better fit to the historical data than the simple models. After testing multiple combinations of exposure variables, smoothing windows, and interaction effects with month, we identified the best complex model to include four-day moving average daily mean temperature (lag 0–3), with an interaction term for month as a categorical variable. This model explained 56.0% of the deviance in daily standardized heat-associated mortality. Replacing mean with minimum and maximum temperature, keeping other parameters consistent, only changed the model diagnostics slightly (54.7 and 55.1% deviance explained, respectively). We did not find strong evidence of an added 'heat wave' effect in our exploration of more complex statistical models, and thus did not investigate the potential impact of heat wave attributes including duration, timing, and intensity (see Anderson and Bell 2011).

Despite the statistical improvement of the complex temperature–mortality models in fitting the historical data, the hindcast for heat-associated mortality in 2016 did not change considerably. Using the complex model with mean temperature, we estimated 82.0 (95% C.I. = 80.0, 84.1) heat-associated deaths in 2016 based on the observed meteorology. Hindcasts based on minimum (80.6 (78.6, 82.6)) and maximum (82.0 (80.2, 83.9)) temperature were similar, and all three were below the historical average of 85.6 heat-associated deaths per year (standardized). The complex model had a similar pattern of seasonal errors as the simple model, with larger errors on hotter days and in hotter months.

The model based on the NWS HeatRisk product did not fit the historical observations as well as the simple and complex temperature–mortality models. This categorical variable explained 33.1% of the deviance in observed heat-associated deaths. Effect estimates were greater than one heat-associated death for the two highest HeatRisk categories (red and magenta), with 1.25 deaths predicted per red day and 1.53 deaths predicted per magenta day. The model prediction for 2016 using HeatRisk was 77.3 (75.7, 78.8) heat-associated deaths, which is below the observed historical average and is the lowest number of heat-associated deaths predicted by any model we examined in this study.

Using widely applied statistical modeling techniques, we were unable to find evidence to suggest that the unusually high number of heat-associated deaths observed in 2016 in Maricopa County, Arizona, was related to observed meteorological conditions. Most models, in fact, suggested that 2016 should have been a year with a normal or below-normal number of heat-associated deaths. If population health risks related to heat are reasonably well approximated by widely applied time series models of the temperature–mortality association (e.g. Barnett et al 2012, Harlan et al 2014, Gasparrini et al 2016), our results indicate that factors other than the weather were predominantly responsible for the surge in heat-associated deaths in 2016.

The study findings serve as a case study for one geographic location and one defined time period and should be interpreted within that context. Replication to other settings with similarly available health records is a necessity and the methodological approach should be highly transferable elsewhere. If these results are more widely applicable, they challenge some existing methods and paradigms in the scientific literature concerning temperature-related health impacts. For instance, many studies project significant increases in heat-associated mortality in the future with models that are largely based on temperature change alone (Huang et al 2011, Hondula et al 2015). Projection-oriented studies typically do not take into account key factors that may act as effect modifiers such as timing and duration (Barnett et al 2012, Gasparrini et al 2016) environmental variables including air quality (Burkart et al 2013, Vanos et al 2014, Willers et al 2016), social attributes (Burkart et al 2014, Kovach et al 2015, Urban et al 2016), health outcomes observed in previous seasons (Rocklöv et al 2009, Ha et al 2011), and mortality displacement (Hajat et al 2005, Yang et al 2012, Saha et al 2013, Qiao et al 2015). Our findings suggest that these factors may be critical to include in assessments of future risk and potential uncertainties associated with applying contemporary models are quite wide. Our results are similar in nature to those of Guo et al (2012), who found high year-to-year variability in the relative risk of heat-related mortality among the elderly in the United States, which also suggests an important role for factors beyond those typically considered in temperature–mortality studies. If successful adaptation to heat was occurring in this community, the historical model should have overestimated the number of heat-associated deaths observed in 2016. Instead, nearly all of the models we generated led to an underestimate. Understanding the extent to which societies will continue to be able to adapt to heat, especially given projected urban- and global-driven climate change, is a critical research need (Åström et al 2017, Gosling et al 2017).

The major finding of our study, that an exceptionally warm summer should not have seen an exceptionally high number heat-associated deaths, may seem counterintuitive. However, conventional summary statistics commonly used to determine monthly and seasonal records are not optimal indicators of the risk of heat-associated deaths. Although 2016 was an abnormally warm summer in Maricopa County based on conventional metrics, the distribution of daily temperatures was relatively favorable from a health perspective. The nonlinearity of the temperature–mortality functions implies that the number of days at the uppermost end of the temperature distribution is especially influential in the predicted annual total of heat-associated deaths; the hottest days are associated with the highest modeled heat-associated mortality. The relative absence of these days in the record for 2016 (figure 1) contributed to a lower total number of hindcast heat-associated deaths for that year. There is some evidence that the timing of heat within the warm season may also be important to consider in modeling heat-related mortality (e.g. Gasparrini et al 2016). In Maricopa County in 2016 the six days with the highest mean temperature were split evenly between the months of June and July, which is consistent with the long-term climatology associated with the regional monsoon pattern that suppresses temperatures in late summer. There were two multi-day periods with anomalously high daily mean temperatures in the month of June, and June was the month with the greatest positive difference in monthly mean temperatures versus the prior ten years. The absolutely high daily mean temperatures in the month of June contributed to a higher hindcast total for annual mortality, but there was no evidence from the historical mortality record to suggest that this relative early-season warmth should have been especially dangerous.

Given that meteorological differences between 2016 and other years in our study period do not appear to account for the increase in heat-associated deaths observed in 2016, we investigated other factors that may be responsible. Because detailed information about the circumstances and exposure histories of decedents are not systematically available, even from an advanced surveillance program like the one operated by MCDPH, it is very difficult to make a strong causal argument based on the health records alone. We offer a few potential explanations drawing from our local observations and expertize. Of the many potential risk factors for heat-associated deaths that are tracked by MCDPH, two stood out with disproportionately high representation in 2016 cases: homelessness and age.

One possible contributing factor was an increase in unsheltered people experiencing homelessness in the region (MAG 2017). Unsheltered people experiencing homelessness are believed to be more vulnerable to environmental heat because of higher exposure and higher physical susceptibility related to pre-existing health conditions including mental health disorders (e.g. Longo et al 2017, Schmeltz and Gamble 2017). According to point-in-time survey data collected by the regional transportation planning association, Maricopa County has seen an increase in the general population experiencing homelessness since 2014 and experienced a 25% increase in unsheltered people experiencing homelessness between 2015 and 2016 (MAG 2017). Given the high degree of vulnerability to heat-associated with homelessness, an increase in the number of unsheltered people experiencing homelessness could increase the number of heat-associated deaths overall. According to MCDPH statistics, the number of heat-associated deaths among individuals experiencing homelessness over the period 2011–2015 ranged from 8 to 13 deaths per year, accounting for 9.5%–16% of total cases. However, in 2016 MCDPH reported 55 heat-associated deaths among people experiencing homelessness, which represented 35.5% of all heat-associated deaths recorded. Further, all heat-associated deaths for those experiencing homelessness were discovered outside, strengthening plausibility that living in unsheltered circumstances would increase a homeless person's vulnerability to heat. Though the causal pathway linking increased numbers of unsheltered people experiencing homelessness and the increase in heat-associated deaths among this population is still uncertain, the statistics suggest an association between shelter status and number of homeless deaths.

An additional population that experienced a high number of heat-associated deaths in the era of the MCDPH surveillance program was the 50–64 age group. In 2016, this age group accounted for 41% of heat-associated deaths. In a typical year, the 50–64 age group accounts for the highest number of heat-associated deaths, and this remained the case in 2016. However, in 2016 the rate of heat-associated deaths was also highest for the 50–64 age group compared, whereas it is usually the 65–74 and 75+ age groups that have the highest rate. When stratified by gender, the differences were more apparent; the mortality rate for 50–64 year old males well exceeded that of both 65–74 and 75+ males in 2016. This is significant as men consistently make up the vast majority of the deaths in this age group (73%–100% between 2011 and 2016) and maintain mortality rates that are at least twice as high as females indicating a potential high-risk group. Although homeless 50–64 year olds made up a significant proportion of the total deaths for 2016, homelessness does not appear to fully account for the unusually high rates in this age group. Given these statistics, it is possible that current intervention strategies overlook this population which may include many who are still of working age and maintain outdoor jobs, putting them at high risk for heat illness or death.

Limitations

Conclusions and recommendations

The authors acknowledge the efforts of colleagues at the Maricopa County Department of Public Health and the National Weather Service Weather Forecast Office in Phoenix, AZ for collecting and archiving the data resources used in this study. HP, DH, and MR were partially supported by the CDC's Climate-Ready Cities and States Initiative, grant no. 1 NUE1EH001318-01-00. MR was also partially supported by the CDC's Environmental Public Health Tracking Program, grant no. 1 NUE1EH001339-01-00. AU was partially supported by the Programme for research and mobility support of starting researchers, Czech Academy of Sciences, project no. MSM100421604 and by the Czech Science Foundation, project no. 18-22125S. The authors also thank Paul Chakalian at Arizona State University and Ken Komatsu at the Arizona Department of Health Services for their useful feedback on drafts of this manuscript.