A Spatial–Temporal Bayesian Model for a Case-Crossover Design with Application to Extreme Heat and Claims Data

Abstract

Epidemiological approaches for examining human health responses to environmental exposures in observational studies frequently address confounding by employing advanced matching techniques and statistical methods grounded in conditional likelihood. This study incorporates a recently developed Bayesian hierarchical spatiotemporal model within a conditional logistic regression framework to capture the heterogeneous effects of environmental exposures in a case-crossover (CCO) design. Spatial and temporal dependencies are modeled through random effects incorporating multivariate conditional autoregressive priors. Flexible frailty structures are introduced to explore strategies for managing temporal variables. Parameter estimation and inference are conducted using a Monte Carlo Markov chain method within a Bayesian framework. Model fit and optimal model selection are evaluated using the deviance information criterion. Simulations assess and compare model performance across various scenarios. Finally, the approach is illustrated with workers’ compensation claims data from New York and Florida to examine spatiotemporal heterogeneity in hospitalization rates related to heat prostration.

1. Introduction

Environmental epidemiological studies frequently use aggregate health data to evaluate the short-term health effects of environmental exposures. For instance, daily counts of deaths, emergency department (ED) visits, and hospitalizations across a city have been correlated with short-term exposure to extreme temperatures [1,2].

The use of aggregate health data is an example of an ecological design, where the actual exposure corresponds to the distribution of individual-level exposure within the at-risk population [3,4]. However, it is impractical for large population-based studies to measure individual environmental exposures directly. Therefore, exposure surrogates that summarize individual-level exposure (e.g., the mean) are used. This approach can lead to exposure misclassification, resulting in biases and incorrect characterization of the uncertainties associated with health effect estimates [5,6].

In extreme climate conditions, extreme heat events (EHEs) are often the primary focus of study. Over short-term periods, there can be significant variability in exposure within the population due to (1) spatial differences in EHEs within the study area and (2) individuals spending time in various environments with different outdoor and indoor conditions. To tackle the first challenge, spatial–temporal models have been developed to estimate pollutant concentrations with fine spatial resolution, ensuring complete spatial–temporal coverage [7,8]. These estimates are then used to calculate short-term exposure values, which are linked to aggregated health outcome data across spatial units of interest (e.g., counties). Studies that account for spatial variability in exposure have been found to yield higher health effect estimates compared to those using data from sparse, stationary monitors, especially for pollutants with high spatial variability (e.g., traffic-related pollutants) [9,10].

Conversely, there has been less focus on addressing both location-based exposure heterogeneity and varying critical time windows simultaneously in ecological studies [11]. Factors contributing to spatial variation can include characteristics of the exposure (e.g., spatial differences in measurement error) and/or attributes of the study population that may be either excluded from the analysis or included as imperfect proxies for key covariates or confounders (e.g., race, ethnicity, and education level) [11,12]. To address spatial variation, Warren and Langlois [13,14] employed geostatistical methods to introduce spatial correlation between lagged regression parameters for different discrete spatial regions, defining correlation based on the distance between region centroids. However, this approach may be inadequate when spatial regions are irregularly shaped or vary significantly in size [15]. Furthermore, using distance-based spatial correlation models can increase computational complexity as the number of unique locations grows, due to the reliance on a Gaussian process with a spatial covariance function and the need to invert large matrices during each iteration of the model-fitting algorithm (typically achieved using Markov chain Monte Carlo sampling techniques) [16].

The study of short-term exposure effects using aggregated data is not a new concept. Some studies allow for heterogeneous effects through stratification or covariates adjustment [17,18], while others directly address these effects using the case-only approach developed by Armstrong [19]. In research on extreme temperatures, stratification has traditionally focused on demographic characteristics such as age, race/ethnicity, and sex. The case-only approach has also been applied to factors like chronic conditions, socioeconomic status, and various census tract characteristics [20,21,22,23]. We propose an extension to the widely used case-crossover (CCO) study design, where heterogeneous exposure–response relationships are modeled using a newly developed Bayesian spatial–temporal model [24]. This approach allows the model to flexibly learn and accommodate potentially complex heterogeneous exposure–response relationships during the fitting process. This study aims to (1) enhance the conventional case-crossover (CCO) design by incorporating a Bayesian spatiotemporal model to capture heterogeneous exposure–response relationships, thus enabling the analysis of environmental exposure effects across varied locations and times; (2) deepen understanding of differential impacts of environmental exposures on health outcomes by accounting for spatial and temporal dependencies; and (3) flexibly adapt to complex, location-specific exposure–response dynamics, allowing for the detection of population-specific variability that traditional models may miss.

