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

Shu Wei Chou-Chen
Centro de Investigación en Matematica Pura y Aplicada
Escuela de Estadística
Universidad de Costa Rica
Costa Rica<EMAIL_ADDRESS>
&Luis A. Barboza
Centro de Investigación en Matematica Pura y Aplicada
Escuela de Matemática
Universidad de Costa Rica
Costa Rica<EMAIL_ADDRESS>

###### Abstract

Respiratory diseases represent one of the most significant economic burdens on
healthcare systems worldwide. The variation in the increasing number of cases
depends greatly on climatic seasonal effects, socioeconomic factors, and
pollution. Therefore, understanding these variations and obtaining precise
forecasts allows health authorities to make correct decisions regarding the
allocation of limited economic and human resources. This study aims to model
and forecast weekly hospitalizations due to respiratory conditions in seven
regional hospitals in Costa Rica using four statistical learning techniques
(Random Forest, XGboost, Facebook’s Prophet forecasting model, and an ensemble
method combining the above methods), along with 22 climate change indices and
aerosol optical depth as an indicator of pollution. Models are trained using
data from 2000 to 2018 and are evaluated using data from 2019 as testing data.
Reliable predictions are obtained for each of the seven regional hospitals.

_K_ eywords Temperature, precipitation, climate change, aerosol optical depth,
statistical learning, forecasting, hospitalization, respiratory diseases.

##  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, specifically the five major
conditions, known as the “big five”, which are chronic obstructive pulmonary
disease (COPD), asthma, acute respiratory infections, tuberculosis, and lung
cancer [11]. 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.

Extensive research indicates that climate influences respiratory health,
resulting in higher rates of hospitalization and mortality in varying seasons
[see e.g., 6, 24, 32, 39, 40]. 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 [15, 19, 25, 26,
28].

Specifically with respiratory diseases, climate factors are found to be
associated with extreme weather [27, 39, 40], and pollution, such as
particular matter (PM) [16, 32]. Airborne particulate matter (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[7]. 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 studies 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 [see e.g. 34, 42]. Moreover, AOD is
found to be associated with socio-economic factors, such as GDP, industry, and
vehicle density [23].

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.
[31].

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[31].

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 [1, 31]. 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 serves 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 aims to provide support to the first level
through a network of 10 main clinics, 13 peripheral hospitals, and 7 regional
hospitals, through the provision of specialized outpatient services,
outpatient and inpatient interventions, hospitalization, and medical-surgical
treatment of medical specialties, and some others of high population demand,
such as 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). A
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 existing literature, most studies are linked to statistical association
between climate and hospitalizations [see e.g., 6, 24, 32, 39, 40].
Additionally, some of them concentrate on forecasting hospitalizations
specifically within emergency departments and respiratory diseases [see e.g.
33, 10, 36]. 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, we analyze
and predict hospital discharges in 7 regions to enable health authorities 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.

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

###  2.1 Data

This section provides an overview of the primary datasets used in our study,
focusing on respiratory hospitalization and input data comprising climate,
pollution, and lagged hospital discharge.

####  2.1.1 Respiratory hospitalization

Weekly data on hospital discharges due to all respiratory diseases (J00-J99
using ICD10 [18]) 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 [3]. For a visual representation of the geographic
distribution of these regions, refer to Figure 1.

Figure 1: Health administrative division of Costa Rica.

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

Daily precipitation estimates from 2000 to 2019, sourced from the Climate
Hazards Group InfraRed Precipitation with Station data (CHIRPs) [12], were
used to measure the land surface rainfall. Simultaneously, 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) [13], 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 indices111See
https://etccdi.pacificclimate.org/list_27_indices.shtml related to climate
extremes, as recommended by the Climate and Ocean: Variability, Predictability
and Change (CLIVAR)222 https://www.clivar.org/ on Climate Change Detection
[20], 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).

