This study utilized clinical records from the daily emergency department admissions of the Polyclinic hospital in the city of Bari during the reference period from 2013–2021. The database of daily hospitalizations recorded under the “main problem” field included patient cases related to the presented pathology upon their arrival at the emergency department. This database was compiled based on 33 codes, representing the specific pathologies as listed in Table 1.

The first phase of data processing and reorganization was carried out by dividing the incoming information by the various examined years. The statistics presented below refer to the years 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, and 2021. In particular, for the year 2013, there were 75,927 emergency department admissions, including 40,265 male patients, 35,032 female patients, and 630 cases of missing data regarding sex. For the year 2014, there were 80,690 total admissions, including 42,554 male patients, 37,127 female patients, and 1009 cases of missing data regarding sex. For 2015, there were 75,334 total admissions, including 40,091 male patients, 34,327 female patients, and 916 cases of missing data regarding sex. For the year 2016, there were 71,550 total admissions, including 38,007 male patients, 32,914 female patients, and 629 cases of missing data regarding sex. For 2017, there were 65,984 total admissions, including 35,130 male patients, 30,343 female patients, and 511 cases of missing data regarding sex. For 2018, there were 65,641 total emergency department admissions, including 34,798 male patients, 30,544 female patients, and 299 cases of missing data regarding sex. For 2019, there were 68,052 total admissions, including 36,034 male patients, 31,442 female patients, and 576 cases of missing data regarding sex. For 2020, there were 43,729 total admissions, including 24,517 male patients, 18,960 female patients, and 252 cases of missing data regarding sex. In 2021, the total number of Emergency Room (ER) admissions amounted to 48,489, including 26,904 male patients, 21,355 female patients, and 230 cases of missing data regarding sex (Table 2, Fig. 1).

For the purpose of this study, only data and statistics both related to cardiovascular pathologies and are strongly correlated with meteorological factors (Table 3).

Analysis of the effects of hot and cold ambient temperatures on mortality and morbidity in the elderly ($>$65 years old) showed that a 1 °C increase in temperature increased cardiovascular mortality (3.44%), while a 1 °C decrease in temperature increased cardiovascular mortality (1.66%) [23]. Therefore, out of the 33 types of pathologies described in Table 1, only those related to CVDs were selected according to the grouping scheme reported in Table 3 and analyzed for the various years examined, as shown in Fig. 2. The effects of climate variations may have different effects on each individual. A first step in identifying whether there are substantial differences among different types of individuals is gender analysis. Table 4 and Fig. 3 show the annual Emergency Room admissions for cardiovascular pathologies by gender from 2013 to 2021. A further and important distinction is made according to the age of the patient. The distribution of cardiovascular admissions in the emergency department is shown in Tables 5a,5b and in Fig. 4 as a percentage of total admissions.

As shown in Fig. 4, for cardiovascular diseases, the age group between 40 and 54 years has the highest number of hospitalizations.

The meteo-climatic parameters considered in this study include: daily mean minimum temperature (Tmin), daily mean temperature (Tmean), daily mean maximum temperature (Tmax), daily mean dew point temperature (Tdewp), daily mean apparent temperature (Tapp), daily mean atmospheric pressure (Patm), and daily mean relative humidity (RH). Air quality parameters considered include CO (carbon monoxide), O${}_3$ (ozone), PM10 (particulate matter), SO${}_2$ (sulfur dioxide), and NO${}_2$ (nitrogen dioxide) (Table 6). The meteorological data for the city of Bari for the reference period 2013–2021 were obtained from the Arpa Puglia website and the Meteonetwork measurement network. Arpa Puglia manages two networks for monitoring activities:

$\bullet$ A dedicated network: 5 automatic stations located at its provincial offices (Bari, Brindisi, Foggia, Lecce, and Taranto).

$\bullet$ A meteorological network supporting the air quality monitoring network currently consisting of 19 stations.

The dedicated Bari station also monitors ultraviolet radiation [33]. Air quality data for the city of Bari from 2013 to 2021 was obtained from the Arpa Puglia website, through the Bari-Caldarola, Bari-CUS, Bari-Kennedy, Bari-Carbonara monitoring stations, and the mobile laboratory. Meteorological data are recorded with a half-hourly frequency for the Arpa Puglia weather stations and with a daily frequency of five minutes and one hour for the Meteonetwork measurement network. The hourly frequency concerns only the years 2020 and 2021. Air quality data were recorded on a daily basis. The apparent temperature and SO${}_2$ (sulfur dioxide) were excluded from the meteorological and air quality parameters, respectively, due to a lack of representative data, which is insufficient for a complete and accurate elaboration of the same.