The structure of the paper is as follows: the Methods section introduces the notation and provides a detailed description of the spatiotemporal models employed, along with the Bayesian algorithm used for parameter estimation and inference. The Simulation Study section describes the comprehensive simulations conducted to evaluate our approach. The Data Application section presents results from the motivating example using workers’ compensation claims data from Florida and New York. Lastly, the Discussion section provides concluding remarks and suggestions for future research.

2. Methods

2.1. Study Design

Case-crossover (CCO) designs are commonly used to analyze health data in environmental epidemiology studies when only case data, such as hospital visits, are available [25]. In a CCO design, each case is matched with a set of controls within a designated referent window, forming a stratum [26]. A widely used strategy for selecting this window is the time-stratified approach, which matches each case to three or four other dates within the same calendar month, ensuring the same day of the week [27]. This approach assumes cases are independent and occur infrequently enough that an individual would not experience the event more than once within the referent window. The time-stratified design is both localizable and ignorable, providing unbiased estimates of regression coefficients when conditional logistic regression [28] is used, enabling the analysis of exposure impacts over time. Figure 1 shows the study design of our real data example.

Our claims data include the date and county information for each claim, using injury codes to identify hospitalizations due to heat prostration. While the New York workers’ compensation claims dataset includes individual-level identifiers, the Florida dataset lacks this information. During the study period, a single individual may have multiple claims, and each year has unique individual identifiers. Additional details about the datasets are provided in the Supplementary Materials. Our data modeling approach is as follows:

$$\begin{array}{cc}\hfill & Y_{ijk}|,\mathsf{\beta }\stackrel{\mathrm{ind}}{\sim }\mathrm{Bernoulli}\{p_{ijk}\},\hfill \\ \hfill & \mathit{log}\{{\displaystyle \frac{p_{ijk}}{1-p_{ijk}}}\}={\mathit{x}}_{ijk}^T\mathsf{\beta }.\hfill \end{array}$$

(1)

For the $k\mathrm{th}$ claim $(k=1,\dots ,K_{ij})$ from the $i\mathrm{th}$ county $(i=1,\dots ,I)$ at the $j\mathrm{th}$ time point $(j=1,\dots ,J)$, ${\mathit{x}}_{ijk}$ denotes a $p\times 1$ vector of covariates for regression, which may include time-dependent or time-varying factors of interest. The vector $\mathsf{\beta }$, also $p\times 1$, contains the regression coefficients associated with the covariates ${\mathit{x}}_{ijk}$.

2.2. Bayesian Spatial–Temporal Model

To capture spatial heterogeneity, we introduce a spatial random effect into the model, defined as

$$log\{{\displaystyle \frac{p_{ijk}}{1-p_{ijk}}}\}={\mathit{x}}_{ijk}^T\mathsf{\beta }+{\omega }_i.$$

(2)

Using the approach outlined in the spatial–temporal Bayesian accelerated failure time models [24], we assume that the random effect ${\omega }_i$ follows a normal distribution $N(0,{\tau }^{-2})$. A conditional autoregressive (CAR) distributionBesag [29] for ${\omega }_i$ is then proposed.

$${\omega }_i|{\mathsf{\omega }}_{-i}\sim N\left({\displaystyle \frac{{\sum }_{i^{\prime }\ne i}{m}_{i{i}^{\prime }}{\omega }_{i^{\prime }}}{{\sum }_{i^{\prime }\ne i}{\omega }_{i{i}^{\prime }}}},{\displaystyle \frac{1}{{\sum }_{i^{\prime }\ne i}{m}_{i{i}^{\prime }}}}{\tau }^{-2}\right),$$

(3)

where ${\mathsf{\omega }}_{-i}$ = $\{{\omega }_1,\dots ,{\omega }_{i-1},{\omega }_{i+1},\dots ,{\omega }_I\}$. It is important to note that $m_{i{i}^{\prime }}$ is defined as 1 if county i is adjacent to county $i^{\prime }$; otherwise, $m_{i{i}^{\prime }}=0$. The joint density function of $\mathsf{\omega }={({\omega }_1,\dots ,{\omega }_I)}^T$ can be expressed as follows:

$$Pr\left(\mathsf{\omega }\right)\propto {\tau }^{-n}exp\left(-{\displaystyle \frac{1}{2{\tau }^2}}{\mathsf{\omega }}^T({\mathit{D}}_{\mathsf{\omega }}-\mathit{C})\mathsf{\omega }\right),$$

(4)

where ${\mathit{D}}_{\omega }=\mathrm{diag}\left\{{\sum }_{i^{\prime }\ne 1}{m}_{1i^{\prime }},\dots ,{\sum }_{i^{\prime }\ne I}{m}_{I{i}^{\prime }}\right\}$ is a diagonal matrix, with the $i\mathrm{th}$ diagonal element representing the total number of counties adjacent to county i. Additionally, $\mathit{C}={\left(m_{i{i}^{\prime }}\right)}_{I\times I}$ denotes the adjacency matrix for all counties included in the study.