Table 1:  Climate change indices related to rainfall and temperature. Variable |  Description
---|---
1\. Tmax_max |  Maximum (across the region) of maximum weekly temperature.
2\. Tmax_mean |  Average (across the region) of maximum weekly temperature.
3\. Tmax_min |  Minimum (across the region) of maximum weekly temperature.
4\. n_Tmax_Q3 |  Average (across the region) of the number of days in which the maximum temperature is higher than the percentile 75 of maximum temperature.
5\. Tmin_min |  Minimum (across the region) of minimum weekly temperature.
6\. Tmin_mean |  Average (across the region) of minimum weekly temperature.
7\. Tmin_max |  Maximum (across the region) of minimum weekly temperature.
8\. n_Tmin_Q1 |  Average (across the region) of the number of days in which the minimum temperature is lower than the percentile 25 of minimum temperature.
9\. amplitude_max_max |  Maximum of the maximum daily amplitude of temperature (difference of daily maximum and minimum temperature).
10\. amplitude_max_mean |  Average of the maximum daily amplitude of temperature.
11\. amplitude_max_min |  Minimum of the maximum daily amplitude of temperature.
12\. amplitude_min_max |  Maximum of the minimum daily amplitude of temperature.
13\. amplitude_min_mean |  Average of the minimum daily amplitude of temperature.
14\. amplitude_min_min |  Minimum of the minimum daily amplitude of temperature).
15\. n_amplitude_Q3 |  Average of the number of days in which the amplitude of temperature is higher than the percentile 75 of daily amplitude of temperature.
16\. n_amplitude_P90 |  Average of the number of days in which the amplitude of temperature is lower than the percentile 25 of daily amplitude of temperature.
17\. precip_max_max |  Maximum (across the region) of the maximum daily precipitation.
18\. precip_max_mean |  Average (across the region) of the maximum daily precipitation.
19\. precip_max_min |  Minimum (across the region) of the maximum daily precipitation.
20\. precip_mean_mean |  Average (across the region) of the average daily precipitation Average (across the region) of the average daily precipitation.
21\. n_precip_max_Q3 |  Average (across the region) of the number of days in which the maximum precipitation is higher than the percentile 75 of daily maximum precipitation.
22\. n_precip_max_P90 |  Average (across the region) of the number of days in which the maximum precipitation is higher than the percentile 90 of daily maximum precipitation.

In relation to pollution data, daily measurements of AOD from 2000 to 2019
were obtained from the MODIS Atmosphere L3 Daily Product (MOD08_D3) [29] 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. Subsequently, we calculated the weighted average
based on the Euclidean distance between the 12 data points and the centroid of
each region $i$. Finally, weekly AOD values were derived by averaging the
daily information.

All precipitation indices were logarithmically transformed due to their
asymmetry. 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 (i.e., 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 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:

$HD_{t}\sim f(\boldsymbol{X}_{t},AOD_{t-10},HD_{t-1},HD_{t-2}),$

where $f(\cdot)$ represents the applied machine learning technique,
$\boldsymbol{X}_{t}$ encompasses all climate covariates at time $t$, and
$AOD_{t-10}$ denotes the aerosol optical depth at time $t-10$.

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

Random Forest

(RF) [2, 17] is a bootstrapped ensemble method that combines results of
regression trees to improve the prediction.

XGBoost

(XGBOOST) [4] is a gradient-boosting algorithm that employs bagging and trains
multiple decision trees in order to produce forecasts.

Facebook’s Prophet Forecasting Model

(PROPHET) [38] is an additive-based model that fits non-linear trends with
multiple seasonality plus holiday factors.

Ensemble method

(ENSEMBLE) is a combination of the predictions of the three aforementioned
models with equal weights.

To train the model, we split the data into training (2001-2018) and testing
sets (2019). 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. After that, we perform model tuning
using the grid search method and minimize the root mean squared error (RMSE).
Once we obtain the optimized hyperparameters (see Table 4 in the Supplementary
section), we predict hospital discharges and their prediction interval using
conformal inference [22] 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 [22]).

Finally, we compute three metrics to compare the predictive performance of the
applied methods for a fixed region. The first one, the Mean Squared Error
($MSE$) is defined as follows:

$\displaystyle MSE=\frac{1}{m}\sum_{t=1}^{m}(HD_{t}-\widehat{HD}_{t})^{2},$

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$) is computed as follows:

$\displaystyle MAE=\frac{1}{m}\sum_{t=1}^{m}|HD_{t}-\widehat{HD}_{t})|,$

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

