Extreme Heat and Emergency Health Impacts in the US (2018-2025)

1. Introduction

Extreme heat events stand out as one of the most deadly weather-related hazards worldwide. Heat-related events caused the most death in the U.S. during 2023, followed by wildfires and tornadoes. The largest number of reported injuries during that year were caused by heat, tornadoes and winter weather (National Oceanic and Atmospheric Administration (NOAA), 2024).

During the period 2022-2025, 36 countries faced national breaking temperature records (Yale Climate Connections, 2025). Five million deaths were attributed to extreme heat and cold globally between 2000-2019 (Chen et al., 2024), and projections indicate that heat-related deaths will increase.

Climate change might decrease cold-related deaths at high latitudes but increases heat-related deaths specially in Tropical Africa and Asia (Climate Impact Lab and United Nations Development Programme, 2024).

Impacts of these events can extend beyond individual increased mortality and morbidity to challenging the resilience and response of the health care system.

There is a wide spectrum of heat-related illnesses which range from dehydration and heat exhaustion, to cardiovascular and respiratory complications, with an uneven effect on vulnerable populations like elderly, young-children and individuals with pre-existing conditions. Prolonged heat exposure can affect the human body at multiple levels including mild physiological stress including dehy-dration, to life-threatening conditions as heat-stroke.

Thiel et al. (2025) identified a broad spectrum of heat-related diseases including electrolyte imbalances, cardiovascular disease and mental health issues. Using recent surveillance data, Vaidyanathan et al. (2024) report substantial increases in heat-related emergency department visits across the United States, with pronounced regional and demographic variability, underscoring the growing public health burden of extreme heat.

Focusing on vulnerable populations, Liss et al. (2017) demonstrate that higher ambient temperatures and heatwave events are strongly associated with increased hospitalizations among older adults, with clear seasonal synchronization between temperature peaks and morbidity. At a national scale, Vaidyanathan et al. (2019) estimate exposure–response relationships between heat index and hospitalizations across 1,617 U.S. counties, revealing significant geographic heterogeneity and showing that health risks can occur even at moderately elevated heat levels. Similarly, Gronlund et al. (2014) find that extreme heat and heatwaves are associated with increased hospital admissions among the elderly, particularly for renal and respiratory conditions, highlighting specific pathways through which heat affects morbidity. Recent work by Xing et al. (2026) provides detailed empirical evidence from 460 communities in Victoria, Australia, demonstrating that heat exposure contributes substantially to both hospital admissions and emergency department visits, with measurable economic and health burdens. Their analysis further shows that a significant proportion of these impacts can be attributed specifically to human-induced climate change, highlighting the role of anthropogenic warming in amplifying heat-related health risks.

Understanding the impacts of extreme heat on health and health services demand would provide a better estimation of future climate-related illness burden. In this paper we focus on evaluating the synchronicity between ambient variables peak timing and emergency department visits peak timing in the large contiguous US states by Health and Human Service (HHS) regions building predictive models for Emergency Department Visits (EDVs) that account for regional impacts of environmental factors as well as geographical factors. In section 2 we provide a description of all data sources used in the analysis. In section 3 we describe the methodology including the different modeling approaches for seasonality quantification and predictive models for EDVs across the 10 HHS regions. Results are provided in section 4 including model diagnostics and predictive skills. Discussion and Conclusions are presented in section 5 with further discussions on the limitations of this approach and future research directions.

2. Data Description

Data for this research was extracted from different sources. Date pre-processing was necessary for consistency and compatibility of different temporal and spatial scales. This section describes the different data sources and data processing for environmental and vulnerability related variables, including ambient maximum temperature and heat index, percentage of forest area, percentage of elderly and non-white population, as well as the health related outcome of this research which is Emergency Department Visits.

2.1. Emergency Department Visits