Let ${\otimes }_k$ denote the referent window that includes observation times j for claim k. Given the observed data $\mathcal{D}=\{(Y_{ijk},{\otimes }_k),i=1,\dots ,I,j=1,\dots ,J,k=1,\dots ,K_{ij}\}$ and spatial random effects $\mathsf{\omega }$, and operating under the case-crossover assumption that ${\sum }_{j\in \otimes }{Y}_{ijk}=1\forall k=1,\dots ,K_{ij},$ the conditional data likelihood is expressed as

$$\mathcal{L}(\mathcal{D}|\mathsf{\omega })=\prod _{i=1}^I\prod _{j=1}^J\prod _{k=1}^{K_{ij}}f(Y_{ijk}|{\omega }_i)=\prod _{i=1}^I\prod _{j=1}^J\prod _{k=1}^{K_{ij}}{\displaystyle \frac{exp({\mathit{x}}_{ijk}^T\mathsf{\beta }+{\omega }_i+{ϵ}_{ijk})}{{\sum }_{j\in \otimes }exp({\mathit{x}}_{ijk}^T\mathsf{\beta }+{\omega }_i+{ϵ}_{ijk})}},$$

(5)

where ${ϵ}_{ijk}$ represents the error term in the regression model, and we assume ${ϵ}_{ijk}\sim N(0,{\sigma }_{ϵ}^2)$.

Given that the workers’ compensation claims data encompass patients enrolled across different years (i.e., temporal cohorts), we extended the application of the multivariate conditional autoregressive (MCAR) prior in our model for further analysis [24,30,31]. Specifically, we utilized the modeling approach proposed by Wang [24], incorporating an additional random effect vector for the temporal cohorts in the $i\mathrm{th}$ county, denoted as ${\gamma }_i={({\gamma }_{i1},\dots ,{\gamma }_{iJ})}^T$. Consequently, the spatial–temporal model can be expressed as follows:

$$\mathit{log}\{{\displaystyle \frac{p_{ijk}}{1-p_{ijk}}}\}={\mathit{x}}_{ijk}^T\mathsf{\beta }+{\mathit{z}}_{ik}^T\mathsf{\xi }+{\mathsf{\eta }}_i^T{\mathsf{\gamma }}_i+{\omega }_i,$$

(6)

where ${\mathit{z}}_{\mathit{ik}}={(z_{i1k},\dots ,z_{iJk})}^T$ is a $J\times 1$ vector with $z_{ijk}$ serving as a dichotomous temporal cohort indicator (1 for yes, 0 for no) for the $k\mathrm{th}$ claim from the $i\mathrm{th}$ county; $\mathsf{\xi }={({\xi }_1,\dots ,{\xi }_J)}^T$ represents the fixed temporal effect; and ${\mathsf{\eta }}_{\mathbf{i}}={({\eta }_{i1},\dots ,{\eta }_{iJ})}^T$, where ${\eta }_{ij}$ is a binary indicator covariate associated with the $j\mathrm{th}$ year-specific random effect ${\gamma }_{ij}$.

Let ${\mathbf{\Delta }}_i={({\omega }_i,{\gamma }_i^T)}^T$, and we assume

$${\mathbf{\Delta }}_i\sim N(\mathbf{0},{\mathbf{\Lambda }}^{-1})$$

$${\mathbf{\Delta }}_i|{\mathbf{\Delta }}_{-i}\sim MVN\left({\displaystyle \frac{{\sum }_{i\ne i^{\prime }}{m}_{i{i}^{\prime }}{\mathbf{\Delta }}_{i^{\prime }}}{{\sum }_{i\ne i^{\prime }}{m}_{i{i}^{\prime }}}},{\displaystyle \frac{1}{{\sum }_{i\ne i^{\prime }}{m}_{i{i}^{\prime }}}}{\mathbf{\Lambda }}^{-1}\right),$$

where $\mathbf{\Lambda }$ represents the hyper-parameters in the MCAR prior. Therefore, the joint density function of $\mathbf{\Delta }=\mathrm{vec}\left\{{\mathbf{\Delta }}_1,\dots ,{\mathbf{\Delta }}_I\right\}$ is

$$Pr(\mathbf{\Delta })\propto {|\mathbf{\Lambda }|}^{{\displaystyle \frac{1}{2}}}exp\left(-{\displaystyle \frac{1}{2}}{\mathbf{\Delta }}^T\left[({\mathit{D}}_{\omega }-\mathit{C})\otimes \mathbf{\Lambda }\right]\mathbf{\Delta }\right).$$