$\displaystyle
IS_{\alpha}=\frac{1}{m}\sum_{t=1}^{m}\left[(U_{t}-L_{t})+\frac{2}{1-\alpha}(L_{t}-HD_{t})\cdot
1_{HD_{t}<L_{t}}+\frac{2}{1-\alpha}(HD_{t}-U_{t})\cdot
1_{HD_{t}>U_{t}}\right],$

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. 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:

$MSE_{rel}=\frac{MSE_{test}}{MSE_{train}},~{}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}IS_{rel}=\frac{IS_{test}}{IS_{train}}.$

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 [for more details, see 30, 43]. 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 [35] and the
following packages `tidymodels` [21], `modeltime` [8], `ranger` [41] ,
`xgboost` [5] `modeltime.ensemble` [9], and `prophet` [37].

##  3 Results

In Table 2, 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.

Table 2: Metric comparison by region and method (all covariates). The best method for each region is shown in bold. Region | Model | $MSE$ | $MAE$ | $IS_{0.05}$ |  $MSE$(train) |  $MAE$(train) |  $IS_{0.05}$(train) | $MSE_{rel}$ | $IS_{rel}$
---|---|---|---|---|---|---|---|---|---
Brunca | RF | 188.55 | 10.18 | 71.79 | 73.10 | 7.00 | 60.09 | 68.78 | 83.69
XGBOOST | 253.34 | 11.74 | 79.12 | 74.65 | 6.94 | 69.71 | 59.13 | 88.10
PROPHET | 237.86 | 13.29 | 63.68 | 55.41 | 6.08 | 59.29 | 45.77 | 93.10
| ENSEMBLE | 205.36 | 11.29 | 70.83 | 58.74 | 6.25 | 34.67 | 55.32 | 48.94
Central Norte | RF | 222.95 | 12.01 | 69.65 | 97.77 | 7.80 | 56.27 | 64.89 | 80.79
XGBOOST | 281.86 | 13.70 | 75.23 | 0.00 | 0.03 | 71.80 | 0.21 | 95.44
PROPHET | 2947.93 | 48.38 | 210.24 | 99.48 | 7.80 | 191.61 | 16.12 | 91.14
| ENSEMBLE | 540.34 | 18.11 | 112.19 | 38.07 | 4.89 | 29.89 | 27.00 | 26.64
Central Sur | RF | 407.37 | 16.63 | 85.27 | 166.16 | 10.01 | 76.44 | 60.21 | 89.64
XGBOOST | 429.42 | 17.28 | 86.63 | 0.00 | 0.03 | 79.24 | 0.17 | 91.47
PROPHET | 660.35 | 21.00 | 110.24 | 131.06 | 9.23 | 89.29 | 43.97 | 81.00
| ENSEMBLE | 408.17 | 17.28 | 89.01 | 57.41 | 6.14 | 35.92 | 35.52 | 40.36
Chorotega | RF | 90.16 | 7.34 | 44.12 | 24.29 | 3.87 | 41.31 | 52.71 | 93.64
XGBOOST | 84.14 | 7.11 | 42.12 | 32.86 | 4.49 | 42.01 | 63.18 | 99.75
PROPHET | 96.64 | 8.09 | 47.00 | 29.61 | 4.22 | 35.05 | 52.15 | 74.58
| ENSEMBLE | 81.11 | 7.18 | 42.19 | 25.60 | 3.94 | 24.81 | 54.92 | 58.81
Huetar Atlántica | RF | 142.02 | 8.77 | 58.77 | 32.12 | 4.36 | 57.38 | 49.73 | 97.63
XGBOOST | 156.60 | 9.13 | 62.70 | 34.71 | 4.50 | 58.33 | 49.28 | 93.03
PROPHET | 237.39 | 12.56 | 66.79 | 22.09 | 3.72 | 63.98 | 29.61 | 95.79
| ENSEMBLE | 160.25 | 9.61 | 61.59 | 25.67 | 3.90 | 24.55 | 40.55 | 39.87
Huetar Norte | RF | 158.97 | 10.13 | 55.13 | 27.52 | 4.12 | 49.76 | 40.70 | 90.26
XGBOOST | 180.60 | 10.85 | 60.19 | 48.21 | 5.42 | 54.86 | 49.94 | 91.13
PROPHET | 163.56 | 10.15 | 58.24 | 34.34 | 4.72 | 55.85 | 46.52 | 95.89
| ENSEMBLE | 161.28 | 10.20 | 56.61 | 33.23 | 4.51 | 27.50 | 44.21 | 48.58
Pacífico Central | RF | 53.64 | 5.89 | 34.90 | 35.48 | 4.73 | 31.01 | 80.21 | 88.87
XGBOOST | 47.99 | 5.46 | 32.65 | 56.81 | 5.83 | 35.97 | 106.65 | 110.16
PROPHET | 85.38 | 8.18 | 32.11 | 34.52 | 4.67 | 32.05 | 57.06 | 99.82
| ENSEMBLE | 55.21 | 6.20 | 31.03 | 37.39 | 4.82 | 28.68 | 77.67 | 92.42