The primary objective of the study was to determine which meteorological and air quality variables have the greatest impact on emergency department admissions for CVDs. The methodology adopted from Telesca et al. [32] (Fig. 5). The first step was to find correlations between meteorological variables and emergency department admissions for CVDs through correlation analysis (see Section 4.2), calculating the Pearson coefficient “r” and p-value. Subsequently, the Random Forest model was applied along with its corresponding metrics to evaluate the model’s performance. If the analysis produced acceptable results (Mean Absolute Error [MAE], coefficient of determination R${}^2$, and error analysis) further analysis was performed to determine the most significant meteorological variables (see Section 4.4). If the results were not acceptable, a data decomposition model (Seasonal and Trend decomposition using Loess (STL) - see Section 4.3) is applied to extract the trend component from the data series. Then, the correlation analysis is repeated. This iterative process eventually leads to the most significant variables for the final correlation.

Random Forest, a powerful ensemble learning method, is built upon decision trees, which are hierarchical structures used for making sequential decisions in assigning data points to classes or predicting continuous values in regression tasks. By combining the predictions of multiple models, Random Forest achieves accurate regression predictions. Each model in the ensemble represents a decision tree, and together, they leverage collective knowledge to enhance the overall predictive capability. To construct decision trees, the algorithm employs bootstrap aggregating, also known as bagging. This technique involves creating multiple bootstrap samples from the original training data. Each sample is obtained by randomly selecting data points with replacement, forming subsets used to train individual decision trees. In addition to bagging, Random Forest introduces randomness in feature selection during decision tree construction. Instead of considering all features at each node, a random subset of features is chosen for splitting. This approach ensures that each decision tree in the ensemble learns and makes predictions based on different aspects of the data, leading to a diverse and robust ensemble. When generating predictions with a Random Forest regression model, the ensemble combines the individual decision tree predictions. A commonly employed approach is to compute the average of the predicted values across all trees. This averaging technique helps mitigate the impact of outliers and noise, yielding more reliable predictions. Random Forest provides a convenient mechanism for estimating the model’s performance without requiring a separate validation set. This mechanism relies on out-of-bag (OOB) samples, which are data points that were left out in each bootstrap sample. By comparing the predictions of these OOB samples with their actual target values, one can evaluate the model’s accuracy. Moreover, Random Forest offers insights into the relative importance of different features in the regression task. By analyzing the impact of feature splits across the ensemble of decision trees, the model calculates variable importance scores. These scores indicate which features contribute the most to the model’s predictive power, providing interpretability and understanding of the underlying data relationships. In brief, there are various theoretical advantages associated with using Random Forest techniques for regression models. These include ensemble learning using decision trees, bagging for robust tree construction, random feature selection to ensure diversity, and the ability to estimate model performance and assess variable importance. Collectively, these aspects contribute to the overall effectiveness, accuracy, and interpretability of Random Forest as a powerful regression modeling technique.

The initial step of the process involved removing records with missing data from the database to ensure data integrity. This was followed by searching for the trends in hospital admissions between 2013 to 2021 to account for possible interference caused by anomalous variations in ER visits during the pandemic years of 2020 and 2021. The exclusion of data from 2020 and 2021 was based on the analysis showing distinct differences in ER visits during those years, as depicted in Fig. 6, which presents the 7-day moving averages normalized with the minimum/maximum (min/max) method to ensure comparability across different years. Thus, only data from the years 2013 to 2019 were considered for further analysis.

Correlation analysis is a bivariate statistical technique that measures the strength of a linear relationship between two variables and calculates this relationship. In order to express quantitatively the intensity of the link between two variables, it is necessary to calculate a correlation coefficient [34]. There are various types of correlation coefficients, but all have some common features such as the value that oscillates between –1 and +1, representing a perfect relationship between variables. While a value of 0 indicates the absence of a relationship [35]. To conduct this analysis, a correlation matrix was used. The correlation matrix is a table in which each cell shows the correlation between two variables. The matrix calculates the Pearson correlation coefficient “r”, one of the most widely used correlation coefficients, for each pair of characteristics.