Similarly, the conditional likelihood of the observed data $\mathcal{D}=\{Y_{ijk},i=1,\dots ,I,j=1,\dots ,J,k=1,\dots ,K_{ij}\}$ is

$$\mathcal{L}(\mathcal{D}|\mathbf{\Delta })=\prod _{i=1}^I\prod _{j=1}^J\prod _{k=1}^{K_{ij}}\left[f\left(p_{ijk}\right)\right],$$

(7)

where $f\left(p_{ijk}\right)$ is the conditional density function, for the $k\mathrm{th}$ claim from the $i\mathrm{th}$ county in the $j\mathrm{th}$ temporal cohort. Denoting ${\mu }^{\prime }\left({\mathit{x}}_{ijk}\right)={\mathit{x}}_{ijk}^T\mathsf{\beta }+{\mathit{z}}_{ik}^T\mathsf{\xi }+{\mathsf{\eta }}_i^T{\mathsf{\gamma }}_i+{\omega }_i$, we have

$$f\left(p_{ijk}\right)={\displaystyle \frac{exp\left({\mu }^{\prime }\left({\mathit{x}}_{ijk}\right)\right)}{{\sum }_{j\in \Omega }exp\left({\mu }^{\prime }\left({\mathit{x}}_{ijk}\right)\right).}}$$

(8)

We consider the non-informative priors for $\mathsf{\beta }$, $\mathsf{\xi }$, and $\mu$, where

$$Pr\left(\mathsf{\beta }\right)\propto 1,Pr\left(\mathsf{\xi }\right)\propto 1,\mu \propto 1.$$

Also, with regard to random effects ${\mathbf{\Delta }}_i$, we select a conjugated prior for $\mathbf{\Lambda }$ with $\mathbf{\Lambda }\sim \mathrm{Wishart}\left(p,\mathit{R}\right)$. To ensure the prior remains vague, p can represent the dimension of $\mathbf{\Lambda },$ and $\mathit{R}$ can be arbitrarily specified as a diagonal matrix $Diag{\{100,\dots ,100\}}_{I\times I}$. Further details on the Bayesian inference using the MCAR prior can be found in the Appendix A.

2.3. Model Selection Criteria

We utilize the deviance information criterion (DIC) to select the best candidate model. Let $\mathsf{\theta }$ represent the complete parameter space of the model. Spiegelhalter et al. [32] introduced the following DIC within the Bayesian framework, which integrates the likelihood and posterior distributions:

$$DIC=2p_D+D\left(\overline{\mathsf{\theta }}\right),$$

(9)

where $\overline{\mathsf{\theta }}$ is the posterior mean of $\mathsf{\theta }$. It is important to note that $D\left(\overline{\mathsf{\theta }}\right)=-2\mathit{log}\mathcal{L}\left(\overline{\mathsf{\theta }}\right)+C$ represents the deviance of the model based on our posterior estimates, with C being a constant in the DIC. Additionally, $p_D=\overline{D}-D\left(\overline{\mathsf{\theta }}\right)$ is the difference between the mean of the deviance under the posterior distribution and the deviance evaluated at the posterior mean.

Inspired by the analysis of the workers’ compensation claims data, we considered the following two candidate models to simultaneously account for both temporal and spatial heterogeneity:

• Model 1 (M1): $\mathit{log}\{{\displaystyle \frac{p_{ijk}}{1-p_{ijk}}}\}={\mathit{x}}_{ijk}^T\mathsf{\beta }+\xi z_{ijk}+{\gamma }_i{z}_{ijk}+{\omega }_i$;
• Model 2 (M2): $\mathit{log}\{{\displaystyle \frac{p_{ijk}}{1-p_{ijk}}}\}={\mathit{x}}_{ijk}^T\mathsf{\beta }+{\mathit{z}}_{ik}^T\mathsf{\xi }+{\mathsf{\eta }}_i^T{\mathsf{\gamma }}_i+{\omega }_i$.

In Model 1 (M1), $z_{ijk}$ represents the year in which the claims were filed, treating the year as a continuous variable. In Model 2 (M2), ${\mathit{z}}_{ik}$ consists of indicator variables for the year when the claims were made. M1 incorporates a random slope to account for the temporal effect, defined as ${\mathbf{\Delta }}_i=\{{\omega }_i,{\gamma }_i\}$. In contrast, M2 includes a spatial–temporal intercept for each cohort and county, expressed as ${\mathbf{\Delta }}_i=\{{\omega }_i,{\gamma }_{ij},i=1,\dots ,I;j=1,\dots ,J\}$. We selected the optimal model for analyzing the workers’ compensation claims data based on the DIC, with a lower DIC indicating a better model fit. The procedures for parameter estimation, inference, and model diagnostics were implemented in R and utilized C functions. The relevant code is available at the GitHub link provided in Appendix B.