Note that 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 $IS_{rel}$ ratio higher than 80%,
with the goal of selecting alternatives with minimal 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 in the testing period across the regions and, furthermore, allow
capturing both high and low-frequency characteristics in the original
discharge series.

Figure 2: Forecast comparison of hospital discharges, by region.

In Figure 2, we present these results. Note that in most regions, the 95%
confidence intervals successfully capture the observed series as well. The
regions that showed greater difficulty in their forecasts were the Huetar
Atlantica and Norte regions, and partially the Pacifico Central region.

Table 3: Relative Importance factors for the best methods (in percentage) by means of the Shapley values. The date covariate corresponds to the time index, needed for the PROPHET algorithm. Covariate | Brunca |  Central Norte |  Central Sur | Chorotega |  Huetar Atlantica |  Huetar Norte |  Pacifico Central
---|---|---|---|---|---|---|---
date | 21.1 |  |  |  |  |  18.4 |  23.2
n_amplitude_Q3 | 11.0 |  |  2.0 |  |  4.4 |  7.5 |  7.6
n_amplitude_P90 | 9.9 |  |  |  |  6.7 |  |
amplitude_max_max | 9.9 |  |  |  |  |  6.1 |  7.4
egreso_2 | 9.3 |  21.6 |  18.9 | 36.1 |  31.0 |  |  8.8
precip_max_min | 9.1 |  |  | 3.1 |  3.1 |  |  6.6
precip_mean_mean | 9.0 |  3.0 |  3.0 |  |  |  7.7 |
Tmax_max | 7.9 |  |  |  |  |  13.1 |
n_precip_max_Q3 | 6.6 |  2.3 |  | 3.0 |  3.6 |  |
Tmax_mean | 6.1 |  |  1.5 |  |  4.0 |  10.9 |
egreso_1 |  |  46.5 |  50.3 | 16.8 |  24.5 |  |
aerosol_10 |  |  15.6 |  16.5 | 17.1 |  14.0 |  |
aerosol |  |  4.3 |  | 4.5 |  |  |
precip_max_max |  |  2.2 |  1.6 | 3.8 |  |  |
precip_max_mean |  |  1.6 |  2.3 |  |  |  6.1 |
Tmin_min |  |  1.5 |  | 2.8 |  |  7.3 |  7.5
Tmin_mean |  |  1.4 |  |  |  |  13.8 |  13.0
amplitude_min_min |  |  |  2.1 |  |  |  |
amplitude_min_mean |  |  |  1.8 |  |  |  |  5.7
amplitude_max_min |  |  |  | 9.3 |  |  9.2 |
n_Tmin_Q1 |  |  |  | 3.5 |  |  |
amplitude_max_mean |  |  |  |  |  5.8 |  |  10.4
amplitude_min_max |  |  |  |  |  2.9 |  |
Tmin_max |  |  |  |  |  |  |  10.0

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.
Also, in Table 3, we show the Shapley values, all in percentage terms to
facilitate interpretation. Note that for those 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. 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.

##  4 Discussion

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. Furthermore,
we compared their predictive capacity in both the training set (2001-2018) and
the testing set (2019) using MSE, MAE, and IS, as well as $MSE_{rel}$ and
$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 the majority of regions, the lagged
variables of hospital discharges and the lagged aerosol variable 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.

As limitations, it is important to note 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.

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.

