Plain Language Summary
In this study, we developed a statistical model to estimate the HRI emergency room visits by two regions (group of Counties) in North Carolina. We obtained daily temperature data assuming a minimal role of humans in greenhouse gas emissions, which is comparable to the pre-industrial revolution. Additionally, climate data was obtained for the next 100-years aggregated into 3 time periods, assuming two greenhouse gas emission scenarios. Using the statistical model that we trained using the actual observations during 2011-2016, we estimated the HRI risk attributable to climate change during 2011-2016 and for the future. Our results suggest that 15% of HRI emergencies during 2011-2016 were attributable to human-caused climate change. The rate of HRI emergencies was significantly higher, assuming a higher greenhouse gas emission scenario than the intermediate scenario, which was consistent across the future.
1 Introduction
Since the mid-20th century, the frequency and duration of hot days have increased globally due to anthropogenic climate change (Hoegh-Guldberg et al., 2018). Global temperature is very likely to increase up to 1.5-1.6°C during 2021-2040, 1.6-2.4°C during 2041-2060, and 1.4-4.4°C during 2081-2100, relative to1850-1900 (IPCC, 2021). These increases in temperature since the mid-20th century are primarily the result of anthropogenic greenhouse gas emissions (IPCC, 2021). Mora, et al., (2017) estimate that roughly 30 percent of the current global population is exposed to extreme heat conditions; this number is expected to increase to 50-75% by 2100 (Mora et al., 2017).
Exposure to extreme heat leads to heat-related morbidity and mortality. Extreme heat outcomes can be characterized as direct (e.g., heat-related illness (HRI)), or indirect (e.g., exacerbation of cardiovascular, respiratory, renal, endocrine, and mental health conditions) (Bell et al., 2018; Ebi, Capon, et al., 2021; Sarofim et al., 2016). Failure to acclimatize during extreme heat conditions can result in HRI ranging from muscle cramps, heat exhaustion, and heat stroke (Danzl, 2018). HRI, if untreated, can lead to life-threatening conditions (LoVecchio, 2016; Nemer & Juarez, 2019). In the United States (US), heat-related fatalities are more common than deaths due to other natural disasters (NOAA, 2019). In the US, there are an average of 702 heat-related deaths per year (Vaidyanathan et al., 2020). There is a potential imbalance between heat-related mortality and morbidity, posing an exponentially higher number of heat-related emergencies. For example, during a four-day heat wave event in North Carolina, there were 556 HRI emergencies compared to 1 heat-related death (NC-DHHS, 2016). The magnitude of the HRI emergencies compared to mortality demonstrate the need to focus on morbidity compared to mortality.
The role of climate change has been associated with the increasing trend of heat-related mortality and morbidity (Bell et al., 2018; Christidis et al., 2019; Ebi, Vanos, et al., 2021). Evidence-based climate detection and attribution play a key role in characterizing the changes in natural climate variability that are attributable to human activities (Ebi et al., 2020; Ebi et al., 2017). There is strong evidence supporting the association between future climate change and mortality (Conlon et al., 2016; Gosling et al., 2017; Guo et al., 2018; Lay et al., 2018). At the same time, attribution of human health risk to anthropogenic climate change is limited to considering mortality as a health outcome (Mitchell et al., 2016; Vicedo-Cabrera et al., 2021).
It is common for attribution analyses to investigate heat-related mortality, providing insight into the magnitude of extreme heat events on the most serious health outcomes. Based on the contrast between the frequency of heat-related mortality and morbidity in North Carolina, we hypothesized that using heat-related mortality would be an underestimate to quantify human health risks associated with climate change. We investigate heat-related morbidity to better understand the scope of human health burden associated with climate change. This study estimated the HRI attributable to anthropogenic climate change. Additionally, future HRI burden attributable to climate change was estimated using the climate projections driven by representative concentration pathways (RCPs). The RCPs are greenhouse gas concentration trajectories adopted by the IPCC that are used for climate modeling (IPCC, 2014). Future climate change is typically represented using four RCP scenarios, RCP2.6 is a stringent mitigation scenario, RCP4.5 and RCP6.0 are intermediate scenarios, and RCP8.5 is a scenario with very high greenhouse gas emissions (IPCC, 2014). This study includes two of the four RCP scenarios, comparing the intermediate emission scenario (RCP4.5) and higher emission scenario (RCP8.5).
The study objectives are 1) to estimate the HRI attributable to the current levels of anthropogenic climate change and 2) to estimate HRI associated with future climate change under two greenhouse gas emission scenarios (RCPs).
2 Materials and Methods
This study includes three analytic components: 1) Modeling and optimization of an epidemiologic model to estimate the rate of HRI emergency department visits, 2) Estimating the HRI burden attributable to current anthropogenic climate change, and 3) Quantifying the HRI burden associated with future climate change.
2.1 Study area
North Carolina has three physiographic regions with distinct climatological profiles: Coastal, Piedmont, and Mountain regions. The Coastal region includes 41 counties, Piedmont consists of 34 counties, and the Mountain contains 25 counties. Due to distinct weather conditions experienced by the population living in these physiographic regions, most of the heat-related research has been conducted using these sub-divisions (NC-DHHS, 2015). The study period includes summer months (May 01-September 30) over five years from 2011 to 2016. Data for 2013 was unavailable and is excluded from the analysis.
2.2 HRI morbidity data
The HRI data was obtained from the North Carolina Department of Health and Human Services (NC DHHS) that has partnered with 124 hospitals to collect statewide emergency department (ED) visit data to provide real-time, electronic public health surveillance, which is stored in the North Carolina Disease Tracking and Epidemiologic Collection Tool (NC DETECT, 2021). Heat-related illnesses were identified using ICD-9 CM codes with E992/E900.0/E900.0/E900; ICD-10 CM codes within T67/X32; and various keywords from the chief complaint and triage notes (Harduar Morano & Waller, 2017; NC-DETECT, 2021). We obtained daily aggregated counts of HRI ED visits. Days with fewer than five HRI cases were suppressed to maintain the confidentiality of patient identifiable information.
Decennial census population data at the county level from 2010 were aggregated to the three regions in North Carolina. Equation (1) was used to estimate the rate of HRI morbidity per study region.
\begin{equation} HRI\ morbidity\ rate=\ \frac{\text{Count\ HRI\ emergency\ department\ visits}}{\text{Population\ at\ risk}}*100,000\ \ \ \ \ \ldots\ (1)\nonumber \\ \end{equation}
2.3 Observed meteorological data:
Daily temperature data (mean (tmean), minimum (tmin), and maximum (tmax)) from the Global Historical Climatology Network-Daily (GHCN-D) database was extracted and aggregated by region. The GHCN-D dataset contains daily temperature measurements based on approximately five stations per region (Houston et al., 2012). Daily temperature measurements were homogenized to account for instrumentation and processing station observations (Rennie et al., 2019). Dew point data was obtained from the Parameter-elevation Regressions on Independent Slopes Model (PRISM) dataset (Model), 2019). The relative humidity (RH), maximum apparent temperature (MAT), US NWS/Steadman’s heat index (NWS_HI & Steadmans_HI), humidex, thermal discomfort index (TDI), and Excess Heat Factor (EHF) were computed by study region (Anderson et al., 2013; Baccini et al., 2008; Castelhano; & Laboclima, 2017; Langlois et al., 2013; PRISM, 2019).
2.4 Natural simulations:
The natural simulations are an estimate of daily maximum temperature (tmax-NS) in the absence of human-caused greenhouse gas emissions. The daily maximum temperature assuming non-anthropogenic climate change was obtained from the dataset developed by the Climate of the 20th Century Plus Detection and Attribution (C20C+D&A) Project (Stone et al., 2019). The C20C+D&A project is built on an ensemble of multiple dynamic models based on the atmosphere-land system. Due to lack of data, we excluded the year 2013 and 2016 in the analysis. We extracted daily maximum temperatures during the summer seasons of 2011, 2012, 2014, and 2015.
2.5 Climate projections
Localized and bias-corrected climate projections were obtained from the Localized Constructed Analogs (LOCA) database (Pierce et al., 2014). The LOCA dataset is statistically downscaled from the Climate Model Intercomparison Project 5 (CMIP5) and corrected for bias using constructed analogs (Pierce et al., 2014). The current study is based on study regions, amounting to more coarse geographies. The use of LOCA data with 1/16º resolution allowed us to assign localized temperature projections to finer geographies. We focused on the Community Climate System Model version-4.0 (CCSM4) and Geophysical Fluid Dynamics Laboratory (GFDL) model outputs as these models were outperformed compared to other climate models in the Southeastern United States (Zhang et al., 2013).
The climate projection dataset contains maximum temperature aggregated for three time periods: 1. Baseline (2011-2016), 2. Mid-century (2036-2065), and 3. Late century (2070-2099). The maximum temperature (tmax-FS) was estimated for each period assuming intermediate (RCP4.5) and high-emission (RCP8.5) scenarios.
2.6 Analysis
The analytic dataset contains the rate of HRI, tmean, tmin, tmax, RH, MAT, NWS_HI, Steadman_HI, humidex, TDI, and EHF at a daily scale (S 1). Additionally, we created a nominal variable to describe seven days of the week (DOW), a binary variable to identify weekend or weekday (Wday), as well as a nominal variable representing month and year.
2.6.1 Evaluating the relationship between heat metrics and HRI morbidity
Spearman correlations were run to determine temperature metrics to include in further analysis. Five of ten variables (tmax, tmean, NWS_HI, Steadmans_HI, and MAT) had a correlation coefficient greater than 55% (S 2). The exploratory analysis suggested a nonlinear relationship between the rate of HRI and heat metrics. The nonlinear relationship between the HRI rate and heat metrics was evaluated using the Generalized Additive Model (GAM) - ’mgcv’ package version 1.8-34 (Wood, 2021) and distributed lag nonlinear model (DLNM) approach (the ’dlnm’ R package version 2.4.5, Gasparrini et al. 2021). GAM is a semi-parametric framework to address the non-linearity using smoothing splines (Dominici & Peng, 2008). The model performance was compared using Akaike Information Criterion (AIC) and R-squared values. The statistical model using tmaxoutperformed to estimate the rate of HRI compared to other heat metrics (equation 2).
\begin{equation} \log\left(\mathbb{E}\left[\text{HRI\ morbidity}\right]\right)=\ \sum_{j=1}^{n}b_{j}\left(t_{\max}\right)\beta_{j}+\ \text{wDay\ }+\ Month\ +\ \text{Year\ }+\ \varepsilon\ \ldots(2)\nonumber \\ \end{equation}
The distributive effect of heat metrics on HRI was estimated using DLNM (Gasparrini, 2011). The DLNM follows an interrupted time-series approach, where daily HRI ED visits were assumed to follow the Poisson distribution and were fit using the GAM, controlling for seasonal effect. The association between heat metrics and HRI was evaluated for up to 5 lag days. Distributed lag nonlinear model (DLNM) approach - ’dlnm’ package version 2.4.5 (Gasparrini et al., 2021). and DLNM allows for an evaluation of the distributed effect of exposure variables on HRI morbidity.
2.6.2 Attributing the burden of HRI morbidity to current anthropogenic climate change:
The statistical model was trained by physiographic regions (Equation 2.) using tmax, time-series variables to estimate the rate of HRI. The model performance metrics were optimum while using three cubic regression splines for tmax for the Coastal region and four for Piedmont. The daily rate of HRI was estimated corresponding to the daily tmax-NS values. The percentage difference between the observed and estimated HRI rates assuming natural simulation was considered as the burden of HRI attributable to climate change. The mean difference between the daily rate of HRI between observed and natural simulation was tested using paired t-test. Additionally, the frequency of hot days between natural simulations and actual observations was compared using the chi-square test. The percent of HRI attributable to anthropogenic climate change is expressed as median percent per year and interquartile range (IQR).
2.5.3 Projecting HRI under future climate change scenarios:
Using tmax as a predictor, we trained a statistical model (Equation 3.) by physiographic regions to estimate the rate of HRI. Future HRI was estimated over three different 30-year periods (baseline, mid-century, and late century), focusing on RCP4.5 and RCP8.5 scenarios. The difference between HRI across two emission scenarios was evaluated using paired t-test and the differences between HRI across the three time periods were assessed using Analysis of Variance.
\begin{equation} \log\left(\mathbb{E}\left[\text{HRI\ morbidity}\right]\right)=\ \sum_{j=1}^{n}b_{j}\left(t_{\max}\right)\beta_{j}\ +\ \varepsilon\ \ldots(3)\nonumber \\ \end{equation}
3 Results
During the study period, 28.81% (219) of Coastal and 28.94% (220) of Piedmont regional observations were suppressed. The suppressed data were imputed with the median value (3) of the suppressed range. The Mountain region was excluded from the study due to poor data quality (50.13% (381) suppressed).
The mean HRI rate was 54.52 per 100,000 and 34.27 in the Coastal and Piedmont regions respectively. The annual HRI rate was consistently higher in the Coastal region than in Piedmont. In both study regions, the rate of HRI was higher (Coastal:40.94% higher; Piedmont: 28.47% higher) during the summer of 2015, compared to the study period (Table 1). The increase in the rate of HRI in 2015 could be due to a 14-day heat event with tmax exceeding 32°C, from June 13th-June 27th, 2015.
3.1 Association between daily maximum temperature and HRI:
The nonlinear association between HRI and tmax was established using the GAM (S 3-A). From the model diagnostics, we observed that about 80% of the deviance in HRI could be explained by equation 2 (S 3-B). Due to a higher number of observations with tmax between 26.7-35°C, there is a narrow residual confidence interval that reflects lower prediction uncertainty.
The results from the DLNM suggest that the risk of HRI significantly increased with tmax of more than 35°C. The HRI risk was higher during the day of exposure than the following days (S 4). When the tmax was recorded as 35°C, the HRI risk declined from about 2 to 1.2 from the day of exposure to the following day, potentially indicating a harvesting – or displacement – effect. It showed non-significant results during lag 2-5 days in the Coastal and Piedmont regions (S 5). The results from the DLNM suggest a negligible distributive effect of daily maximum temperature exposure on HRI morbidity.
As the HRI risk was higher during the day of exposure than in the latent period, further analysis is based on the same-day exposure-outcome relationship. As the primary goal of the current work was to build a prediction model rather effect estimation, the GAM (Equations 2, 3) was used to estimate the rate of HRI.
3.2 Burden of HRI attributable to anthropogenic climate change:
Over the four years studied (2011, 2012, 2014, and 2015), the frequency of hot days was about 30% higher in both the Coastal and Piedmont regions, among the actual observations than natural simulations (p- value < 0.001) (Figure 1). We observed a significant reduction in the mean rate of HRI morbidity in the Coastal (0.08 per 100,000; p- value < 0.001) and Piedmont (0.05 per 100,000; p- value < 0.001) assuming natural scenario than actual observations. In the Coastal region, 13.40% (IQR: -34.90,95.52) of the HRI morbidity is attributable to anthropogenic climate change and 16.39% (IQR: -35.18,148.26) in the Piedmont region. Based on our attribution analysis, about 83 HRI ED visits per summer season (152 days) in the Coastal region and 85 in Piedmont could be attributed to anthropogenic climate change.
3.3 Burden of HRI morbidity in contexr of future climate change:
Aggregate tmax values for the baseline, mid-century, and late-century were estimated assuming intermediate emission (RCP4.5) and higher emission (RCP8.5) scenarios using the CCSM4 and GFDL-ESM2M model outputs. In both the Coastal and Piedmont regions, we observed a significant increase in the median HRI.
In the Coastal region, during the mid-century, we observed up to 31.45% increase in the median HRI assuming a higher emission scenario compared to intermediate (p- value < 0.001). In the late century, the median HRI increased up to 78.65%, assuming higher emission scenario, compared to the intermediate scenario (p- value < 0.001). Additionally, assuming the intermediate emission scenario, the median HRI increased up to 53.01% during the mid-century and up to 67.98% in late century, compared to the baseline period (p- value < 0.001). Similarly, assuming a higher emission scenario, the median HRI increased up to 68.77% during the mid-century and up to 116.31% in the late century, compared to the baseline (p- value < 0.001).
In the Piedmont region, during the mid-century, the median HRI increased up to 24.17% assuming higher emission, compared to intermediate (p- value < 0.001). In the late century, the median HRI increased up to 65.85% assuming higher emissions, compared to intermediate emission scenario (p- value < 0.001). Additionally, assuming intermediate emissions, the median HRI increased up to 55.89% during the mid-century and up to 75.59% in late century, compared to the baseline (p- value < 0.001). Assuming higher emission scenario, the median HRI increased up to 77.28% during mid-century and up to 110.35% in late century, compared to the base line (p- value < 0.001).
4 Discussion
We identified that anthropogenic climate change attributed to higher frequency of hot days in North Carolina over a five-year period from 2011-2016. In North Carolina, about 15% of observed heat-related illness (HRI) during this five-year period is attributable to anthropogenic climate change. Additionally, there is a significant increase in HRI associated with current and future climate change.
Our findings that explore HRI are consistent with existing literature that attributes heat-related mortality to a changing climate. We observed a significant decline in the frequency of hot days assuming natural simulations than the actual observations, similar to Mitchell et al. (2016) and Vicedo-Cabrera et al. (2021). The HRI rate was significantly higher during the actual observations than in natural simulations. About 15% of HRI in North Carolina is attributable to anthropogenic climate change during the study period. The HRI attributable to anthropogenic climate change could be translated to an average of up to 85 HRI emergency department visits per physiographic region per summer season. Similarly, Vicedo-Cabrera et al. (2021) reported an average of about 23 heat-related deaths per year in 6 major cities in North Carolina attributable to anthropogenic climate change (Vicedo-Cabrera et al., 2021). Our results suggest that heat-related morbidity is 3.69 times higher than the heat-related mortality rate reported by Vicedo-Cabrera et al. (2021). The higher number of HRI could be due to the difference in the total number of cities included in our study (all the area covered under Coastal and Piedmont regions) compared to Vicedo-Cabrera et al., 2021, which included 6 major cities in North Carolina. Additionally, the natural simulation data used in this study is from the C20C+D&A project (Stone et al., 2019), whereas their simulation runs were obtained from the Detection and Attribution Model Intercomparison Project (DAMIP).
The use of heat-related mortality as a metric to estimate the burden of anthropogenic climate change could be an underestimate. In North Carolina, the annual mean heat-related mortality rate is about 200 times less than the rate of heat-related ED visits (0.11 (1997-2001) vs. 22.2 (2007-2012) per 100,000 persons) (Mirabelli & Richardson, 2005; Sugg et al., 2016). This study estimated the health burden of anthropogenic climate change using heat-related ED visits to minimize underestimation.
Along with the current climate risk attribution, we estimated the HRI associated with the future climate change. This study discussed the burden of HRI associated with climate change by comparing the HRI rate assuming intermediate and high emission scenarios. We observed a significant increase (up to 78.65% in Coastal and 65.85% in Piedmont regions) in HRI assuming a higher emission scenario (RCP8.5) compared to the intermediate emission scenario (RCP4.5). Additionally, the median rate of HRI significantly increased during the mid (up to 68.77% in Coastal and 77.28% in Piedmont) and late century (up to 116.31% in Coastal and 110.35% in Piedmont) compared to the baseline during both the emission scenarios. Similar results were reported by Lay et al. (2018), who estimated an increase in HRI emergencies by 32% in 2050 and 79% in 2090, assuming RCP8.5 compared to the RCP4.5 scenario. Kingsley et al. (2016) reported a 20% increase in HRI assuming the RCP8.5 scenario and attributable to climate change. Our results (up to 31.45% increase in HRI during mid-century and 78.65% in the late century) are similar to the HRI changes reported by Lay et al. (2018) and Kingsley et al. (2016).
Lay et al. (2018) estimated the attributable cost of heat on morbidity by exploring employer-based health insurance claims database of people under the age of 65. Discussing HRI in terms of cost often provides compelling insights that would effectively advocate policy change but were associated with limitations. The health data being used by Lay et al. (2018) excluded the most vulnerable population groups, such as the unemployed and elderly. In comparison, our study did not restrict vulnerable population groups from North Carolina.
The evidence-based findings from our study discussing HRI attributable to climate change, play a key role in public health education and preparedness that are relevant to extreme temperature exposure. Translating our results into public health action by developing community scale risk mitigation plans could substantially minimize the HRI risk. Additionally, population vulnerabilities such as age, gender, comorbidities, household type, income, nature of the employment and daily activity, are known to interact or mediate with temperature exposure in exacerbating HRI risk (Ebi, Vanos, et al., 2021). Along with population vulnerabilities, community build characteristics could mediate the HRI risk associated with climate change. Certain phenomenon such as the urban health island and heat dome effect, were identified to be driving factors associated with extreme heat exposure disparities by geography (Henderson et al., 2022; Tuholske et al., 2021). These phenomena are typically driven by neighborhood characteristics such as land use and land cover. Further studies discussing human health risks attributable to the current and future climate change, by considering population vulnerabilities and neighborhood characteristics could address these gaps in the research.
In this study, the future projections of HRI were estimated using static population (2010 decennial census). The objective of this study is to estimate the percent change in HRI over time, rather than presenting an absolute count of future HRI ED visits. Few studies adjusted for future population growth (Lay et al., 2018; Martinez et al., 2016) to describe the results based on absolute counts to estimate the cost associated with hospitalizations and ED visits. Due to data limitations, the human health burden associated with future climate change in this study was discussed using the percent change in the rate of HRI. We did not calculate the excess number of morbidities associated with future climate change, which is essential for the cost estimation. Additionally, changes in population characteristics across North Carolina physiographic divisions could influence our study results.
5 Conclusions
This work adds strong evidence quantifying the human health risk associated with current and future climate change in the Southeastern United States. This study estimated about 15% of the HRI in North Carolina during the study period is associated with anthropogenic climate change. Additionally, a substantial increase in HRI assumed a high emission scenario compared to an intermediate emission scenario. Our findings suggest that anthropogenic climate change has a significant effect on human health. Our results indicate that there is a need for evidence-based public health interventions that are informed by attribution.