Forecasting hospital discharges for respiratory conditions in Costa Rica using climate and pollution data

1.   Introduction

Respiratory diseases are conditions that affect the organs and tissues within the lungs and airway systems, leading to difficulty in breathing. These illnesses are currently the leading contributors to the global burden of diseases when assessed through disability-adjusted life-years. Moreover, the rising healthcare costs associated with these diseases are placing a growing strain on the economies of nations worldwide, particularly due to the "big five" conditions: Chronic obstructive pulmonary disease, asthma, acute respiratory infections, tuberculosis, and lung cancer ^[1].

To maintain a balanced health economy, it is crucial to understand the distribution and the incidence of these diseases to optimize healthcare systems, ensure efficient resource allocation, and enhance healthcare access and quality. Therefore, our focus is to provide a reliable forecast for respiratory hospitalizations at the hospital level in Costa Rica by using potential climatic factors as inputs in predictive models, enabling health authorities to anticipate future hospital needs. Specifically, three different competitive predictive models (Random Forest, XGBoost, and PROPHET) along with an ensemble method are used and compared.

Extensive research indicates that climate influences respiratory health, resulting in higher rates of hospitalization and mortality in varying seasons ^[2,3,4,5,6]. Furthermore, it is well known that human activities significantly impact the short-term and long-term climate by releasing greenhouse gases and other pollutants into the atmosphere, leading to more extreme weather events in recent years. As a consequence, researchers worldwide have been studying the effects of pollution and climatological factors on human health in different regions of the world ^{[7,8,9,10,11]}.

Specifically with respiratory diseases, climate factors are associated with extreme weather ^[5,6,12], and pollution, such as particular matter (PM) ^[4,13]. Airborne PM is a complex mixture of different chemical species originating from numerous sources. These particles can be a product of combustion, suspension of soil materials, suspension of substances from the sea, and can also be formed by chemical reactions in the atmosphere ^[14]. Furthermore, elevated concentrations of PM not only substantially heighten health risks but also contribute significantly to the challenges associated with climate change. Therefore, the monitoring of PM is essential for comprehending and mitigating air pollution, aiming to enhance public health. Nevertheless, the measurement of PM presents challenges in developing countries due to the constrained monitoring infrastructure, resulting in gaps in data analysis both temporally and spatially.

On the other hand, an aerosol is defined as a stable suspension of solid and liquid particles in gas. It is possible to obtain related measurements with the facilitation of remote sensing techniques. Aerosol optical depth (AOD), a satellite-derived metric that quantifies the presence of aerosols across the entire atmospheric column, is a good option. Several researchers around the world have found that the driving factors of AOD are related to urban and economic development, agricultural activities, industrialization, landscape aggregation, and regional transportation ^[15,16]. Moreover, AOD is associated with socio-economic factors, such as GDP, industry, and vehicle density ^[17].

Costa Rica, a country with an area of 51,179 $Km^2$, is located in Central America with an estimated population of approximately 5,003,402 inhabitants in 2018. It is a stable democratic country with a socioeconomic development model based on the opening and economic liberalization, upholding human rights, and ensuring universal basic access to essential goods and services, such as education and health system, while preserving extensive natural resources. ^[18].

Although Costa Rica is renowned for its biodiversity, with most environmental policies focusing on the protection and conservation of natural resources, unregulated land use has been observed, resulting in significant environmental impacts. This issue is not limited to the Great Metropolitan Area, the country's most urban and densely populated region, but is also prevalent in other parts of the national territory. Consequently, this form of human development can intensify vulnerability to natural disasters and extreme weather conditions^[18].