## References

  * [1] M. Baird.  (In)equity and primary health care: The case of Costa Rica and Panama.  International Journal of Social Determinants of Health and Health Services, 53(2):122–129, 2023.  PMID: 36726329.
  * [2] L. Breiman.  Random forests.  Machine learning, 45(1):5–32, 2001.
  * [3] CCSS.  Estadística en salud.  Available online at https://www.ccss.sa.cr/estadisticas-salud (accessed on October 16, 2023).
  * [4] T. Chen and C. Guestrin.  Xgboost: A scalable tree boosting system.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, page 785–794, New York, NY, USA, 2016. Association for Computing Machinery.
  * [5] T. Chen, T. He, M. Benesty, V. Khotilovich, Y. Tang, H. Cho, K. Chen, R. Mitchell, I. Cano, T. Zhou, M. Li, J. Xie, M. Lin, Y. Geng, Y. Li, and J. Yuan.  xgboost: Extreme Gradient Boosting, 2023.  R package version 1.7.5.1.
  * [6] Y. Chen, D. Kong, J. Fu, Y. Zhang, Y. Zhao, Y. Liu, Z. Chang, Y. Liu, X. Liu, K. Xu, C. Jiang, and Z. Fan.  Associations between ambient temperature and adult asthma hospitalizations in Beijing, China: a time-stratified case-crossover study.  Respiratory Research, 23(1):38, Feb. 2022.
  * [7] I. Colbeck and M. Lazaridis.  Aerosols and environmental pollution.  Naturwissenschaften, 97(2):117–131, Feb. 2010.
  * [8] M. Dancho.  modeltime: The Tidymodels Extension for Time Series Modeling, 2023.  R package version 1.2.8.
  * [9] M. Dancho.  modeltime.ensemble: Ensemble Algorithms for Time Series Forecasting with Modeltime, 2023.  R package version 1.0.3.
  * [10] Y. de Souza Tadano, H. V. Siqueira, and T. A. Alves.  Unorganized machines to predict hospital admissions for respiratory diseases.  In 2016 IEEE Latin American Conference on Computational Intelligence (LA-CCI), pages 1–6, 2016.
  * [11] Forum of International Respiratory Societies.  Respiratory diseases in the world. realities of today – opportunities for tomorrow.  Technical Report 9, Mar. 2014.
  * [12] C. Funk, P. Peterson, M. Landsfeld, D. Pedreros, J. Verdin, S. Shukla, G. Husak, J. Rowland, L. Harrison, A. Hoell, and J. Michaelsen.  The climate hazards infrared precipitation with stations–a new environmental record for monitoring extremes.  Scientific data, 2(1):150066–150066, 2015.
  * [13] C. Funk, P. Peterson, S. Peterson, S. Shukla, F. Davenport, J. Michaelsen, K. R. Knapp, M. Landsfeld, G. Husak, L. Harrison, J. Rowland, M. Budde, A. Meiburg, T. Dinku, D. Pedreros, and N. Mata.  A high-resolution 1983–2016 tmax climate data record based on infrared temperatures and stations by the climate hazard center.  Journal of Climate, 32(17):5639 – 5658, 2019.
  * [14] T. Gneiting and A. E. Raftery.  Strictly Proper Scoring Rules, Prediction, and Estimation.  Journal of the American Statistical Association, 102(477):359–378, mar 2007.
  * [15] Y. Guo, A. Gasparrini, B. G. Armstrong, B. Tawatsupa, A. Tobias, E. Lavigne, M. d. S. Z. S. Coelho, X. Pan, H. Kim, M. Hashizume, Y. Honda, Y. L. Guo, C.-F. Wu, A. Zanobetti, J. D. Schwartz, M. L. Bell, A. Overcenco, K. Punnasiri, S. Li, L. Tian, P. Saldiva, G. Williams, and S. Tong.  Temperature Variability and Mortality: A Multi-Country Study.  Environmental Health Perspectives, 124(10):1554–1559, Oct. 2016.
  * [16] G. B. Hamra, N. Guha, A. Cohen, F. Laden, O. Raaschou-Nielsen, J. M. Samet, P. Vineis, F. Forastiere, P. Saldiva, T. Yorifuji, and D. Loomis.  Outdoor particulate matter exposure and lung cancer: a systematic review and meta-analysis.  Environmental Health Perspectives, 122(9):906–911, Sept. 2014.
  * [17] T. Hastie, R. Tibshirani, and J. Friedman.  The Elements of Statistical Learning: Data Mining, Inference, and Prediction.  Springer Series in Statistics. Springer, 2nd edition, 2009.
  * [18] ICD-10-CM.  J00-J99: Diseases of the respiratory system.  Available online at https://www.icd10data.com/ICD10CM/Codes/J00-J99 (accessed on September 10, 2022).
  * [19] H. Kan, R. Chen, and S. Tong.  Ambient air pollution, climate change, and population health in China.  Environment International, 42:10–19, 2012.
  * [20] T. R. Karl, N. Nicholls, and A. Ghazi.  CLIVAR/GCOS/WMO Workshop on Indices and Indicators for Climate Extremes Workshop Summary, pages 3–7.  Springer Netherlands, Dordrecht, 1999.
  * [21] M. Kuhn and H. Wickham.  Tidymodels: a collection of packages for modeling and machine learning using tidyverse principles., 2020.
  * [22] J. Lei, M. G’Sell, A. Rinaldo, R. J. Tibshirani, and L. Wasserman.  Distribution-free predictive inference for regression.  Journal of the American Statistical Association, 113(523):1094–1111, 2018.
  * [23] L. Li and Y. Wang.  What drives the aerosol distribution in Guangdong - the most developed province in Southern China?  Scientific Reports, 4(1):5972, Aug. 2014.  Number: 1 Publisher: Nature Publishing Group.
  * [24] S. Lin, M. Luo, R. J. Walker, X. Liu, S.-A. Hwang, and R. Chinery.  Extreme high temperatures and hospital admissions for respiratory and cardiovascular diseases.  Epidemiology (Cambridge, Mass.), 20(5):738–746, Sept. 2009.
  * [25] C. Liu, R. Chen, F. Sera, A. M. Vicedo-Cabrera, Y. Guo, S. Tong, M. S. Coelho, P. H. Saldiva, E. Lavigne, P. Matus, N. Valdes Ortega, S. Osorio Garcia, M. Pascal, M. Stafoggia, M. Scortichini, M. Hashizume, Y. Honda, M. Hurtado-Díaz, J. Cruz, B. Nunes, J. P. Teixeira, H. Kim, A. Tobias, C. Íñiguez, B. Forsberg, C. Åström, M. S. Ragettli, Y.-L. Guo, B.-Y. Chen, M. L. Bell, C. Y. Wright, N. Scovronick, R. M. Garland, A. Milojevic, J. Kyselý, A. Urban, H. Orru, E. Indermitte, J. J. Jaakkola, N. R. Ryti, K. Katsouyanni, A. Analitis, A. Zanobetti, J. Schwartz, J. Chen, T. Wu, A. Cohen, A. Gasparrini, and H. Kan.  Ambient Particulate Air Pollution and Daily Mortality in 652 Cities.  New England Journal of Medicine, 381(8):705–715, Aug. 2019.  Publisher: Massachusetts Medical Society _eprint: https://doi.org/10.1056/NEJMoa1817364.
  * [26] J. Lou, Y. Wu, P. Liu, S. H. Kota, and L. Huang.  Health Effects of Climate Change Through Temperature and Air Pollution.  Current Pollution Reports, 5(3):144–158, Sept. 2019.
  * [27] F. Makrufardi, A. Manullang, D. Rusmawatiningtyas, K. F. Chung, S.-C. Lin, and H.-C. Chuang.  Extreme weather and asthma: a systematic review and meta-analysis.  European Respiratory Review, 32(168), June 2023.  Publisher: European Respiratory Society Section: Reviews.
  * [28] S. G. E. K. Miraglia, P. H. N. Saldiva, and G. M. Böhm.  An Evaluation of Air Pollution Health Impacts and Costs in São Paulo, Brazil.  Environmental Management, 35(5):667–676, May 2005.
  * [29] MODIS Atmosphere Science Team.  Modis/terra aerosol cloud water vapor ozone daily l3 global 1deg cmg, 2017.
  * [30] C. Molnar.  Interpretable machine learning. a guide for making black box models explainable.  Book, 2020.
  * [31] Organización Panamericana de la Salud.  Perfil del sistema y servicios de salud de Costa Rica con base al marco de monitoreo de la Estrategia Regional de Salud Universal.  Feb. 2019.  Place: San José.
  * [32] H. Pedder, T. Kapwata, G. Howard, R. N. Naidoo, Z. Kunene, R. W. Morris, A. Mathee, and C. Y. Wright.  Lagged Association between Climate Variables and Hospital Admissions for Pneumonia in South Africa.  International Journal of Environmental Research and Public Health, 18(12), June 2021.  Publisher: Multidisciplinary Digital Publishing Institute (MDPI).
  * [33] J. Peng, C. Chen, M. Zhou, X. Xie, Y. Zhou, and C.-H. Luo.  Peak outpatient and emergency department visit forecasting for patients with chronic respiratory diseases using machine learning methods: Retrospective cohort study.  JMIR Medical Informatics, 8(3):e13075, 3 2020.
  * [34] S. Provençal, P. Kishcha, A. M. d. Silva, E. Elhacham, and P. Alpert.  AOD distributions and trends of major aerosol species over a selection of the world’s most populated cities based on the 1st Version of NASA’s MERRA Aerosol Reanalysis.  Urban climate, 20:168, June 2017.  Publisher: NIH Public Access.
  * [35] R Core Team.  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria, 2023.
  * [36] I. N. Soyiri, D. D. Reidpath, and C. Sarran.  Forecasting asthma-related hospital admissions in london using negative binomial models.  Chronic Respiratory Disease, 10(2):85–94, 2013.  PMID: 23620439.
  * [37] S. Taylor and B. Letham.  prophet: Automatic Forecasting Procedure, 2021.  R package version 1.0.
  * [38] S. J. Taylor and B. Letham.  Forecasting at scale.  American Statistician, 72, 2018.
  * [39] Q. A. Tran, V. T. H. Le, V. T. Ngo, T. H. Le, D. T. Phung, J. D. Berman, and H. L. T. Nguyen.  The Association Between Ambient Temperatures and Hospital Admissions Due to Respiratory Diseases in the Capital City of Vietnam.  Frontiers in Public Health, 10, 2022.
  * [40] H. Tsangari, A. K. Paschalidou, A. P. Kassomenos, S. Vardoulakis, C. Heaviside, K. E. Georgiou, and E. N. Yamasaki.  Extreme weather and air pollution effects on cardiovascular and respiratory hospital admissions in Cyprus.  Science of The Total Environment, 542:247–253, Jan. 2016.
  * [41] M. N. Wright and A. Ziegler.  ranger: A fast implementation of random forests for high dimensional data in C++ and R.  Journal of Statistical Software, 77(1):1–17, 2017.
  * [42] W. Zhang, Q. He, H. Wang, K. Cao, and S. He.  Factor analysis for aerosol optical depth and its prediction from the perspective of land-use change.  Ecological Indicators, 93:458–469, Oct. 2018.
  * [43] E. Štrumbelj and I. Kononenko.  Explaining prediction models and individual predictions with feature contributions.  Knowledge and Information Systems, 41, 2014.

