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.