Correlation analysis was conducted utilizing the p-value as a tool for significance testing of the null hypothesis. The p-value represents the probability of obtaining a specific set of observed values assuming the null hypothesis is true, indicating the correctness of our assertion with minimal error. A significantly small p-value indicates that an extreme observed result would be highly improbable under the null hypothesis. This statistical measure establishes the reliability of the correlation values obtained. Acceptable hypotheses for the input variable set were defined as having p-values less than 0.01 and r-values greater than or equal to 0.45 [36]. In the context considered, none of the input variables of the model meet the above conditions, as evidenced by the correlation matrix (Fig. 7). Therefore, the decomposition model will be used.

Many forecasting methods are based on the idea that if there is a systematic pattern, it can be identified and separated from any random fluctuations by smoothing methods of the historical series data. The smoothing effect is the removal of random disturbances, once the pattern is known, in order to project it into the future for forecasting purposes. Decomposition models are mainly divided into additive time series models and multiplicative time series models. The additive model assumes that the effects of each component are independent of each other and that each component is expressed in absolute terms. An additive model is appropriate when the amplitude of the seasonal oscillation does not vary with the level of the series. The error can take positive or negative values, while the neutral value is expressed with the value 0, meaning it does not affect the series. The multiplicative model assumes that the effects of each component on the evolution of the phenomenon are interrelated based on the absolute magnitude of the trend component and that the other components are expressed proportionally. A multiplicative model is suitable when the seasonal fluctuation varies, increases, or decreases proportionally with the variation of the series level. The error can only take non-negative values and has a neutral value of 1. STL is a statistical method that allows the decomposition of time series data into three components: trend, seasonality, and residual. The trend component reflects the long-term variation of the series. A trend exists when there is a persistent increasing or decreasing direction in the data. The seasonal component reflects seasonality, i.e., the variation of data that occurs at specific regular intervals of one year or less. The residual component describes casual or irregular influences. It represents the residuals of the series after the other components have been removed [37].

Let $Y_v$, $T_v$, $S_v$e $R_v$, for v = 1 to n, represent the data, trend component, seasonal component, and residual component, respectively (Eqn. 1).

(1)$$\begin{array}{c}\hfill Y_v=T_v+S_v+R_v\hfill \end{array}$$

The STL consists of a sequence of smoothing operations, each of which, with one exception, employs the same smoother: locally weighted regression, or locally estimated scatterplot smoothing (Loess). Local regression or Loess involves determining, for each point in the initial data set, the coefficients of a low-degree polynomial to perform regression on a subset of data and calculate the value of this polynomial for the given point. The coefficients of the polynomial are computed using weighted least squares, which gives more weight to points near the point whose response is being estimated and less weight to more distant points. This method is used to fit the entire curve and accurately decompose the time series [38]. Fig. 8 shows an example of the decomposition, using the STL model, of the daily average temperature from 2013 to 2021 examined in this thesis. Four graphs are shown:

$\bullet$ The first is the original series;

$\bullet$ The second is the trend;

$\bullet$ The third is the seasonality;

$\bullet$ The fourth is the residual.

As for our study, the STL decomposition model was applied to the original data after verifying the stationarity of the series using the Dickey-Fuller method [39]. The Dickey-Fuller test is a statistical method used to check whether a time series is stationary or not, i.e., if its statistical properties do not change over time. Thus, a modified series was obtained, represented by the trend values of the original data. Then, by recalculating the correlation matrix, significantly higher Pearson coefficient values were obtained compared to the original data series (Fig. 9). These results allow the application of machine learning models to produce an effective simulation model of the trend of hospital admissions over time.