3. Simulation Study

3.1. Simulation Set-Ups

The comprehensive simulation study for this modeling approach is detailed in our previous work [24]. In this section, we concentrate on models that account for both temporal and spatial effects. We replicated the structure of the workers’ compensation claims data and performed extensive simulations across different scenarios. The data were generated using the following models:

• Scenario 1 (S1): $\mathit{log}\{{\displaystyle \frac{p_{ijk}}{1-p_{ijk}}}\}=x_{1,ijk}+x_{2,ijk}+{\omega }_i+{\gamma }_{ij}$;
• Scenario 2 (S2): $\mathit{log}\{{\displaystyle \frac{p_{ijk}}{1-p_{ijk}}}\}=x_{1,ijk}+x_{2,ijk}+{\omega }_i+{\mathit{z}}_{ik}^T\mathsf{\xi }$, where $\mathsf{\xi }=(0,0.5,-0.5,0.6,-0.8)$.

In this context, $x_{1,ijk}$ is a continuous variable generated from a standard normal distribution, $N(0,1)$; $x_{2,ijk}$ is a binary variable that follows a Bernoulli distribution, $Bernoulli\left(0.5\right)$. Both $x_{1,ijk}$ and $x_{2,ijk}$ vary across different claims, cohorts, and counties. The vector ${\mathit{z}}_{ik}=(z_{i1k},\dots ,z_{iJk})$ includes $z_{ijk}$ as the indicator for the $j\mathrm{th}$ temporal cohort. In Scenario S1, ${\omega }_{ij}={\omega }_i+{\gamma }_{ij}$ is generated iteratively, accounting for the dependency between adjacent cohorts, with ${\omega }_{ij}=2{\omega }_{i,j-1}+\zeta$, where $\zeta$ is drawn from a standard normal distribution, $N(0,1)$.

To simulate environmental exposure and produce comparable data, we also treated $x_{1,ijk}$ as a county-level risk factor in each scenario to assess our methods. In other words, $x_{1,ijk}$ varies only between counties and not among individual claims, meaning that $x_{1,ijk}=x_{1,i}$.

Given that Florida has 67 counties, we mirrored this in our setup, with $i=1,\dots ,67$. To streamline computational complexity, we assumed five temporal cohorts per county $(J=5)$. The number of claims, $K_{ij}$, for each county and temporal cohort was determined based on the proportion of cases in the actual workers’ compensation claims data, ensuring that each county had a minimum of five claims per day. To summarize the results, we generated 1000 Monte Carlo datasets for each scenario, and for each simulated dataset, we fitted the following models (M1 and M2) using a Bayesian framework:

• M1: $\mathit{log}\{{\displaystyle \frac{p_{ijk}}{1-p_{ijk}}}\}={\beta }_1{x}_{1,ijk}+{\beta }_2{x}_{2,ijk}+\xi z_{ijk}+{\gamma }_i{z}_{ijk}+{\omega }_i$;
• M2: $\mathit{log}\{{\displaystyle \frac{p_{ijk}}{1-p_{ijk}}}\}={\beta }_1{x}_{1,ijk}+{\beta }_2{x}_{2,ijk}+{\mathit{z}}_{ik}^T\mathsf{\xi }+{\mathsf{\eta }}_i^T{\mathsf{\gamma }}_i+{\omega }_i$.

It is important to note that ${\mathsf{\eta }}_i$ and ${\mathsf{\gamma }}_i$ are defined in Section 2.3, and the matrices ${\mathit{D}}_{\omega }$ and $\mathit{C}$ (the adjacency matrix) used in the MCAR and CAR priors were generated specifically for the counties in Florida. For parameter estimation, the posterior mean of each parameter was derived from 2000 MCMC samples, resulting in a total of 2500 posterior samples after discarding the first 500 samples as burn-in. While additional MCMC samples may be required depending on model convergence diagnostics, our numerical studies indicated that satisfactory results were achieved under these conditions (see results below). In summary, for the Monte Carlo replications, we report the bias (Bias), standard deviation (SD) of parameter estimates, mean squared error (MSE), and model selection probabilities based on the DIC.

3.2. Simulation Results

Here, we present the summary statistics for the estimates of the primary parameters ${\widehat{\beta }}_1$ and ${\widehat{\beta }}_2$ across various scenarios. Detailed results for all scenarios are available in

Table 1. Summary of the estimation results for the scenarios. Par: parameters; Bias: the bias of Monte Carlo average of parameter estimates; SD: Monte Carlo standard deviation of parameter estimates; MSE: mean squared error.

			S1			S2
	Par	Bias	SD	MSE	Bias	SD	MSE