## Supplementary

Table 4: Optimal hyperparameters found using grid search for each method and region. For the RF and XGBOOST methods, the number of trees was held constant at 1000. Method | Hyperparameter | Brunca |  Central Norte |  Central Sur | Chorotega |  Huetar Atlantica |  Huetar Norte |  Pacifico Central
---|---|---|---|---|---|---|---|---
RF | mtry | 25 |  20 |  17 | 15 |  12 |  20 |  22
min_n | 40 |  22 |  30 | 30 |  40 |  40 |  37
XGBOOST | mtry | 12 |  14 |  14 | 12 |  12 |  10 |  10
min_n | 40 |  8 |  8 | 40 |  40 |  35 |  35
tree_depth | 13 |  11 |  11 | 13 |  13 |  12 |  12
learn_rate | 0.007 |  0.078 |  0.078 | 0.007 |  0.007 |  0.003 |  0.003
loss_reduction | 0.001 |  0.018 |  0.018 | 0.001 |  0.001 |  0.000 |  0.000
sample_size | 0.898 |  0.828 |  0.828 | 0.898 |  0.898 |  0.589 |  0.589
PROPHET | prior_scale_changepoints | 1.000 |  1.778 |  1.000 | 1.000 |  3.162 |  1.778 |  1.000
prior_scale_seasonality | 3.162 |  1.778 |  10.000 | 1.000 |  1.778 |  5.623 |  1.000