Regarding to its health system, since the 1940s Costa Rica prompted its national health system and social investment by creating the Costa Rican Social Security Fund (CCSS, from its Spanish acronym Caja Costarricense de Seguro Social). The basis of the health system in Costa Rica is universal and public, which covers 95% of the Costa Rican population ^[18,19]. The CCSS carries all health care functions in the country and divides the territory into 7 Health Regions (región sanitaria in Spanish), and at the same time, it is divided into 104 Health Areas (AS, from its Spanish acronym áreas de salud) and then, 1045 Basic Team for Comprehensive Health Care (EBAIS, from its Spanish acronym Equipos Básicos de Atención Integral de Salud). Each AS caters to a population ranging from 15,000 to 40,000 individuals in rural areas and 30,000 to 60,000 individuals in urban areas. Subsequently, these ASs are subdivided into 1,045 segments, each overseen by an EBAIS. Each segment, in turn, serves approximately 4,000 people.

The health service in Costa Rica operates on three distinct levels of attention, interconnected through referral and counter-referral mechanisms. The initial level of service acts as the primary entry point to the healthcare system, providing essential services within each EBAIS at the community level. Additionally, there are two possibilities in which a person can access the health system: (1) an insured person can be attended by a private doctor who is registered by the CCSS and access reference to more specialized attention, and (2) an insured person is attended by a doctor hired by the company. The second level supports the first level through a network of 10 main clinics, 13 peripheral hospitals, and 7 regional hospitals. It provides specialized outpatient services, interventions, hospitalization, and medical-surgical treatments for medical specialties and high-demand areas like dermatology, urology, and ophthalmology, all of low complexity. Finally, the third level includes outpatient and inpatient services of greater complexity and specialization that require high technology and level of specialization. It is provided through 3 national hospitals and 6 specialized hospitals; the area of influence of this level covers the territory of several provinces.

Due to this health administrative flow, the accurate diagnostics of respiratory conditions that require hospitalization are registered in hospital discharge events (second level), instead of hospital entrance (first level). As a consequence, the limitation of this data lies in the presence of a lagged and unknown effect, necessitating consideration in the analysis. Beyond the delayed impact associated with hospital admissions and discharges, there is also a temporal lag effect pertaining to climate and pollution variables.

In the literature, most related studies are linked to statistical association between climate and hospitalizations ^{[2,3,4,5,6,20]}. Additionally, some concentrate on forecasting hospitalizations specifically within emergency departments and respiratory diseases ^{[21,22,23,24,25,26,27,28,29]}.

Our main goal is to predict hospital discharges due to respiratory diseases at the secondary healthcare level of the system, in the hope of enhancing the country's healthcare economy. Specifically, by utilizing pollution and climatic indices, we aim to obtain accurate forecasts of respiratory hospitalizations in seven regions by comparing four different predictive models. Our goal is to provide easily implemented models to health authorities in order to make accurate predictions and informed budget decisions in the near future. To the best of our knowledge, this approach is the first predictive study undertaken in the region with these distinctive characteristics. Moreover, this study may serve as a potential tool to be replicated in other countries with similar diverse climatic characteristics around the world.

This paper is divided as follows: Section 2 describes the data and the implemented methodologies. Section 3 presents the analysis and forecasting results. Finally, Section 4 discusses the prediction performance, limitations, and future work.

2.   Materials and methods

This section first provides an overview of the primary datasets used in our study. Then, the statistical methods applied to predict the hospital discharge are outlined.

2.1.   Data

We describe the main variable of interest, respiratory hospitalizations, followed by the input data comprising climate and pollution data, along with their sources.

2.1.1.   Respiratory hospitalization

Weekly data on hospital discharges due to all respiratory diseases (J00-J99 using ICD10 ^[30]) from 2000 to 2019 for seven regions (Brunca, Central Norte, Central Sur, Chorotega, Huetar Atlántica, Huetar Norte and Pacífico Central), were obtained from the Health Statistics Area of the CCSS. The choice of using 'discharge' instead of 'hospital entrance' is based on the accurate registration of diagnostics at the second level of the health system in Costa Rica. According to the official statistics from the CCSS, the average hospitalization duration in regional hospitals in 2019 due to respiratory diseases was 5.97 days ^[31]. For a visual representation of the geographic distribution of these regions, refer to Figure 1.

2.1.2.   Input data: Climate, pollution, and lagged hospital discharge

We consider distinct climatic indices related to precipitation and temperature extremes, as well as the aerosol index and previous weekly hospital discharges for each Health Region as input variables to predict weekly hospital discharges.