M1	${\beta }_1=1$	−0.010	0.012	0.001	−0.004	0.012	0.000
	${\beta }_2=1$	−0.005	0.011	0.001	−0.001	0.022	0.000
M2	${\beta }_1=1$	−0.054	0.092	0.012	−0.003	0.012	0.000
	${\beta }_2=1$	−0.023	0.032	0.003	−0.002	0.022	0.001

. The bias values for ${\widehat{\beta }}_1$ and ${\widehat{\beta }}_2$ under M1 were the lowest among all candidate models; however, M1 displayed the poorest performance with the highest bias, especially in Scenario S1. In terms of mean squared error (MSE), M1 continued to show the best performance, yielding the smallest values across the different scenarios. Importantly, M2’s performance was similar to that of M1, demonstrating minimal bias and variability in Scenario S2.

We also assessed model fitting for scenarios where $x_1$ was the only risk factor that varied across counties; the results are presented in

Table 2. Summary of the estimation results for the scenarios with the county-level risk factor $x_{1,i}$. Par: parameters; Bias: the bias of Monte Carlo average of parameter estimates; SD: Monte Carlo standard deviation of parameter estimates; MSE: mean squared error.

			S1			S2
	Par	Bias	SD	MSE	Bias	SD	MSE
M1	${\beta }_1=1$	−0.084	0.131	0.024	−0.020	0.080	0.007
	${\beta }_2=1$	−0.006	0.030	0.001	−0.002	0.023	0.001
M2	${\beta }_1=1$	−0.092	0.148	0.030	−0.009	0.101	0.010
	${\beta }_2=1$	−0.045	0.063	0.006	−0.002	0.023	0.001

. Compared to the findings in Table 1, the bias of ${\widehat{\beta }}_1$ in M1 was significantly higher, particularly in Scenario S2, when compared to M2, despite the mean squared errors (MSEs) being similar. Regarding model selection across different scenarios, the selection probabilities based on the DIC criterion are outlined in

Table 3. Summary results of the selection probabilities for different models under different scenarios. $x_{1,ijk}$ is a subject-level risk factor varied across counties and temporal years; $x_{1,i}$ is a county-level risk factor.

	${\mathit{x}}_{1,\mathit{ijk}}$		${\mathit{x}}_{1,\mathit{i}}$
	M1	M2	M1	M2
S1	0.93	0.07	0.91	0.09
S2	0.10	0.90	0.22	0.78

. While Model M1 was generally preferred, in Scenario S2—characterized by spatial and temporal heterogeneity with a consistent linear temporal trend within counties—Model M2 recorded the highest selection rate. Additional simulation results for a spatial model that does not account for the temporal cohort effect variation (M3) can be found in Appendix C.

4. Data Application

4.1. Data Source

