## 2 Materials and Methods  
### Data Preparation Description  
The data used were those of a previously conducted study [16] on air quality and atmospheric composition, acquired by Copernicus Atmosphere Monitoring Service (CAMS) ([[https://ads.atmosphere.copernicus.eu/#](https://ads.atmosphere.copernicus.eu/#)]([https://ads.atmosphere.copernicus.eu/#](https://ads.atmosphere.copernicus.eu/#))!/search?text=&type=dataset&keywords=((%20%22 Product%20 type%20 Reanalysis%22%20)%20 AND%20(%20%22 Variable%20 domain:%20 Atmosphere%20(composition)%22%20)%20 AND%20(%20%2022 Spatial%20 coverage:%20 Europe%22%20)%20 AND%20(%20%22 Temporal%20 coverage:%20 Past%22 %20)) accessed on 1 November 2023). In particular, the pollutant data provided by CAMS are derived from the median ensemble of nine numerical air quality models [20; 21; 22; 23; 24; 25; 26; 27]. In this study, which examined the relationship between air pollution and mortality for endocrine nutritional and metabolic diseases, specific pollutant concentrations such as CO, NO, NO\({}_{2}\), O\({}_{3}\), PM10, PM2.5, and SO\({}_{2}\) were analyzed using the Pythonlibrary Geopandas ([[https://geopandas.org/en/stable/](https://geopandas.org/en/stable/)]([https://geopandas.org/en/stable/](https://geopandas.org/en/stable/))) (version 0.14.1, 1 November 2023). To enhance data consistency, additional variables related to pollution and socio-economic factors were included. Figure 1 presents an outline of the data preparation workflow. The dataset thus constructed consisted of 107 Italian provinces, each described by 28 environmental and socio-economic variables. In this study, we calculated the standardized mortality ratio (SMR), which represents the ratio of observed deaths in a given area to the number of the expected ones. The expected mortality count was estimated through the age mortality rates from a reference population (in this case, the national population) to the mean annual population, stratified by age groups, within each province. Figure 2 illustrates the spatial distribution of the SMR for endocrine, nutritional, and metabolic diseases across Italian provinces in 2019. Figure 3 provides a visual overview of the analytical process carried out in this study. The first step of this work was to create a well-structured dataset of the Italian provinces with multiple types of data: the air pollution data were collected from the CAMS platform, while the socio-economic indices, such as life quality distribution and density of urban green areas, and the mortality data related to endocrine, nutritional, and metabolic diseases were obtained from the database of ISTAT. During the preliminary phase of the analysis, the multicollinearity of the dataset's variables was assessed by calculating Variation Inflation Factor (VIF). All variables with a VIF greater than 10 were iteratively removed from the dataset. Figure 1: Dataset preparation workflow. We identified groups of neighboring provinces (clusters) that were used separately in the training and testing phases to prevent data leakage between spatially correlated observations. These clusters of provinces were determined using a combination of methods. The analysis involved utilizing a hierarchical clustering approach on the matrix containing Euclidean distances calculated between provinces. Additionally, the semivariogram graph of the standardized mortality ratios (SMRs) was employed to determine the optimal stopping point for the hierarchical clustering process. The spatial range derived from the semivariogram, which indicates the distance at which data points are considered to no longer exhibit correlation, was visually estimated. This was performed by identifying the point where the plot leveled off (the plateau), indicating that the variance no longer increased significantly with distance. To minimize spatial bias during the assessment of the regression performance, we implemented a leave-one-cluster-out cross-validation approach.