The emergency department visit (EDV) data was obtained from the Centers for Disease Control and Prevention’s (CDC’s) National Environmental Public Health Tracking Network (https://ephtracking.cdc.gov/DataExplorer/) under the heat-related illness (HRI) content. Data is available at a daily rate from 1/1/2018 to 12/31/2025 and for each of the ten Health and Human Services (HHS) regions (see Figure 1). EDV was expressed as HRI-associated visits per 100,000 emergency department visits.

Because EDV data was only available by HHS region, the temperature and heat index data were aggregated to this level. Furthermore, because the EDV data did not go further back than the 2018 (other than 12/31/2017, which was excluded for simplicity), temperature data was limited to this range despite being available at prior years.

2.2. Maximum Temperature Data

The temperature data was obtained from the Physical Sciences Laboratory of the National Oceanic and Atmospheric Administration (NOAA) (https://psl.noaa.gov/data/gridded/data.cpc.globaltemp.html). The maximum surface temperature in Celsius at a daily rate from 1/1/2018 to 12/31/2025 was used. Because the data was gridded in $0.5\times 0.5$ degree latitude and longitude squares, it needed to be aggregated to the HHS region level.

Shape-files from the United States Census Bureau (https://www.census.gov/geographies/mapping-files/time-series/geo/carto-boundary-file.html) were used at the 500 km resolution level. The daily temperature for each region was calculated by averaging the gridded values, weighted by the proportion of the area in the region to account for grid overlapping between two regions.

2.3. Heat Index

Heat Index (HI) data were obtained from the CDC’s National Environmental Public Health Tracking Network, specifically from the “HRI” section and the “Historical Temperature and Heat Index” subsection. These data are available at a daily resolution at the county level for the period January 1, 2018 through December 31, 2022, and were additionally aggregated to the HHS regional level.

The heat index integrates both, the air surface temperature and the relative humidity, and it is frequently used to measure heat exposure (Anderson et al., 2013). For calculating the Heat Index, CDC uses a National Weather Service’s modified version of the Rothfusz regression1.

$$\begin{array}{cc}\hfill HI=& -42.379+2.04901523T+10.14333127RH-0.22475541T\times RH\hfill \\ \hfill & -6.83783\times {10}^{-3}{T}^2-5.481717\times {10}^{-2}R{H}^2\hfill \\ \hfill & +1.22874\times {10}^{-3}{T}^2\times RH+8.5282\times {10}^{-4}T\times R{H}^2\hfill \\ \hfill & -1.99\times {10}^{-6}{T}^2\times R{H}^2\hfill \end{array}$$

(1)

where $HI$ is heat index (°F); T is the air temperature (°F) and $RH$ is the Relative humidity (%). The raw data are collected daily at the county level across all 3,143 US counties and sourced from the CDC website.

To align with the 10 HHS regions, we combined the data in two steps:

• First, compute the average heat index for each state:

  $$\begin{array}{cc}\hfill {\mathrm{HI}}_{\mathrm{state}}=& {\displaystyle \frac{1}{n_c}}{\displaystyle \sum _{i=1}^{n_c}}{\mathrm{HI}}_{{\mathrm{county}}_i}\hfill \end{array}$$

  (2)

  where $n_c$ is the number of counties in the state.
• Then, to better reflect the population distribution, we calculate the population-weighted average for each HHS region:

  $$\begin{array}{cc}\hfill {\mathrm{HI}}_{\mathrm{region}}=& {\displaystyle \frac{{\sum }_{j\in \mathrm{region}}{\mathrm{HI}}_{{\mathrm{state}}_j}\times {\mathrm{Pop}}_j}{{\sum }_{j\in \mathrm{region}}{\mathrm{Pop}}_j}}\hfill \end{array}$$

  (3)

where $Po{p}_j$ is the total population for state j and $H{I}_{statej}$ is the heat index for state j calculated from equation 2.

2.4. Vulnerability-Related Variables

Vulnerability-related data were obtained from the CDC’s National Environmental Public Health Tracking Network under the “HRI” (heat-related illnesss) ection, “Vulnerability and Preparedness: Heat” subsection. These variables capture demographic and environmental characteristics that may influence population susceptibility to HRI. Three variables were included in the analysis:

• elderly_percentage: Percentage of the population aged 65 or older living alone. This variable is collected annually at the state level from 2009 to 2021.
• nonwhite_percentage: Percentage of the population identifying as a race other than White. This variable is measured every decade at the county level from 2000 to 2010.
• forest_coverage: Percentage of land area covered by forest. This variable is measured every four years at the state level from 2001 to 2021.

The spatial distribution of the three variables by state and aggregated to HHS regions is presented in Figure 2, Figure 3 and Figure 4 for selected years.

Elderly population is highly concentrated in the states of Iowa, Missouri, Nebraska and Kansas (HHS Region 7), while the highest proportion of forested areas concentrates in the North East states of Maine, New Hampshire, Vermont, Massachusets and Rhode Island (HHS Region 1). In the case of non-white population the highest concentration occurs in the South West states of California, Nevada and Arizona (HHS Region 9).

3. Methodology

A first step in the analysis was to quantify the seasonality of the main response variable (EDV) and the proposed predictors: maximum temperature and heat index. The aim was to compare the features of the seasonal cycles (peak time and amplitude) and their regional variability, as well as to assess the potential synchronicity among them.

The detection of a potential lag effect between EDV and temperature, and between EDV and heat index, was also assessed by calculating the sample cross-correlation function between the response and the predictors after removing the trend and seasonal components using a seasonal–trend decomposition approach.

Finally, predictive models for EDV were fitted by incorporating the overall observed seasonal effect of the response, along with the effects of the predictors, which included either heat index or temperature as well as vulnerability-related variables. Regional variability in the strength of this dependence was accounted for through the inclusion of additional random effects in the model. Sequential model comparison was conducted as additional predictors were added by evaluating model predictive accuracy using a training–testing data split, while model fit was assessed using the Bayesian Information Criterion (BIC).

3.1. Seasonality Quantification

Variables EDV, temperature and heat index are strongly seasonal across all HHS regions; therefore seasonality quantification using useful and data informed metrics becomes an important initial step to describe the temporal patterns of these variables. For EDV, temperature, and heat index, seasonality metrics were calculated following Naumova et al. (2022).

Let $Y_{tj}$ be the response variable (EDV) at time t for region j, and $X_{tj}$ be the temperature or heat index at time t and region j. A log transformation was applied to EDV since it is a count variable.

We initially model the response as a purely seasonal model with one harmonic:

$$\begin{array}{cc}\hfill Y_{tj}=& {\beta }_{0j}+{\beta }_{sj}sin\left(\varphi t\right)+{\beta }_{cj}cos\left(\varphi t\right)+{\epsilon }_{tj}\hfill \end{array}$$

(4)

where $\varphi ={\displaystyle \frac{2\pi }{365.25}}$; ${\beta }_{0j}$, ${\beta }_{sj}$ and ${\beta }_{cj}$ are the unknown parameters for the one-harmonic periodic regression; and ${\epsilon }_{jt}$ are random normal errors with zero mean and constant variance within region.

Model coefficients are estimated by region using maximum likelihood, giving us regression coefficient estimates, their standard errors and covariances for each region. From these estimates we calculate:

• Phase angle $\psi =arctan({\beta }_{sj}/{\beta }_{cj})$
• $\mathrm{Var}\left(\psi \right)={\displaystyle \frac{{\beta }_{cj}^2{\sigma }_{sj}^2+{\beta }_{sj}^2{\sigma }_{cj}^2-2{\sigma }_{{\beta }_{cj}{\beta }_{sj}}{\beta }_{sj}{\beta }_{cj}}{{({\beta }_{cj}^2+{\beta }_{sj}^2)}^2}}$
• Peak timing $P_T=\left\{\begin{array}{cc}\psi /\varphi \hfill & \mathrm{if}{\beta }_{sj},{\beta }_{cj}>0\hfill \\ (\psi +\pi )/\varphi \hfill & \mathrm{if}{\beta }_{cj}<0\hfill \\ (\psi +2\pi )/\varphi \hfill & \mathrm{otherwise}\hfill \end{array}\right.$
• $\mathrm{Var}\left(P_T\right)={\displaystyle \frac{\sqrt{\mathrm{Var}\left(\psi \right)}}{\varphi }}$
• Amplitude $\gamma =\sqrt{{\beta }_{cj}^2+{\beta }_{sj}^2}$
• $\mathrm{Var}\left(\gamma \right)={\displaystyle \frac{{\beta }_{cj}^2{\sigma }_{sj}^2+{\beta }_{sj}^2{\sigma }_{cj}^2+2{\sigma }_{{\beta }_{cj}{\beta }_{sj}}{\beta }_{sj}{\beta }_{cj}}{{\beta }_{cj}^2+{\beta }_{sj}^2}}$
• Intensity $\theta =e^{\gamma }$

These metrics are then summarized by regions and represented graphically to assess spatial variability within their geographical setting.

3.2. Emergency Department Visits (EDV) Modeling

Enhancing effects of extreme air temperature and heat index are expected to produce amplified responses in the corresponding health-related outcomes. Under these assumptions, EDV may be subject to an overall seasonal effect that can be modified by localized seasonal effects, as well as by the impacts of temperature or heat index, which may also vary across regions due to local conditions.

Two sets of models are proposed to represent EDV outcomes, using either temperature or heat index as potential predictors. In total, six generalized linear mixed-effects models (GLMMs) were considered. Let $X_{tj}$ denote the temperature or heat index at time t for region j, and let ${\mu }_{tj}=E\left(Y_{tj}\right)$ denote the expected value of EDV. The link function is denoted by $g(·)$.

In these models, parameters $b_s$ and $b_c$ are the corresponding fixed effects regression coefficients for the one-harmonic sine and cosine functions, while ${\beta }_{sj}$ and ${\beta }_{cj}$ are corresponding random effects with values changing across regions. ${\beta }_x$ is the fixed effect regression coefficient for the external climatic predictors (maximum temperature or heat index) while ${\beta }_{xj}$ is the corresponding random effect with values changing across regions. Parameters ${\beta }_{0j}$ represents a random intercept or baseline value which changes across regions. The random effects are assumed to follow mean-zero distributions with variance components to be estimated from the data.

Two different sets of models are proposed to represent EDV outcomes using temperature or heat index as potential predictors. Six generalized linear mixed effects models (GLMM) were considered. Let $X_{tj}$ be the temperature or heat index at time t for region j. The models are:

1:

$g\left({\mu }_{tj}\right)=b_{s}sin\left(\varphi t\right)+b_{c}cos\left(\varphi t\right)+{\beta }_{0j}$

2:

$g\left({\mu }_{tj}\right)=b_{s}sin\left(\varphi t\right)+b_{c}cos\left(\varphi t\right)+{\beta }_{sj}sin\left(\varphi t\right)+{\beta }_{cj}cos\left(\varphi t\right)+{\beta }_{0j}$

3:

$g\left({\mu }_{tj}\right)=b_{s}sin\left(\varphi t\right)+b_{c}cos\left(\varphi t\right)+{\beta }_{0j}+b_x{X}_{tj}$

4:

$g\left({\mu }_{tj}\right)=b_{s}sin\left(\varphi t\right)+b_{c}cos\left(\varphi t\right)+{\beta }_{0j}+b_x{X}_{tj}+{\beta }_{xj}{X}_{tj}$

5:

$g\left({\mu }_{tj}\right)=b_{s}sin\left(\varphi t\right)+b_{c}cos\left(\varphi t\right)+{\beta }_{sj}sin\left(\varphi t\right)+{\beta }_{cj}cos\left(\varphi t\right)+{\beta }_{0j}+b_x{X}_{tj}$

6:

$g\left({\mu }_{tj}\right)=b_{s}sin\left(\varphi t\right)+b_{c}cos\left(\varphi t\right)+{\beta }_{sj}sin\left(\varphi t\right)+{\beta }_{cj}cos\left(\varphi t\right)+{\beta }_{0j}+b_x{X}_{tj}+{\beta }_{xj}{X}_{tj}$

where $g(.)$ is the link function and ${\mu }_{tj}=E\left[ED{V}_{tj}\right]$ is the expected value of EDV for time t and region j.

Models 1 and 2 exclude the effect of an environmental predictor and are purely seasonal, while the remaining models 4 to 6 include the effect of an environmental variable (maximum temperature or heat index) on the EDV response. Models 2, 5 and 6 include both a fixed global seasonal effect and a random seasonal effect that varies across regions. Models 4 and 6 include both a fixed global effect of X over EDV, and a random effect of X over EDV that changes across regions. While model 1 considers a random intercept which modifies the baseline values of EDV from region to region, Model 2 also accounts for regional modifications of the seasonal effects over the global seasonal signal, by adding random effects on the sine and cosine components.

All models were compared in terms of predictability, by estimating the model parameters for a training set, and estimating the the mean squared error (MSE) on a testing set, and in terms of goodness of fit. Different train-test split data percentages were compared to evaluate the sensitivity of the prediction errors to the different split schemes.

4. Results

The seasonality analysis of the main time series in this research is presented in Section 4.1. As an initial step, the time series are treated as independent across the different HHS regions, and separate models are fitted for each region. Exploration about the delay between the environmental time series and the EDV response time series is presented in Section 4.2. Following this analysis, we present the EDV predictive models that include maximum temperature as the main external covariate in Section 4.3 (maximum temperature models), followed by the EDV predictive models that use the heat index as an external covariate in Section 4.4 (heat index models).

4.1. Seasonality Analysis

The one-harmonic seasonal models fitted to the observed time series of EDV, maximum temperature and heat index are shown in Figure 5, Figure 6 and Figure 7. The model (red lines) is able to reproduce the overall seasonal behavior for all time series at each HHS region, more specifically, the peak timing, but they are not able to reproduce the more extreme events (blue dots). Representation of regional differences in amplitude and inter-annual variability requires additional model components or covariates.

Seasonality metrics for EDV, temperature, and heat index are shown in Table A1, Table A2, and Table A3, respectively. These results show an average peak timing across regions of 194, 200.4 and 203.4 days from the beginning of the year for EDV, maximum temperature and heat index respectively. The EDV peak occurs earlier than the peak timings for both maximum temperature and heat index. There is a difference of nearly one week between the average peak timing of EDV and that of maximum temperature. This difference increases to nearly 10 days on average when comparing the average peak timing of EDV with that of the heat index.

Figure 8, Figure 9, Figure 10, and show maps of the estimated peak times for EDV, temperature, and heat index, respectively. The earliest peak timing for EDV occurs in HHS Region 5 (Illinois, Indiana, Michigan, Minnesota, Ohio, and Wisconsin), while the latest peak occurs in Region 9 (Arizona, California, and Nevada) among the continental states. This pattern is not necessarily synchronized with the earliest and latest peak times for maximum temperature and heat index.

An inspection of detailed regional differences between the response and the seasonal effects reveals that a maximum difference of two weeks occurs in Region 5 (Illinois, Indiana, Michigan, Minnesota, Ohio, Wisconsin) when comparing EDV and heat index peak timing, while a minimum difference of one-quarter of a day between peak times occurs in Region 6 (Arkansas, Louisiana, New Mexico, Oklahoma, Texas). Similar extreme differences in peak timings are observed between EDV and maximum temperature for the same regions. These regional differences are accounted for through random effects that represent the region-to-region variability of the seasonal cycle model representation.

4.2. Lag Detection Between EDV and Environmental Predictor Variables

A common way to quantify the time delay between two time series and their association strength is by calculating the sampling cross-correlation functions (CCF) and looking the lag values where the highest sample cross-correlation occurs (Runge et al., 2014). To minimize the effect of the autocorrelation of each time series on this calculation, we applied the STL decomposition (Seasonal–Trend decomposition using LOESS) to remove seasonal and non-seasonal trends in the data (Cleveland et al., 1990).

STL decompositions for EDV, temperature, and heat index are shown in Figure A1, Figure A2, and Figure A3, respectively.

Cross-correlations between the STL remainders for EDV and temperature and EDV and heat index are shown in Figure A4 and Figure A5, respectively.

In most cases, except for the CCF between EDV and maximum temperature residuals in Region 5, the maximum cross-correlation values occurred at lag = 0. This indicates that no significant time delay between the response and the climatic variable was detected, suggesting an immediate heat response. In view of these results, lagged predictors were not included in the proposed models described in the following section. Models that include maximum temperature as a predictor are referred to as temperature models and are discussed in Section 4.3, whereas models that include the heat index as a predictor are referred to as heat index models and are discussed in Section 4.4.

4.3. Temperature Models

Temperature models were evaluated using both the full data (to the end of 2025) and data limited to the heat index time range (until the end of 2022). To test the sensitivity of different train-test data splits, different percentages were considered for each data range. Figure 11 compares the MSEs for models 1 to 6, including temperature as a predictor in models 3 to 6.

As shown in Figure 11, models 3, 4, 5, and 6 had much lower MSEs than models 1 and 2. The lowest MSEs occurred at a 70-30 train-test split for both data ranges, with model 4 and model 6 being optimal for the data ending in 2025 and 2022, respectively. However, the difference in MSEs was very small, and all four models including temperature had a similar performance.

The best two models, models 4 and 6, were used to make predictions in the testing set. Figure 12 shows the predictions for these two models fo all HHS regions. Both models suggest distinct regional differences in the associations between EDV and maximum temperature.

4.4. Heat Index Models

Heat index models were evaluated using data that ends at the end of 2022, as shown in Figure 13. Similar to the temperature models, models 3, 4, 5, and 6 had much lower MSEs than models 1 and 2 (which did not include heat index as a predictor). Though the four models with heat index show similar predictive performance at a 70-30 train-test split, model 5 had the lowest MSE.

Figure 14 shows the predictions for the model with the best predictive performance (model 5). Results are shown for all HHS regions.

In contrast with the best predictive models selected in the temperature case, model 5 suggests a more homogeneous dependence of EDV on the environmental predictor (heat index) across regions, excluding a random effect for the dependence between the mean EDV and the heat index. An estimated regression coefficient ${\widehat{b}}_x=0.1893$ implies an overall 20% increase in the mean EDVs rate per unit increase in the heat index.

Both temperature and heat index models are able to adequately capture the seasonal patterns of EDV. Extreme observations generally lie within the 95% prediction intervals but remain more difficult to reproduce, possibly due to the presence of overdispersion, which is not accounted for by the Poisson model.

Overdispersion can be accounted for by using the Negative Binomial 2 family as the probability distribution for the response. Model comparisons between the Poisson and Negative Binomial 2 (NB2) families were conducted. The results are presented in the following section.

4.5. Fitting Performance vs. Prediction Accuracy for Temperature and Heat Index Models

The Bayesian Information Criterion (BIC) is a model selection tool to evaluate how well a model captures the overall distribution of the observed EDV counts when it used within the context of model comparison. Lower values of BIC indicate better goodness of fit after penalizing for model complexity. In this case, the data exhibit clear overdispersion, meaning that the variance exceeds the mean. The Poisson model imposes equality between the mean and variance, which leads to an inadequate representation of the variability in the data, particularly in the tails. In contrast, the NB2 model relaxes this restriction by allowing the variance to exceed the mean, enabling it to better capture the observed dispersion. As a result, the NB2 model provides a substantially higher likelihood by fitting the full distribution more accurately. This improvement in likelihood outweighs the modest complexity penalty in BIC, yielding lower (better) BIC values as shown in Figure 15 and Figure 16 for temperature and heat index models respectively.

Mean Squared Error (MSE), on the other hand, evaluates how accurately a model predicts the conditional mean of EDV and does not incorporate information about the variance structure. From this perspective, the simpler Poisson model can have an advantage because it directly estimates the mean without introducing an additional dispersion parameter. This can result in more stable point predictions and slightly smaller prediction errors. In contrast, the NB2 model’s added flexibility allows for greater variability around the mean, which may increase the variability of point predictions and lead to a marginally higher MSE. Consequently, the Poisson model can achieve better predictive performance in terms of MSE, despite providing an inferior fit to the overall distribution. Figure 17 and Figure 18 show a comparison of the MSEs between the Poisson family and Negative Binomial 2 family for the temperature and heat index models respectively.

In this work, we focus on improving model predictability by treating the Poisson model as a preferable alternative to the NB2 model.

4.6. Intra-Class Correlation Coefficient

The intra-class correlation coefficient (ICC), also referred to as the variance partition coefficient, measures the proportion of total variability attributable to the random effects associated with the data clusters (i.e., region-to-region variability). Fixed effects influence the ICC indirectly by reducing the unexplained variance components used in its calculation.

To quantify the contribution of regional variability in emergency department visits (EDV), including seasonal components and covariate effects, we computed ICC values both before adjusting for model covariates (unadjusted ICC) and after including covariates such as seasonal effects and temperature or heat index (adjusted ICC).

Table 1 presents the ICC estimates for the best-performing temperature and heat index models under both the Poisson and Negative Binomial type 2 (NB2) families. For the Poisson models, unadjusted ICC values are relatively low (approximately 5–10%), indicating modest baseline clustering and that only a small fraction of the total variability is attributable to differences between regions prior to covariate adjustment.

In contrast, the adjusted ICC values for the Poisson models are very close to 1. This indicates that, after accounting for seasonal and meteorological covariates, nearly all of the remaining unexplained variability is attributable to between-region differences. In other words, the residual variation becomes highly structured at the regional level, with very little within-region variability remaining once covariates are included.

For the Negative Binomial models, adjusted ICC values remain substantially lower (approximately 0.10–0.14). This reflects the additional observation-level variability introduced by the NB2 dispersion parameter, which accounts for overdispersion beyond the Poisson assumption. The increased within-region heterogeneity inflates the total variance, thereby reducing the proportion attributable to between-region variation and leading to smaller ICC values even after covariate adjustment.

4.7. Residual Analysis

Residual analysis is a critical step in regression model diagnostics, used to detect model misspecification, including incorrect functional forms, distributional assumptions, or the presence of influential outliers. For generalized linear mixed models (GLMMs) with non-normal responses, a simulation-based residual approach is recommended over traditional Pearson and deviance residuals (Feng et al. (2020)). In this framework, residuals are defined as simulation-based quantile residuals taking values in the interval $[0,1]$. They represent the position of each observed value within the cumulative distribution of responses simulated from the fitted model, i.e., $r_i=F_{model}\left(y_i\right)$

This simulation-based quantile residual approach, implemented in the R package DHARMa (Hartig (2024)), was used for residual diagnostics. Figure 19 shows the DHARMa residual plots for the best temperature model (Model 5) fitted with a negative binomial (NB2) distribution. The left panel presents a quantile–quantile (QQ) plot against the expected uniform distribution, where each point represents a scaled residual and the red diagonal line indicates the theoretical expectation under correct model specification; systematic deviations from this line indicate potential model misspecification. The right panel displays residuals versus predicted values (rank-transformed), where the residuals should be uniformly distributed and centered around 0.5 with no discernible pattern under a well-specified model. A notable feature is the presence of several residuals equal to 1, indicating that the corresponding observed values lie at or beyond the upper tail of the simulated distribution (i.e., nearly all simulated responses are smaller than the observations), suggesting that the model may underestimate extreme values.

Analogous to Figure 19, Figure 20 presents the DHARMa residual plots for the best heat index model (Model 5) using the NB2 distribution family.

Analogous results are observed for the heat index model (Figure 20). The diagnostic tests yield consistent conclusions for both models. The Kolmogorov–Smirnov (KS) test for uniformity is significant, indicating that the residuals deviate from the expected uniform distribution. In contrast, the dispersion test is not significant, suggesting that both models adequately account for overdispersion. However, the outlier test is significant, indicating an excess number of observations falling outside the simulated range, which is consistent with a lack of fit in the upper tail of the distribution.

Taken together, these diagnostics suggest that, although the negative binomial models successfully address overdispersion, they fail to adequately capture extreme observations and possibly other structural features of the data. The significant deviation from uniformity further indicates residual misspecification, potentially arising from unmodeled nonlinear effects, temporal dependence, or missing covariates. While the NB2 specification provides an improvement over the Poisson model, additional model refinement is required to better represent the distributional characteristics of heat-related emergency department visit counts.

5. Discussion

This study examined the relationship between extreme heat indicators and emergency department visits (EDVs) for heat-related illnesses across the ten U.S. Health and Human Services regions over the period 2018–2025. By combining seasonal analysis with generalized linear mixed-effects models, we provide insights into the temporal dynamics, regional variability, and predictive performance of climate–health relationships at a national scale.

Our results confirm that both maximum temperature and heat index are important predictors of heat-related EDVs. Models incorporating these environmental variables consistently demonstrated lower prediction error compared to models based solely on seasonal components. Among the two indicators, models using the heat index achieved slightly better predictive performance and goodness of fit than those using temperature. This finding is consistent with the construction of the heat index, which integrates both temperature and relative humidity to better represent human heat exposure.

A key result of this analysis is the immediate heat effect on the response variable EDV. Cross-correlation analyses showed that the maximum correlation between EDVs and both temperature and heat index occurred at lag zero in most regions. This indicates that heat-related health impacts are largely contemporaneous with exposure, emphasizing the importance of real-time monitoring and rapid response to extreme heat conditions at this spatial and temporal scales. These results align with the results presented in Xing et al. (2026) where an immediate effect of heat response was observed and persisted for 21 days.

The seasonal structure of EDVs further highlights important dynamics in heat-related health outcomes. The average EDV peak occurs around day 194 of the year, while maximum temperature and heat index peak later, around days 200 and 203, respectively. The observed differences of approximately 6 days for temperature and 9 days for heat index, suggest that EDVs begin increasing before environmental variables reach their maximum values. This pattern aligns with previous findings reported in the literature, where increased heat-related morbidity precedes peak temperature conditions, likely reflecting early exposure effects and population vulnerability conditions (Vaidyanathan et al., 2019), (Liss et al., 2017).

Regional heterogeneity was also observed in the relationship between temperature and EDVs. Models incorporating random effects indicated that the strength of this association varies across regions, reflecting differences in environmental conditions and regional characteristics. In contrast, the best-performing models suggested that the relationship between heat index and EDVs is more consistent across regions, indicating that the heat index may capture more uniform physiological responses to heat exposure.

Model comparison revealed a trade-off between predictive accuracy and distributional fit. While Poisson models achieved lower mean squared error (MSE), indicating better predictive performance, negative binomial (NB2) models provided improved goodness of fit by accounting for overdispersion in the data. Residual diagnostics further indicated that, despite addressing overdispersion, both model classes exhibit some degree of misspecification, including deviations from uniformity and the presence of extreme observations outside expected modeled data ranges. These results suggest that additional model refinement, such as alternative distributional assumptions or more flexible model structures, may be necessary.

The intra-class correlation coefficient (ICC) analysis provides additional insight about the sources of EDVs’ variability. Unadjusted ICC values were relatively low (approximately 5–10% for the Poisson family and 4% for the NB2 family), indicating that region-to-region differences alone explain a limited portion of total sources of variability. However, after adjusting for seasonal and environmental covariates, most of the remaining unexplained variability is attributed to regional differences in the Poisson case. For the NB2 model case, additional variability is accounted for by overdispersion, and only a moderate fraction of the remaining unexplained variability (11%-14%) after fixed effects adjustment is attributed to regional differences.

First, the analysis is constrained by the availability of EDV data, which spans only from 2018 onward, limiting the evaluation of longer-term trends. Second, spatial aggregation at the HHS regional level may mask finer-scale variability and limit the detection of localized effects. Third, vulnerability-related variables were excluded due to lack of statistical significance in the models, which may reflect limitations in data resolution rather than absence of real effects. Considerations on different population strata as sex, age and socio-economic status might reveal additional differences in heat responses.

6. Conclusions

In summary, this study demonstrates that environmental heat indicators, particularly heat index, are valuable predictors of heat-related emergency department visits across the United States. Models incorporating temperature or heat index substantially outperform purely seasonal models, underscoring the importance of including environmental drivers in predictive frameworks.

The finding that EDV peaks occur earlier than the peaks of temperature and heat index has important public health implications. It suggests that interventions such as heat-health warning systems and preparedness strategies should be implemented prior to the highest levels of heat exposure, rather than in response to peak conditions.

From a methodological standpoint, mixed-effects modeling provides a useful framework for capturing both shared seasonal trends and regional variability. However, further work is needed to address model limitations related to overdispersion, residual structure, and extreme responses.

Future research should focus on improving spatial and temporal resolution of both health and environmental data, expanding the temporal coverage of EDV records, and enhancing access to integrated datasets. Improved data availability will be critical for better understanding lagged effects, vulnerability factors, and regional disparities in heat-related health outcomes.

Overall, this study contributes to the growing evidence linking extreme heat to increased healthcare demand and provides a basis for developing regionally adaptive public health strategies aimed at mitigating the impacts of a changing climate on population health.