Random Forest is the AI model used to simulate the trend of emergency department visits for the two considered pathologies, and was also estimated the most influential meteorological variables. The Random Forest is a supervised machine learning algorithm, a special classifier consisting of a set of simple classifiers, called decision trees, represented as independent and identically distributed random vectors. Decision trees are common supervised learning algorithms. Decision trees start with a basic question to determine an answer, and these questions make up the decision nodes of the tree that serve to divide the data. Each question helps to make a final decision, which is indicated by the leaf node. When multiple decision trees form an ensemble in the Random Forest algorithm, they predict more accurate results, particularly when the individual trees are not correlated with each other. An ensemble refers to a set of decision trees, and their predictions are aggregated to identify the most popular result [40]. The Random Forest combines the results of all these decision trees to obtain a single result. Each decision tree is composed of a data sample drawn from a training set with replacement, which means that individual data points can be chosen more than once. From this training sample, one-third is set aside as test data, also known as an OOB sample. Another instance of randomness is injected through feature begging or the random subspace method, adding greater diversity to the data set and reducing the correlation between decision trees. Depending on the type of problem, the determination of the forecast varies. For a regression task, the individual decision trees will be averaged, while for a classification task, the majority vote will give the prediction class. Finally, the test sample is used for cross-validation, finalizing the prediction. Cross-validation is a statistical technique that allows the use of both training and test data alternately. There are several cross-validation methods, and the K-fold was used in this paper. K-fold divides the data into k different subsets and k-1 are used for the training phase and the last, remaining subset, for the test phase. The error is then calculated on the observations of the subsets kept out. This procedure is repeated k times by choosing a different subset and obtaining k estimates of the test error. The final estimate will be an average of these values [41].

Table 7 reports the values of the MAE and R${}^2$ metrics for the sample extracted during training and testing phases, as well as the same metrics in the average value of the cross-validation subsets.

The values obtained through cross-validation are analogous to those determined for the random sample, highlighting the model’s ability to avoid overfitting phenomena, demonstrating the accuracy of the model. Subsequently, in Fig. 10, it is shown in the left graph how the trend of the data concerning cardiovascular pathologies predicted by the Random Forest model, described in blue, has an identical trend to the actual data.

The accuracy of the model is highlighted in Fig. 10. The blue points represent both the real and simulated data, and for the most part, they fall within the error bands (Fig. 11). Most data points fall within $\pm$ 10% of the error bands and only a few values fall within $\pm$ 20% of the error bands. In particular, 2443 out of 2542 data points (total number of data points, 96.10% of the total number of ER admissions), fall within the $\pm$ 5% error band (Table 8).

Since the Machine Learning simulation model has been shown to provide accurate results, we can proceed with the use of Feature Importance techniques through the SHapley Additive exPlanations (SHAP) method to identify the most influential characteristics in determining the trend of daily hospital admissions for CVDs, according to the block diagram of the methodology (Fig. 5). Feature Importance shows how the most characteristic variables for CVDs, Fig. 12, are atmospheric pressure with a percentage of 37%, minimum temperature with 19%, and carbon monoxide (CO) with a percentage of 18%. These three variables represent about 74% of the model’s predictability, with atmospheric pressure being predominant.

The SHAP model integrated and analyzed cumulative datasets derived from the previous models in this study. This data was visualized using bee swarm plot diagrams, as illustrated in Fig. 13. While the features are ranked by their contributed weight to the prediction, with consideration given to the range and value (i.e., $\pm$) in the model.

The SHAP method provides insights into both the global and local contributions of each variable, with the results presented through an importance ranking visualized through a bee swarm scatter plot. Within this representation, the horizontal axis signifies the SHAP value, while the color of the point indicates the intensity of the values (blue for low values, red for high values) of each variable.

The SHAP model identified the variables i with the most pronounced impact on hospital admissions. These variables are arranged in descending order of importance as follows: Atmospheric Pressure (P_atm), Minimum Temperature (Tmin), and Carbon Monoxide (CO) concentration. In Fig. 13, the correlation between low values of P_atm and Tmin and their correspondence to elevated SHAP values (indicating significance) is readily discernible. This graphical representation underscores their substantial role in enhancing the effectiveness of the simulation model. Conversely, the behavior of CO exhibits the opposite pattern, increasing in importance for higher values.

The remaining variables were found to exhibit less significance. These include: Dew Point (Tdewp), Relative Humidity (RH), Ozone (O${}_3$), and Mean Temperature (Tmean). On the other hand, Tmax provides a nearly negligible contribution, regardless of its values. This comprehensive analysis allows for a nuanced understanding of the varying impacts of these variables on hospital admissions.