Daily precipitation estimates from 2000 to 2019, sourced from the Climate Hazards Group InfraRed Precipitation with Station data (CHIRPs) ^[32], were used to measure the land surface rainfall. Furthermore, the Geophysical Research Center (CIGEFI) of the University of Costa Rica supplied maximum and minimum daily temperatures across the country, employing the same estimation procedure as the Climate Hazards Center Infrared Temperature with Stations (CHIRTS) ^[33], with an extension of the time period until 2019. Both data sources offer high spatial resolution at 5km by 5km.

After reviewing the 27 core indices^* related to climate extremes, as recommended by the Climate and Ocean: Variability, Predictability, and Change (CLIVAR)^† on Climate Change Detection ^[34], we adapted them for tropical countries and computed 22 climate change indices for each region related to cumulative rainfall estimates and temperature for each region $i = 1, ..., 7$ and week $t = 1, ..., T$ (See Table 1).

^*See https://etccdi.pacificclimate.org/list_27_indices.shtml

^†https://www.clivar.org/

Specifically for the temperature indices, we eliminated those indices related to extreme winter conditions (numbers of days in which maximum and minimum temperatures are below 0℃). Moreover, since we are dealing with weekly data, those indices related to the number of days greater than the 90th percentile and below the 10th percentiles, respectively, are zero for the majority of some periods of the year. To maintain the variability of these indices, we modified them to the 75th and 25th percentiles to retain the indices' variability throughout the year. We also incorporated weekly maximum, minimum, and average of daily amplitude, as well as the number of days in which the amplitude is greater than the 75th and 90th percentiles, to account for significant temperature changes within a day.

On the other hand, for precipitation indices, we eliminated the number of consecutive dry days and consecutive wet days since Costa Rica has distinct rainy and dry seasons, and these indices are basically zero in the dry season and 7 days a week in the rainy season. Furthermore, we incorporated the maximum, average, and minimum of maximum daily precipitation across the region, and modified the yearly day count with precipitation greater than the 95th and 99th percentiles to weekly day counts with precipitation greater than the 75th and 90th percentiles to account for the variability of the input variable.

In relation to pollution data, daily measurements of AOD from 2000 to 2019 were obtained from the MODIS Atmosphere L3 Daily Product (MOD08_D3) ^[35] from the National Aeronautics and Space Administration (NASA). We obtained daily data points at 12 locations across the country using a spatial resolution of $1 \times 1$ degree grid. The health regions do not naturally align with the $1\times 1$ grid, making it necessary to interpolate the variables to estimate a representative value for each region at a fixed point in time. Therefore, we calculated the weighted average based on the Euclidean distance between the 12 data points and the centroid of each region $i$. We believe this is a reasonable method for estimating a representative value for each region, making the best possible use of the geographic information available. Finally, weekly AOD values were derived by averaging the daily information.

The selection of the above environmental covariates is based on the results cited in the introduction, where previous studies have found a significant relationship between extreme weather data, pollutants, and hospital admissions. Due to the lack of admission information in the CCSS, we consider the dependent variable of hospital discharges to be a fairly close proxy for the former variable. All precipitation indices were logarithmically transformed due to their asymmetry where the main purpose of this transformation is to ensure stability in calculating the importance factors of these variables, rather than to improve the predictive power of the models. Furthermore, we computed the sample cross-correlation function between all covariates and hospital discharges. We observed that the highest correlation between all covariates and hospital discharges occurs in the same week, that is, at $0$ lag, except for AOD, where the highest sample cross-correlation is detected with a lag of $10$ weeks, namely AOD_10, thus we incorporate this lagged covariate instead of the same-week AOD in the model. Finally, hospital discharges from both the previous week and the week before are used in the model.

2.2.   Statistical Learning Methods

We employed a supervised statistical learning approach to forecast hospital discharges based on the environmental information in each region. The prediction model for hospital discharges at time $t$ ($HD_t$) in a fixed region $i$ is defined as follows:

