Can we use a machine learning approach to predict the impact of heatwaves on emergency department attendance?

Perth is the capital city of Western Australia (WA) (figure 1). It has a Mediterranean climate with hot and dry summers and cool, wet winters because of the domination of subtropical high-pressure systems [33]. The study was conducted in the Perth metropolitan area and covered 10 years (2006–2015) for the warm months (November–April).

Zoom In Zoom Out Reset image size

2.2.1. Population data

2.2.2. EDA data

2.2.3. Heatwave data

Heatwave was measured as Excess Heat Factor (EHF). It was initially available in gridded data with 5 × 5 kilometre pixels in a Network Common Data Form format provided by the Australian Bureau of Meteorology [5]. The data was converted into daily raster layers then was extracted from the layers and aggregated to (SA3) [35] using ESRI ArcMap software (Version 10.5).

An EHF value was calculated by combining the effects of excess heat (the long-term temperature difference categorized by every locality's climatology of heat, EHI_sig) and heat stress (the short-term (preceding 30 days) temperature difference based on recent thermal acclimatization, EHI_accl). EHF formulated as

$$\begin{eqnarray*}&&{\bf{E}}{\bf{H}}{\bf{F}}={{\rm{EHI}}}_{{\rm{sig}}}\times \,{\rm{\max }}\left(1,{{\rm{EHI}}}_{{\rm{accl}}}\right)\end{eqnarray*}$$

where EHI_sig = (DAT₀ + DAT₁ + DAT₂)/3 -DAT 95_th percentile, and

$$\begin{eqnarray*}&&{{\rm{EHI}}}_{{\rm{accl}}}=\left({{\rm{DAT}}}_{0}+{{\rm{DAT}}}_{1}+{{\rm{DAT}}}_{2}\right)/3-\left({{\rm{DAT}}}_{-1}+\ldots +{{\rm{DAT}}}_{-30}\right)/30\end{eqnarray*}$$

Average daily temperature (DAT) is the average of the maximum temperature and the next day's minimum temperature [5]. Thus, EHF provided a relative measure of the load, intensity, spatial distribution, and duration of a heatwave event.

In this study, a day was deemed as a heatwave day if the EHF value was greater than zero. Then the possible lag effects of heatwave on health were assessed. A cut-off value of 80% was applied to define a heatwave day for Perth, which means that if 80% of the study area (SA3s) had an EHF value >0 on a day, the whole Perth metropolitan area would be counted as a heatwave day.

2.2.4. Air quality data

ML was implemented using RStudio (Version 1.1.463). Several predictive models described below were evaluated for their accuracy. The whole data set was randomly divided into a 70:30 split (70% for training models and 30% for testing cross validation). Further validation of EDA was conducted by training the optimised model using 2006 to 2014 (nine years) data to predict 2015 (one year) data. The predicted values were then compared with the observed 2015 data. Finally, the optimised model was applied to the data for young children.

To assess the effect of heatwaves and air quality on EDA, all potential risk factors were included in the models, including all-age group, socioeconomic status (SES) measured by the socioeconomic index for area (SEIFA) categories (1 = disadvantaged and most disadvantaged area, 2 = middle area, 3 = second least disadvantaged and the least disadvantaged area) [38], daily heatwaves events (heatwave/non-heatwave), air pollutants (low/middle/high levels), landscape fire events, week periods (weekend/weekday) and public holidays (yes/no).

The Baseline Model does not include independent variables (predictors or risk factors). It only includes dependent variable (EDAs) in the training data to predict the number of EDAs. The Multiple Linear Regression model assumes that EDA has a linear function of risk factors. The Decision Tree Models [39] make no assumptions about the relations between the dependent and independent variables. Since the internal cross-validation is done at the same time as growing a tree, the result has an error measurement that is used to find the optimal number of splits for decision nodes and terminal nodes. The Full Decision Tree includes all nodes and splits, while the Pruned Tree removes trees that are non-critical, and only keeps the optimal number of splits in the model to enhance efficiency and predictive accuracy by reducing overfitting.

A Random Forest (RF) Model is an ensemble-learning method for which a multitude of decision trees are constructed at the training time and outputting the mean prediction of the individual trees to explain the relationships between EDA and risk factors [40]. A maximum of four variables (mtry, number of variables randomly sampled at each split) were selected for each node of a decision tree based on the default rule of 1/3 variables set in the RandomForest package. A total of 500 decision trees (ntree) were modeled [40, 41]. The best model was determined based on the maximum model coefficient of determination (R²) and the lowest error values (the root mean square error (RMSE) and the mean absolute error (MAE)). Finally, the actual and predicted values were compared to obtain the R² value.

The percentage increase in mean square error (%IncMSE) is a relative measure for ranking the importance of all risk factors. This measure was derived based on how much the accuracy decreased when a variable was excluded from the model [41]. Variables with larger positive %IncMSE values would have a significant impact on predicted EDA. The higher the score, the greater the contribution of a variable to predicting the EDA.

The Geographical Random Forest (GRF) Model is an extension of a RF that can address the spatial heterogeneity of the model. A newly introduced package, SpatialML [42, 43] was used to develop the GRF models with the local weight value set to one and the global weight value set to zero. An adaptive kernel for the GRF models was used. There were 21 SA3s in the study area, so 20 bandwidths were selected to consider 20 neighbors' weights for each SA3 modeling.

The GRF model was applied to young children as an exploration. The sum of air pollutants and EHF values for the years 2006 to 2015 was computed to derive a single value for each SA3. For example, to derive a PM_2.5 value for one SA3, the sum of 10 years of PM_2.5 concentrations for this SA3 was computed. The prediction for EDA was performed on the validation datasets with 30% of randomly selected data.

Several models were compared using ML approaches to estimate optimised model with the testing dataset.

Table 3 presents the model goodness of fit results for the all-age group and the young children's group. Among all the models examined, the RF models outperformed other models with the lowest RMSE and MAE values.

Results from table 4 display the importance ranking of risk factors for the all-age (model 1) and young children (model 2) groups, respectively. In model 1, age was the most important risk factor while in model 2, SES was the most important risk factor for predicting EDA.

For air quality measures, the top risk factors contributing to EDA were PM_2.5, PM₁₀ and fires in both models. Heatwaves (EHF) also contributed positively to predicting EDA in both models.

Further validation was performed using the RF model with the 2006–2014 EDA data and predicted EDA for 2015. Construct validation was then conducted by comparing results between the predicted and actual 2015 EDA data (figure 2). Adjusted by the population, the goodness of fit measure R² for EDA counts was 0.953, indicating the model fitted the data extremely well.

Zoom In Zoom Out Reset image size

Further analysis was conducted on the geographical variation of the risk factors by the GRF model for the 0–4 age group. This GRF model fitted extremely well (R² = 0.975) with low errors (RMSE: 0.00013, MAE: 0.00011). Overall, SES (SEIFA) remained the most important risk factor with the highest ranking, followed by heatwaves (EHF), and air pollutants.

The impacts of SEIFA, EHF, and PM_2.5 for predicting EDA in different SA3 areas were presented in figure 3. Overall, these three risk factors were relatively more important in the Perth southern suburbs than in other areas for young children, indicating a spatial interaction between them.

Zoom In Zoom Out Reset image size

The results from the RF model indicated that demographic factors such as age and SES were the most important predictors of EDA. The air quality measures were also important predictors, including particulate matters (PM_2.5 and PM₁₀). Fires are an important source of short-term particulate matters air pollution [44–46]. Heatwaves positively contributed to the prediction model for EDA.

In the past decade, very few studies used ML approaches to predict environmental factors such as heatwaves [47, 48] or air quality [49, 50]. Up until now, only one study in China [48] used an RF model and 3 years of data to predict the effect of heatwaves on heatstroke (R² = 0.70). But it did not compare different ML algorithms and didn't include air quality as a risk factor. Previous heatwave-related studies mainly used traditional epidemiological methods for prediction such as generalized estimating equations [28], Poisson regression [51, 52], and logistic regression models [53]. When there are complexes of independent variables or the method needs to rely on distribution patterns (e.g., autocorrelation or multicollinearity), the methods are limited. Such autocorrelations often exist for daily data on weather and air pollutants. As the RF modeling does not require a specific data distribution pattern, it can be applied to any type of data.

Another advantage of using the RF model is that it can determine the importance ranking of all risk factors in the prediction model. The RF model was built on the outcomes from multiple decision tree models, which can learn the patterns from the training data to predict future events more accurately. In this study, 500 cross-validations were performed to achieve the best outcome. In traditional epidemiological studies, age was identified as one of the important demographic factors [13, 54]. Using RF models, the current study demonstrated that age was the most important risk factor for predicting EDA. Apart from the elderly, children, especially young children are vulnerable to the effects of heatwaves due to their immature organ systems, poor adaptive capacities, and ability to produce more metabolic heat as compared to adults [55, 56]. In addition, children have a reduced ability to lose heat through sweating, and this can cause some severe outcomes such as convulsions and disorientation [57].

The SES was identified in this study as the second most important risk factor for predicting EDA in the all-age group model and the most important risk factor in the stratified young child model. A few studies observed that SES was a risk factor for heatwave-related mortality [58–62]. For example, a study in Hong Kong [58] used a Poisson model and showed that SES was an important risk factor for estimating the effect of heatwaves on mortality. People living in low socioeconomic districts might influence by the inequalities in intracity urban ambient air ventilation caused by the affordability of air conditioning and building infrastructure conditions. A case crossover study conducted in Adelaide [63] used multivariate logistic regression models and found that living alone (OR = 2.41, 95% CI: 1.32–4.40), with a lower socioeconomic status (OR = 2.10, 95% CI: 1.09–4.04) and no private health insurance (OR = 1.82, 95% CI: 1.05–3.16) increased the risk of heat related hospitalizations. Very few studies assessed the impact of SES on EDA. Using a generalised additive model, one study conducted in Brisbane [64] found that heat-related EDA was increased in most disadvantaged areas compared with the least disadvantaged areas for the 65–74 year age group during heatwaves. The Chinese study [48] used the RF method also found that SES was important for heatstroke prediction.

In the global models, air quality measures ranked high among the risk factors as compared to heatwaves, indicating the importance or sensitivity of air quality to EDA. Among all air pollutants, particulate matter was the most important risk factor. Currently, nearly no studies have used ML approaches to assess the impact of air quality and heatwaves on EDA for vulnerable populations, especially for children. A few epidemiological studies observed the significant adverse effects of heatwaves and poor air quality on health by other methods such as cross-correlation, generalised linear model, and generalised estimating equation model [12, 28, 65, 66]. Our study demonstrated that air quality should be included when estimating the effect of heatwaves on health as air pollution is often worse during a heat-wave event, especially the effect of particulate matter [12, 67, 68]. Moreover, our study showed that PM_2.5 from landscape fires was an important risk factor for the increased air pollution and EDA. This is consistent with a study conducted in Sydney, which found that air pollution arising from fire led to increased EDA [46]. A recent study by our team on landscape files in Perth also demonstrated that exposure to PM_2.5 during files was associated with increased overall EDA, EDA for cardiovascular diseases, and acute respiratory infection. Our study also found that there were more wildfire events during the summer season, indicating that heatwave was a potential contributing factor to wildfires [44]. During fires, more particulates can be generated. This creates a smog effect, making the air uncomfortable for breathing, which in turn leads to an increase in the use of emergency health services due to an increase in respiratory distress [69].

It is important to consider if risk factors such as climate conditions, social and environmental factors have any geographic variations for targeted intervention [70]. However, including spatial components in the modeling for estimating heatwave impacts on EDA was rarely reported. Although RF is a highly flexible algorithm, it does not account for spatial heterogeneity. The GRF, an extension of the RF algorithm, was recently introduced by Georganos, et al [43] as a novel geographic weighted RF approach to assess local variations of outcomes. To our knowledge, this is the first study using GRF to assess the impact of heatwaves, air quality, and other risk factors on health and to predict geographic variations in the EDA service demand for children in WA.

The goodness of fit test revealed that the GRF model fitted the data relatively well (R² = 0.975). The SES and EHF were the two most important risk factors for predicting EDA demand for young children living in Perth's southern suburbs compared to other areas. While the SES remained the most important risk factor, heatwave was the second important risk factor. This is different from the global RF model and indicates that heatwave has strong spatial variations across the study area (figure 3). We also found that there were spatial interactions between EHF, SES, and PM_2.5 for young children. This means the afore-mentioned risk factors contributed more to EDA in the similar locations. Using a Poisson regression model, a study conducted in Shanghai reported that with a 10 μg m⁻³ increase in concentrations of PM_2.5 (56.3 μg m⁻³ average daily concentration), the risk of emergency and outpatient visits for cardiovascular disease increased by 0.50% during the warm seasons [71]. However, a study conducted in Sydney did not observe a significant positive relationship between maximum temperatures, PM_2.5 (9.34 μg m⁻³ average hourly median concentration) and illness (fever and respiratory illness) for children younger than 6 years when the ARIMA time series model was applied [72]. The different outcomes from these two studies may be related to different levels of air pollutants exposure and possibly the sensitivity of their analytic methods.

This study has several strengths. Firstly, EDA was used to assess the acute health effects of heatwaves, which was more sensitive than a mortality indicator. Secondly, a ML approach was used to identify the optimised models (i.e., RF and GRF) to predict EDA during anomalies in weather and air quality with satisfying results. Thirdly, the gridded method was used for heatwave measurement instead of averaging station-based temperature data or using single-station data estimation. It makes the estimation of heatwave intensities (EHF) more accurate than those traditional measures. Moreover, the ML approach and RF model applied in this study are capable of modeling nonlinearity and interactions among different risk factors, which improves the model's goodness of fit. It also provides useful information on the ranking of the importance of all risk factors. Finally, the GRF model can provide information on local variations of risk factors. Such information can be of great value in prioritizing and allocating resources to communities in at-risk areas more accurately.

There are some limitations of this study. Firstly, this study only explored a few models. Future studies can include more models for assessment. Secondly, the current GRF method is not able to allocate large amounts of time series data directly into the model to assess temporal and spatial modeling simultaneously. Improving such an advanced method will be of great value for future studies. There is no doubt that more studies are needed to confirm the accuracy and reproducibility of the models produced in this study.