The results of this study provide valuable insights into the relationship between environmental factors and cardiovascular health. Our findings suggest that atmospheric pressure, minimum temperature, and carbon monoxide are the most important factors in predicting the number of hospital admissions for CVD. These results are consistent with previous studies that have identified these factors as risk factors for CVD [6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 20]. Previous studies have consistently shown that extreme temperatures increase mortality from CVD [9, 10, 11, 12, 13]. This association between temperature and CVD mortality is observed for both low and high temperatures [22]. Chronic exposure to cold or heat can impair cardiovascular function, consequently increasing the risk of heart attack, arrhythmias, thromboembolic diseases, and heat-induced sepsis [23]. Furthermore, changes in ambient temperature have been shown to contribute to cardiovascular mortality by causing increased blood pressure, blood viscosity, and heart rate [23]. Seasonal variations in CVD pose a significant health challenge [24], particularly for populations residing in milder climates, who may be less adapted or prepared for extreme climate changes. A majority of studies have reported a winter increase in cardiovascular anomalies and cardiac death in the northern hemisphere, where temperatures are exceedingly cold [12, 25]. In light of this evidence, it is evident that temperature plays a crucial role in the prevalence of CVDs and should be considered a vital parameter when developing predictive models for hospital admissions related to these conditions.

Our results not only emphasize the importance of temperature but also highlight the significance of atmospheric pressure and carbon monoxide concentration in predicting hospital admissions for CVDs. These findings are consistent with existing literature. A connection between small changes in low levels of ambient carbon monoxide concentrations and cardiovascular hospitalizations has been demonstrated [42]. Exposure to carbon monoxide at concentrations found in heavy tobacco smokers or individuals with significant occupational exposure has been shown to play a role in the pathogenesis of CVD [43]. Short-term exposure to ambient carbon monoxide has been associated with an increased risk of CVD hospitalizations [42], as well as increased rates of cardiovascular and coronary heart diseases in major Chinese cities [44]. Ambient carbon monoxide levels have also been positively correlated with hospital admissions for CVDs [45]. Moreover, carbon monoxide exposure has been shown to have detrimental effects on exercise performance in subjects with coronary artery disease, further demonstrating its impact on myocardial ischemia [46]. In addition to carbon monoxide, high atmospheric pressures have been associated with an increased number of cardiac arrest admissions [47]. This evidence suggests that both atmospheric pressure and carbon monoxide concentration are crucial factors to consider when developing predictive models for hospital admissions related to CVDs. Given the growing body of research, it is crucial to incorporate temperature, atmospheric pressure, and carbon monoxide concentration into predictive models for hospital admissions due to CVDs, as these factors play a significant role in the onset and exacerbation of cardiovascular conditions.

Recent research has provided significant contributions to understanding and addressing the negative impacts of climate change on CVD. One such study [48] proposes a healthcare service recommendation framework based on an embedded user profile model. This personalized approach could be extended to the domain of CVD, enabling targeted healthcare service suggestions based on meteorological conditions and air quality. An innovative approach, presented in [49] and [50], utilizes deep learning through convolutional neural networks or ResNet-based classification models to identify the presence of diseases, such as COVID-19, from radiological images. This methodology could also be adapted for early diagnosis of CVD, correlating specific radiological signs.

Furthermore, Yu and colleagues [51] have devised an approach that aims to acquire insights into the causality of diseases and foster a more comprehensive understanding of the interconnections between symptoms and diseases. This framework shows potential in revealing causal connections between meteorological factors, air quality, and CVDs. Collectively, these studies collectively contribute to enhancing our understanding of the complex relationships between environmental variables and cardiovascular health.

Future studies should prioritize establishing alert thresholds based on these parameters, with additional efforts dedicated to expanding the scope of the model to encompass a wider array of diseases and other geographical regions. To enhance the model’s predictive prowess, validation procedures and potential refinements will be investigated. This could involve integrating data sourced from multiple hospitals situated across different geographical areas and climates. The ultimate goal is to gain access to the comprehensive regional database of Puglia hospitals to validate the model’s generalizability. This step will serve to reinforce its reliability and broaden its applicability.

In conclusion, there is a strong relationship between daily Emergency Room admissions and the variation of some weather and air quality parameters. Future studies should be directed toward pinpointing and defining alert thresholds rooted in these parameters.