where $f(\cdot)$ represents the applied machine learning technique, $\boldsymbol{X}_t$ encompasses all 22 climate covariates at time $t$ (see Table 1), and $AOD_{t-10}$ denotes the AOD at time $t-10$. The variables $HD_{t-1}$ and $HD_{t-2}$ represent the lags of one and two weeks respectively of the discharge series. These lags are considered as predictors because they allow to capture the non-stationarity of the original series.

For $f(\cdot)$, we apply four different machine learning methods:

${\bf{Random}} \;{\bf{Forest}}$ (RF) ^[36,37] is a bootstrapped ensemble method that combines results of regression trees to improve the prediction.

${\bf{XGBoost}}$ (XGBOOST) ^[38] is a gradient-boosting algorithm that employs bagging and trains multiple decision trees in order to produce forecasts.

${\bf{Facebook's}}\; {\bf{Prophet}}\; {\bf{ Forecasting }}\; {\bf{Model}}$ (PROPHET) ^[39] is an additive-based model that fits non-linear trends with multiple seasonality plus holiday factors.

${\bf{Ensemble}}\; {\bf{method}}$ (ENSEMBLE) is a combination of the predictions of the three aforementioned models with equal weights.

We selected these methods because they are widely known in the literature and extensively used in the modeling and prediction of admission variables or hospital discharges; see ^{[40,41,42,43,44,45,46,47,48]} for examples of recent uses. We acknowledge the potential for employing hybrid or more complex methodologies (as exemplified in ^[49,50]); however, for this article, we opt for straightforward and purely predictive techniques. We favor simplicity in implementation and interpretation, particularly considering the potential utilization of our methods by non-expert users of the CCSS. For example, any of the models used in this article are implemented in R libraries such as tidymodels ^[51], modeltime ^[52], as well as in scikit-learn ^[53] in Python.

To train the model, we split the data into training (2001-2018) and testing sets (2019). We selected the testing period with the intention of encompassing a complete cycle (annual) of observed behavior in regional discharges, where there was no presence of external shocks such as the COVID-19 pandemic that began in the year 2020. The training process involves a systematic approach, where the training set is divided into sliding windows. Each sliding window comprises a 2-year analysis set and a 1-year assessment set, accounting for the temporal annual seasonality inherent in the data. Following this, we proceed with model tuning utilizing the grid search method employing a Latin hypercube approach with 20 combinations of hyperparameters, all based on the default tuning procedures provided by the $\texttt{tune_grid}$ function in the ${\texttt{tidymodels}}$ package. Subsequently, we refine this process through a grid search over a regular grid with a smaller number of combinations. In both instances, our objective remains consistent: to minimize the mean squared error (MSE). Once we identified the optimized hyperparameters based on both the aforementioned criterion and the tuning procedure (see Table 2), we predict hospital discharges and their prediction interval using conformal inference ^[54] in 2019 to compare with the testing set. The use of conformal inference ensures that prediction intervals are free from any assumptions about the distribution of our outcome, and it also allows for achieving interesting coverage properties even for small samples (see ^[54]).

Finally, we compute three metrics to compare the predictive performance of the applied methods for a fixed region. The first one, MSE ^[55] is defined as follows:

where $m$ is the number of weeks in the testing period, $HD_t$ is the observed hospital discharge at week $t$, and $\widehat{HD}_t$ is the estimated hospital discharge according to each model at time $t$. The Mean Absolute Error ($MAE$) ^[55] is computed as follows:

and finally, the Interval Score (IS) ^[56] is defined as

where $U_t$ and $L_t$ are the upper and lower limits of the prediction interval of $(1-\alpha) \%$, respectively.

MSE and MAE evaluate the point forecasts with the observed hospital discharges. On the other hand, $IS_{\alpha}$ assesses a $(1- \alpha) \%$-prediction interval by comparing its upper and lower limits against the observed values. This metric is more comprehensive than MSE for evaluating the models' predictive capacity when uncertainty is summarized through a predictive interval ^[56]. While it is not the only measure that achieves this, it is one of the simplest to compute and easiest to interpret. Finally, to quantify the degree of overfitting in each of the models, we compare the MSE and the IS obtained in the training set versus those obtained in the testing period through the relative MSE and relative IS as follows:

We also calculate the contribution of each covariate in predicting the dependent variable. In particular, we are interested in quantifying the global contribution through the aggregation of local contributions based on Shapley values. In other words, we calculate the importance of each covariate through its mean absolute Shapley value using the ${\texttt{vip_shap}}$ function from the ${\texttt{vip}}$ R package ^[57]. This allows us to compare which covariates have a greater impact on the prediction process for each specific region.

All calculations were performed in the statistical software R ^[58] and the following packages ${\texttt{tidymodels}}$ ^[51], ${\texttt{modeltime}}$ ^[52], ${\texttt{ranger}}$ ^[59], ${\texttt{xgboost}}$ ^[60] ${\texttt{modeltime.ensemble}}$ ^[61], and ${\texttt{prophet}}$ ^[62]. The block diagram depicted in Figure 2 provides a concise overview of the methodological steps outlined in this article. The code used in this application, as well as the computational details of this article, can be found at https://github.com/shuwei325/pollution_hosp.

3.   Results

In Table 3, we provide a comparison of the four employed methods using the metrics detailed in Section 2. It is crucial to highlight that our objective is to assess the predictive performance of the four models during the testing period. This assessment involves comparing the observed outcomes with the projected ones, utilizing the set of covariates observed within the same period.

Due to differences in the behavior of the series in each region, it is confirmed that there is no method that is optimal for all cases. Likewise, the ensemble method is not optimal according to the calculated metrics. In order to choose the proposed models, we sought to have the smallest IS in the testing period, and we also ensured an $\text{IS}_{rel}$ ratio higher than 80%, with the goal of selecting alternatives with potential overfitting. The other metrics were used to verify that the selected methods are also competitive compared to alternatives. In general terms, all selected methods successfully predict the testing period across the regions and capture both high and low-frequency characteristics of the original discharge series in a descriptive manner, as seen in Figure 3.

Figure 3 illustrates that the 95% predictive intervals successfully capture the observed series in most regions, except for Huetar Atlantica and Norte regions, and partially in the Pacifico Central region, which showed greater difficulty in their forecasts.

For these last regions, the PROPHET method constitutes an interesting alternative, since this method also allows incorporating a time index that can explain those non-stationary components not being explained by the covariates. Its flexibility in modeling trends, seasonal components, and trend changes makes it particularly well-suited for these complex-behaved regions, as it was purposefully designed to accommodate such complexities (see ^[39]). This method has been used in situations where the complexity of the patterns observed in the series makes the use of classical time series models not as effective (see ^[63,64]). For example, the $\text{IS}_{rel}$ reflects the model's ability to provide accurate predictions in the testing set, even in the Huetar Atlantica region where it is slightly outperformed by the RF. This is largely due to its balanced adaptability to the series. Also, in Table 4, we show the Shapley values, all in percentage terms to facilitate interpretation.

It is worth noting that the PROPHET method performs better in three regions (Brunca, Huetar Norte, and Pacifico Central), where the date is the most important feature, while lagged HD and AOD are not. It is worth noting that for other four regions where the RF and XGBOOST methods are the best, the covariates with the greatest contribution to the prediction of hospital discharges tend to be the lagged discharge variables (orders 1 and 2) along with the lagged AOD variable. This indicates that for these regions, the climate and pollution features, as well as lagged HD, are important inputs to predict HD, while for those that PROPHET are best models, temporal information such as nonstationarity and seasonality provides better information to forecast the HD.

On the other hand, depending on the region, there are environmental variables that contribute to predictions to varying degrees once the contribution of lagged variables has been taken into account. For example, it is noteworthy that for the Chorotega and Huetar Norte regions, the amplitude between maximum and minimum temperature is important, as well as in the latter region, the regional minimum of minimum temperature and the regional maximums of maximum temperature.

In general terms, a considerable number of the covariates we used show significant contribution in the prediction models. The difficulty in determining an optimal fitting method may be due to the high variability in weekly discharge data, and the existence of complex seasonality components among different health regions, which may largely stem from the combination of epidemiological cycles of respiratory diseases. These situations of high variability and seasonality in hospital admission data (which we are approximating through discharges) have already been evidenced in previous studies ^{[23,25,65,66]}. It is also worth mentioning that as we are choosing models where IS$_{rel}$ is less than 80%, we then directly minimize the occurrence of overfitting, which as far as we know, was not directly addressed in these same articles previously cited.

4.   Conclusion and recommendations

To summarize, we used climatic and pollution data to model and predict weekly hospital discharges due to respiratory conditions, a leading contributor to the global healthcare burden in terms of costs, at regional hospital levels in Costa Rica. Four statistical learning approaches: RF, XGBOOST, PROPHET, and the ENSEMBLE method, were applied to each of the seven regions. According to the predictive metrics for each region, reliable forecasts are obtained, and specific conclusions can be drawn. In Brunca, Huetar Norte, and the Pacífico Central region, PROPHET obtained better forecasts, and the temporal input was detected as the most important predictor. For other regions (Central Norte, Central Sur, Chorotega, and Huetar Atlántico), where other methods obtained better results, temporal information was not included as a predictor, and other climatic indices and lagged hospitalization and AOD were found to be the most important factors.

These conclusions are obtained from comparing their predictive performance in both the training set (2001-2018) and the testing set (2019). Initially, the training procedure utilized sliding windows of a 2-year training period and a 1-year assessment period, allowing for the optimization of the model's hyperparameters while considering yearly seasonality. Subsequently, their predictive capacity was evaluated in the testing set using metrics such as MSE, MAE, and IS, as well as $\text{MSE}_{rel}$ and $\text{IS}_{rel}$ to assess overfitting. We conclude that the optimal model for each region depends on specific environmental and pollution variables, which do not align for all cases. However, at least for most regions, the lagged variables of hospital discharges and the lagged AOD are of the highest relative importance. In other cases, the non-stationary component is strong enough for a time-indexed model to be preferable. Despite all of the above, the selected model for each region demonstrates its ability to generate reliable forecasts. Although no previous predictive studies were found in Costa Rica, we can conclude that the MSE and MAE metrics of our study are lower compared to a similar study of emergency admissions in the city of Santiago, Chile (see Table 3 in ^[25]), which reports MSE and MAE for their forecasts.

As limitations, it is important to recognize that other significant factors, such as the socio-economic conditions of each region, could serve as potential predictors for hospitalizations due to respiratory conditions. However, many socio-economic factors, including GDP, population density, and healthcare facilities, among others, are not available on a weekly basis or remain constant over time. We argue that, given the established association in the literature between social factors and AOD, we indirectly incorporate this information into the model, recognizing it as a proxy for the anthropogenic effect on climate.

Due to the organizational structure of the CCSS, where the decentralization of services through regional hospitals and clinics by health region is encouraged, decision-making at the hospital level in Costa Rica is crucial. This study enables health authorities to anticipate and visualize the hospital requirements for serving the population with respiratory problems in Costa Rican Health Regions. Additionally, it empowers authorities to administer and allocate limited human and economic resources effectively and efficiently, fostering a balanced national health economy.

As future work, two directions may be considered. First, since our focus is on forecasting respiratory hospitalization using climatic and pollution data, the implicit consideration of nonlinear associations in these models raises questions about the non-linear relationships and long/short-term associations among them, which need further analysis in different microclimates in Costa Rica. Second, while climatic and pollution data are key variables for forecasting respiratory hospitalization, they are not the only factors influencing respiratory health. Socioeconomic data, albeit limited, could be directly or indirectly included to enhance predictions.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Luis A. Barboza and Shu Wei Chou Chen received funding through V.I. C3195, V.I. C4196, and V.I. C3175. We also extend our gratitude to the Health Statistics Area of the CCSS and CIGEFI, UCR, for providing the data used in this manuscript.

Conflict of interest

The authors declare there is no conflict of interest.