The source of injury and illness outcome data used in this research are the state of Florida and New York workers’ compensation claims datasets provided by the state workers’ compensation boards of Florida (https://myfloridacfo.com/division/wc/ accessed on 1 July 2023) and New York (https://www.wcb.ny.gov/ accessed on 1 July 2023).These claims datasets represent all types of illnesses and injuries, and include work-related injury and illness claims registered and accepted by the respective state’s workers’ compensation board from 1 January 2000 to 31 December 2020. The information that can be distributed to the public by workers’ compensation boards is subject to state statutes, with an aim of keeping claims information unidentifiable to a specific person. Workers’ compensation claims data represent perhaps the only source of data on worker illness and injury for a large state population. For each claim record, information includes the date of injury or illness, the county where the injury or illness occurred, the nature of the injury or illness, the industrial classification, and other characteristics dependent on state statutes such as age and gender for New York claims.

The nature of injury codes for each claim was used to identify certain heat-related illnesses and non-heat-related illnesses. While Florida uses the National Council on Compensation Insurance (NCCI) coding scheme, New York has used the Workers’ Compensation Insurance Organization (WCIO) coding scheme since 2013, and formerly used the Occupational Injury and Illness Classification System (OIICS) coding system from 2010 to 2015.

In Florida, there were approximately 1.7 million claims with indemnity benefits processed from 1 January 2000 through 31 December 2020. In New York, there were approximately 4 million claims from events that occurred between 1 January 2000 through 31 December 2020, and these claims included those with indemnity benefits as well as medical only, non-compensated to date, and canceled claims. About 0.05% of the claims in the New York dataset and about 0.2% of the claims in the Florida dataset could be categorized as heat-related illness claims.

Daily meteorological data on maximum temperatures were obtained from the U.S. National Center for Environmental Information through the National Oceanic and Atmosphere Agency (National Climatic Data Center https://www.noaa.gov/education/resource-collections/data/climate accessed on 1 January 2023) for the baseline period, 1980 to 1999, and the study period, 2000 to 2019. The data included daily maximum temperature (Tmax) for each county and date combination. For each county, a calendar day-specific 95th percentile threshold was calculated from the baseline data (1980–1999) using a 31-day window centered on the particular calendar day of interest, i.e., day in which an illness occurred. For example, the 95th percentile Tmax threshold for May 15 was determined based on all Tmax values for that county spanning May 1–May 31 during the 1980–1999 period (20 years × 31 days/year = 620 daily Tmax values). The 95th percentile maximum temperature values were referred to as the Extreme Temperature Threshold 95th percentile (ETT95). This exposure metric captured the frequency of extreme events and is relevant in the context of climate change. Extreme temperatures may exist for only one day, but could also take place for multiple days as a heat wave, which may then have a different effect on illnesses and injuries. We also explored one-day lag exposure periods to determine if there is a lag that occurs between exposure and injury. Considering the computation time complexity of the Bayesian approach, we applied our proposed models to the New York and Florida workers’ compensation claims data for 2010–2020. We focus on the 3103 heat prostration illness claims reported in the 2010–2020 New York and Florida workers’ compensation claims data.

4.2. Results

We explored spatial–temporal models that account for both spatial and temporal heterogeneity for this data application. For each model, 1000 samples were drawn during the burn-in period, followed by an additional 2000 samples drawn from the posterior distribution.

The results for fixed-effect parameters are summarized in

Table 4. Result summary for the data application of the heat prostration hospitalization using workers’ compensation claims data under M1 and M2 models. EST: parameter estimate; CL: credible limit.

Covariate	New York		Florida
	M1	M2	M1	M2
	EST (95%CL)	EST (95%CL)	EST (95%CL)	EST (95%CL)
Age
$>=$65	REF	REF	REF	REF
<65	0.039 (0.035, 0.041)	0.040 (0.036, 0.043)
Gender
Male	REF	REF	REF	REF
Female	0.057 (0.009, 0.097)	−0.092 (−0.155, −0.013)
EHE
Yes	REF	REF	REF	REF
No	−0.011 (−0.015, −0.007)	−0.022 (−0.029, −0.012)	−0.012 (−0.017, −0.004)	−0.018 (−0.023, −0.011)
Year(Cont)	0.044 (0.039, 0.052)	-	0.075 (0.068, 0.086)	-
Year = 2010	-	REF	-	REF
Year = 2011	-	−0.441 (−1.663, −0.116)	-	−0.036 (−0.107, 0.078)
Year = 2012	-	0.134 (−0.067, 0.143)	-	−0.114 (−0.157, 0.022)
Year = 2013	-	0.192 (−0.097, 0.331)	-	−0.112 (−0.165, 0.037)
Year = 2014	-	−0.429 (−0.877, −0.121)	-	−0.103 (−0.165, 0.050)
Year = 2015	-	0.095 (−0.570, 0.781)	-	-0.053 (−0.159, 0.110)
Year = 2016	-	0.416 (−0.026, 0.741)	-	0.052 (−0.030, 0.191)
Year = 2017	-	1.526 (0.761, 1.891)	-	0.277 (0.193, 0.407)
Year = 2018	-	1.495 (0.754, 1.901)	-	0.267 (0.193, 0.403)
Year = 2019	-	2.214 (1.514, 2.996)	-	0.619 (0.536, 0.765)
Year = 2020	-	5.375 (2.499, 7.816)	-	1.271 (1.018, 1.553)
Model Diagnosis
DIC	61,202.31	69,965.20	34,235.64	29,876.01

. Based on the DIC criterion, M2 was the best candidate model. Accordingly, a significantly higher number of claims was observed for those who were male, aged 65 and above, and had exposure to EHEs. A higher chance of hospitalization was observed for claims on EHE days. For instance, the average hospitalization of cases who were not exposed to EHEs was $0.01=1-\exp (-0.011)$ (95% CI 0.01–0.07) times less than that of those who were. In addition, the county-level risk factor, precipitation, was found to have a significant effect on mortality risk only in the case of M3. This finding agrees with our simulation studies, where M3 showed superior performance in parameter estimation compared with the other alternatives in the presence of county-level risk factors. In particular, the average heat prostration hospitalization of men who lived in the higher-humidity area was $1.0\%(1-\exp (-0.011)$ (95% CI 0.78–1.5%) more than that of men in lower-humidity areas. We observed similar results while using one-day lag exposure, as shown in

Table 5. Result summary for the data application of the heat prostration hospitalization using workers’ compensation claims data under M1 and M2 models for lag one EHEs. EST: parameter estimate; CL: credible limit.

Covariate	New York		Florida
	M1	M2	M1	M2
	EST (95%CL)	EST (95%CL)	EST (95%CL)	EST (95%CL)
Age
$>=$65	REF	REF	REF	REF
<65	0.038 (0.034, 0.040)	0.041 (0.035, 0.042)
Gender
Male	REF	REF	REF	REF
Female	0.056 (0.007, 0.098)	−0.090 (−0.154, −0.015)
EHE
Yes	REF	REF	REF	REF
No	−0.013 (−0.016, −0.008)	−0.023 (−0.029, −0.013)	−0.014 (−0.013, −0.004)	−0.018 (−0.023, −0.011)
Year(Cont)	0.046 (0.037, 0.052)	-	0.077 (0.066, 0.087)	-
Year = 2010	-	REF	-	REF
Year = 2011	-	−0.443 (−1.661, −0.116)	-	−0.036 (−0.107, 0.078)
Year = 2012	-	0.135 (−0.063, 0.144)	-	−0.114 (−0.157, 0.022)
Year = 2013	-	0.193 (−0.096, 0.331)	-	−0.113 (−0.165, 0.037)
Year = 2014	-	−0.431 (−0.876, −0.121)	-	−0.104 (−0.165, 0.050)
Year = 2015	-	0.096 (−0.570, 0.781)	-	−0.051 (−0.154, 0.110)
Year = 2016	-	0.416 (−0.026, 0.741)	-	0.052 (−0.030, 0.191)
Year = 2017	-	1.530 (0.777, 1.891)	-	0.277 (0.193, 0.407)
Year = 2018	-	1.496 (0.754, 1.901)	-	0.261 (0.193, 0.403)
Year = 2019	-	2.213 (1.514, 2.996)	-	0.620 (0.536, 0.765)
Year = 2020	-	5.374 (2.499, 7.816)	-	1.272 (1.018, 1.553)
Model Diagnosis
DIC	61,221.31	69,966.20	34,235.33	29,876.00

.

5. Conclusions

This study examined various models to determine the optimal approach for analyzing spatiotemporal heterogeneity in workers’ compensation claims data, specifically for heat prostration hospitalizations in New York and Florida. The findings highlight that Model M2, which incorporates both spatial and temporal random effects with time treated as a continuous variable, displayed superior performance in most scenarios, as indicated by low bias and variability in parameter estimates. This model’s accuracy in capturing the complex interactions between risk factors and exposure timing suggests practical value for informing heat-related public health interventions. Additionally, when strong hierarchical structures of risk factors across counties and temporal variability were present, Model M1 proved advantageous, offering flexibility in interpreting localized and temporally specific risks. This study’s use of the deviance information criterion (DIC) further underscored the robustness of these models in selecting appropriate frameworks for accurate, interpretable results.

Previous research on spatial and temporal modeling has largely focused on aggregate data analyses, often at the county or administrative level, or has incorporated either spatial or temporal correlations without adequately capturing both. This study advances current methodologies by integrating spatiotemporal dependencies within a Bayesian hierarchical framework, thereby addressing a critical gap for localized, population-level analyses of environmental exposures in health data. Unlike traditional methods, which may miss nuanced exposure–response relationships in ecological health studies, this approach leverages random effects for both spatial and temporal variability, aligning with emerging trends favoring multilevel Bayesian models in environmental epidemiology. By using continuous-time treatment, this study enhances upon previous static models, providing a nuanced understanding of variable exposure impacts across regions and time frames.

The study has several limitations. First, it does not address potential lag effects of extreme heat events (EHEs), which may impact health outcomes over delayed intervals. Additionally, the hierarchical structure of the data—although considered in the modeling—may still present challenges due to the potential misclassification of risk factors across population clusters. Another limitation is the computational intensity of Bayesian approaches, particularly with high-dimensional data and complex parameter structures, which can introduce longer processing times and potential convergence issues. Furthermore, the study’s reliance on workers’ compensation data, while valuable for heat prostration insights, may not capture broader population-level trends in heat-related hospitalizations. Although the exposure metric reflects the frequency of extreme events and is relevant in the context of climate change, we acknowledge the potential for misclassification of outcomes in claims data.

Future research could explore models that incorporate lag effects, such as distributed lag models (DLMs) or spatially varying Gaussian process models for critical windows (SpGPCWs), to better capture the delayed impacts of EHEs on health outcomes. Additionally, the adoption of advanced algorithms like blocked Gibbs sampling or slice sampling within the Bayesian framework could help reduce computational demands, making these models more accessible for extensive data applications. Research expanding this approach to other health conditions or datasets, like the Surveillance, Epidemiology, and End Results (SEER) program or the National Program of Cancer Registries, could provide broader applications for spatiotemporal modeling in health outcomes associated with environmental exposures. These developments could significantly improve risk prediction and public health planning in response to environmental stressors.