Explaining deep learning models for ozone pollution prediction via embedded feature selection

Available online 21 March 2024 1568-4946/© 2024 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Explaining deep learning models for ozone pollution prediction via embedded feature selection M.J. Jiménez-Navarroa,∗, M. Martínez-Ballesterosa, F. Martínez-Álvarezb, G. Asencio-Cortésb aDepartment of Computer Languages and Systems, University of Seville, ES-41012, Seville, Spain bData Science and Big Data Lab, Pablo de Olavide University, ES-41013, Seville, Spain A R T I C L E I N F O Keywords: Pollution Time series forecasting Ozone Feature selection Deep learning XAI A B S T R A C T Ambient air pollution is a pervasive global issue that poses significant health risks. Among pollutants, ozone (O3) is responsible for an estimated 1 to 1.2 million premature deaths yearly. Furthermore, O3adversely affects climate warming, crop productivity, and more. Its formation occurs when nitrogen oxides and volatile organic compounds react with short-wavelength solar radiation. Consequently, urban areas with high traffic volume and elevated temperatures are particularly prone to elevated O3levels, which pose a significant health risk to their inhabitants. In response to this problem, many countries have developed web and mobile applications that provide real-time air pollution information using sensor data. However, while these applications offer valuable insight into current pollution levels, predicting future pollutant behavior is crucial for effective planning and mitigation strategies. Therefore, our main objectives are to develop accurate and efficient prediction models and identify the key factors that influence O3levels. We adopt a time series forecasting approach to address these objectives, which allows us to analyze and predict future O3behavior. Additionally, we tackle the feature selection problem to identify the most relevant features and periods that contribute to prediction accuracy by introducing a novel method called the Time Selection Layer in Deep Learning models, which significantly improves model performance, reduces complexity, and enhances interpretability. Our study focuses on data collected from five representative areas in Seville, Cordova, and Jaen provinces in Spain, using multiple sensors to capture comprehensive pollution data. We compare the performance of three models: Lasso, Decision Tree, and Deep Learning with and without incorporating the Time Selection Layer. Our results demonstrate that including the Time Selection Layer significantly enhances the effectiveness and interpretability of Deep Learning models, achieving an average effectiveness improvement of 9% across all monitored areas. 1. Introduction Approximately 4.2 million premature deaths worldwide are related to air pollution [1]. Among these, 1 to 1.2 million respiratory deaths in adults per year are attributable to tropospheric ozone (O3) exposures [2]. As a greenhouse gas, ozone plays a critical role in climate warming and is a significant indicator of photochemical pollution due to its strong oxidizing properties [3]. Furthermore, ozone adversely affects crop yields and corrodes many metallic and non-metallic materials [4]. The formation of ozone involves primary pollutants emitted, such as nitrogen oxides (NO𝑥= NO +NO2) and volatile organic compounds (VOCs) in the presence of short-wavelength solar radiation. Urban areas are hotspots for human exposure to air pollution, mainly from road traffic [5]. In these areas, ozone precursors are emitted from many sources, with traffic emissions being the major contributor. Urban ∗Corresponding author. E-mail addresses: [email protected] (M.J. Jiménez-Navarro), [email protected] (M. Martínez-Ballesteros), [email protected] (F. Martínez-Álvarez), [email protected] (G. Asencio-Cortés). traffic is estimated to contribute to ambient nitrogen dioxide (NO2) concentrations ranging from 3% to 56% [6]. This circumstance leads urban agglomerations to experience episodes of high ozone concentration on sunny days with high temperatures and traffic density, particularly in cities in the Mediterranean basin. Ozone concentrations reach a maximum peak in summer and a minimum peak in winter [7] while tending to be higher on weekends compared to weekdays at urban sites (aka ozone weekend effect) [8]. Typically, ozone also exhibits a diurnal profile of concentrations, with high noontime concentrations due to higher solar radiation and temperature. Accurate air quality forecasts in urban environments have become crucial tools for early warning, protecting public health, and developing mitigation strategies. Therefore, predicting the temporal evolution of ozone concentrations in urban locations [9–11] is a priority to ensure quality of life, provide https://doi.org/10.1016/j.asoc.2024.111504 Received 31 July 2023; Received in revised form 1 February 2024; Accepted 17 February 2024 2 M.J. Jiménez-Navarro et al. accurate information to the population of affected areas and alert them when there are exceptionally high levels [12]. In addition, understanding the behavior of the pollutants and the influent variables has been extensively studied [13–15]. Several models have been proposed for time series forecasting [16], achieving remarkable results in various problems. However, there is a growing demand for accurate predictions and computationally efficient and interpretable solutions. Computation resources can be reduced using hardware- and/or software-centric techniques [17], promoting the development of green AI [18]. Various approaches exist to improve model interpretability in artificial intelligence, such as global/local model-agnostic or example-based methods. In this context, feature selection [19] is a unique technique that offers benefits in reducing computational costs and improving interpretability. This technique aims to identify a subset of relevant features that are sufficient or more informative than the original set of features. This technique reduces computational resources regarding memory and execution time, and the solution becomes better understood by focusing on the relevant identified features. Furthermore, in time series forecasting problems, providing information about the relevant features along the temporal dimension is essential. In this work, we study a dataset collected from six areas of Seville and Jaen provinces, Spain. Our study focuses on models incorporating a global interpretability method, which embeds the feature selection mechanism within the model design instead of other feature selection approaches. Thus, we can determine future behavior and relevant features once the model is obtained. Examples of models that incorporate feature selection include tree-based methods [20] (e.g., Decision Tree or CART) and regularization-based methods [21] (e.g., Lasso or L1). Neural networks may be a suitable approach for time series forecasting [22]. However, two main issues remain open: computational cost and interpretability. The Temporal Selection Layer (TSL) proposed in Jimenez-Navarro et al. [23,24] addresses both of these issues by transforming a neural network into a model with embedded feature selection which attempts to improve the efficacy and interpretability while reducing computational costs1. Therefore, this study presents the following contributions: •An improvement of the originally proposed TSL by adapting regularization for multi-step forecasting problems. •Conducted a real-world scenario comparison of Decision Tree (DT), L1, Feed-Forward (FF), and Long Short-Term Memory (LSTM) neural networks. The selected models embed a feature selection process and are known for their simplicity and efficiency. This study is outlined as follows. Relevant and previous work are reviewed in Section 2. In Section 3, the proposed improvement is introduced along with the experimental workflow. The main results of the performance of Section 4and their comparison with the rest of the algorithms are discussed. Finally, Section 5summarizes the main conclusions drawn from this study. 2. Related work Air pollution is a well-known problem with several proposed solutions, but not all include a feature selection process. In this section, we analyze several approaches to feature selection for neural networks to justify their effectiveness in a wide range of problems. The approaches have been categorized into three types [25] of feature selection approaches: filter, wrapper, and embedded. The filter approach relies on a heuristic to assess the importance of each feature and then establishes a threshold to remove irrelevant 1All the code needed to reproduce the results has been included in the following repository: https://github.com/manjimnav/TSLayer-Ozone, along with detailed instructions. ones. Although this approach is often efficient, it does not guarantee the selection of the most relevant features. It may not necessarily correspond to the importance within the assigned model, which limits the reliability in interpreting how the model uses these selected features. The wrapper approach employs the model itself as a heuristic, retraining it with different feature sets. This approach can yield optimal feature selections for a specific model depending on the process and may enhance interpretability by aligning the features more closely to the model. However, it comes with the drawback of potentially being computationally intensive, making it challenging to scale for a larger number of features and complex models. The embedded approach usually offers the best approach, combining the strengths of the filter and wrapper methods by efficiently selecting features that are both optimal for the model and interpretable. Furthermore, this feature selection approach usually does not require hyperparameter tuning in contrast to other approaches, simplifying the training process and reducing computational costs. However, implementing these methods involves developing an algorithm that incorporates the feature selection procedure into its design. 2.1. Filter In their study, Khemphila et al. [26] utilize information gain scoring to rank features in a heart disease dataset. A threshold-based feature selection is applied, and a neural network is trained with the selected data. Unlike their approach, our work focuses on embedded methods, wherein the selected features are directly related to the trained model. Lui et al. [27] construct an ensemble of neural networks grouped into threes. Each group employs different selection methods, including PCA, RankSum, and t-test. The group outputs are averaged, and majority voting is used for final predictions. This method significantly impacts efficacy but performs poorly due to indiscriminate feature selection. Yan [28] adopts an AdaBoost model to identify relevant features and trains a neural network using these features. While effective for stock index prediction, this method shares a typical challenge with filter methods: the feature selection process remains separate from the model, potentially leading to suboptimal feature choices for the neural network. Niu et al. [29] employ the non-linear feature selection algorithm RReliefF. The relevant features this algorithm determines are then used as input for neural networks, resulting in considerable performance improvements across all scenarios. 2.2. Wrapper Verikas et al. [30] introduce a feature ranking method based on cross-entropy error constraints. They recursively remove features from least to most important during neural network training. Our paper focuses on efficacy and performance, aiming to minimize computational resources during feature selection. Kabir et al. [31] propose an incremental feature selection and model-building method. They use feature correlations to divide features into groups and incrementally add them during training if they improve error rates. Tong et al. [32] use a genetic algorithm in which individuals represent models, leading to high memory consumption for numerous populations. The computational cost can become prohibitive with large models or datasets, making it unsuitable for our paper’s scope. Alshammari et al. [33] utilize a neural network and the Mayfly metaheuristic algorithm to optimize the number of relevant features. In contrast, our method proposes an efficient and effective embedded feature selection approach, avoiding the need for a significant increase in computational resources. 3 M.J. Jiménez-Navarro et al. Fig. 1. Summary of experimentation workflow. The data for each station 𝑆𝑖is transformed, used for the training, and finally evaluated. During the preprocessing, the data is divided in different folds returning the test 𝑋𝑡𝑒𝑠𝑡 𝑖,𝑦 , corresponding to the year 𝑦, and training dataset 𝑋𝑡𝑟𝑎𝑖𝑛 𝑖,𝑦 with the remaining years for station 𝑖. Then, the missing data is inputted, the data is standardized using the 𝑋𝑡𝑟𝑎𝑖𝑛 𝑖, and the lagged information is built. Furthermore, the training process is repeated for each preprocessed data 𝑋𝑡𝑟𝑎𝑖𝑛 𝑖,𝑦 , selected model 𝑀𝑗, and sampled hyperparameters 𝐻𝑗 𝑘from the Bayesian optimization for 𝐾iterations. Finally, the model is evaluated for the 𝑋𝑡𝑒𝑠𝑡 𝑖,𝑦 to obtain the results. 2.3. Embedded Cancela et al. [34] introduce the E2E-FS technique for integrating feature selection within neural networks. This method employs an additional loss function strategy to effectively filter out irrelevant features. In contrast, our method automatically determines the number of features through an approximation of the Heaviside function to ensure differentiability without modifying the loss function. Da Costa et al. [35] propose a feature selection method that enhances the stability of feature rankings applied to neural networks. However, our focus is automatically selecting the best subset of features to improve efficacy or achieve comparable results with fewer features. Our approach does not rely on a two-step process involving input layer weight analysis and subsequent feature ranking. Yuan et al. [36] present an embedded feature selection method for neural networks using the point-centered convolutional neural network. While effective in identifying moldy peanuts using hyperspectral images, their approach assigns continuous weights, making determining the relevance of specific features challenging. Zhang et al. [37] propose incorporating feature selection within a neural network framework using the Group Lasso penalty. Although effective for Feed-Forward layers, this approach limits its use in other models. On the contrary, our research offers a versatile and generalpurpose approach that applies to various neural network architectures. 3. Methodology This section details the methodology proposed to evaluate the different approaches. Section 3.1 shows the proposed models employed in the comparison. Section 3.2 describes TSL and its background. Section 3.3 shows the improvement applied to TSL in this work. Section 3.4 shows the data used during the experimentation and its properties. Section 3.5 provides an overview of how the experimentation workflow is performed. Fig. 1 summarizes the workflow followed during the experimentation to obtain the results. 3.1. Models For comparative purposes, we established the following criteria for selecting the models: •Models should have efficient inference time. •Models should be interpretable, meaning they should indicate relevant features. •Simplicity is preferred, as discussed in Section 4, where we observed that simpler models yield remarkable results for this problem. Firstly, we propose two simple and very efficient models that also provide an embedded feature selection. Decision tree (DT) and Lasso (L1) models were chosen due to their well-known simplicity and high efficiency. In this case, both models use a direct approach for multi-step time series forecasting. Feed-Forward neural networks (FF) and Long Short-Term Memory (LSTM) [38] were used by employing the direct approach to produce multi-step forecasting. Note that neural networks are not usually considered simple in other problems. However, in our experience, neural networks do not need to be particularly deep and wide to solve time series forecasting problems [39,40]. Therefore, we decided to employ this model type as it accomplishes the specified criterion in this context. 3.2. Time selection layer The FF and LSTM neural networks are employed in two ways: including and not including the TSL. For instance, when the TSL is used, the model is named TLSTM in the case of LSTM and TFF for FF. This layer is formally defined as the element-wise multiplication of the binarized weights and the input tensor, defined in the following equation: 𝑇 𝑆𝐿(𝑋𝑀𝑥𝐷) = 𝐻(𝑊𝑀𝑥𝐷)◦𝑋𝑀𝑥𝐷 (1) where ‘‘◦’’ represents the Hadamard product, 𝑊represents the window size of the time series, and 𝐷is the dimension size. Moreover, an approximation to the Heaviside function (𝐻) serves as a gate that indicates if the feature is selected. Using backpropagation, TSL is trained as an additional layer at the top of any neural network. When a weight in the feature selection mask is set to zero, it eliminates the contribution of the corresponding feature in the forward step. Weight regularization is essential in this layer, applied as a Lasso regression to the feature selection mask to prevent undesirable outcomes such as not selecting any feature. This regularization is calculated as the sum of a regularization factor multiplied by the number of selected features. This regularization is then included in the loss of the model. 3.3. Regularization improvement in time selection layer We present an extension of our previous work [23] with an additional improvement. Our proposal introduces TSL, an extra layer incorporating embedded feature selection into a neural network. TSL layer was previously used with remarkable results in time series forecasting. However, it was not utilized for multi-step forecasting. We hypothesize that forecasting larger periods requires more lagged information to reproduce future behavior effectively. Therefore, the need for additional lagged features becomes apparent as the horizon increases. 4 M.J. Jiménez-Navarro et al. Table 1 Table with datasets relevant properties employed in all areas. Property Aljarafe Asomadilla Bermejales Torneo Ronda del Valle Location Seville Cordova Seville Seville Jaen Area type Suburban Suburban Urban Urban Urban Emission source Background Background Background Traffic Background Features TMP, WD, WP, CO, NO2,𝑃 𝑀10 and O3 Sample frequency Hourly Sample period 2006–2023 Ozone unit μg∕m3 To accommodate this, we modified the regularization term in TSL using the following formula: ∑𝑁 𝑖𝜎𝐹𝑖 𝐻 where 𝑁represents the number of features, 𝜎∈Rrepresents the regularization term, 𝐹𝑖∈ 0,1represents whether the feature is selected (1) or not (0), and Hrepresents the number of future steps to predict. 3.4. Data Hourly O3multivariate time series were obtained annually from 2006 to 2023. Data were collected from four monitoring stations in the metropolitan area of Seville and one in Jaen. These monitoring sites were classified according to the type of area (U-Urban, S-Suburban) and the predominant emission source (B-Background, T-Traffic). The sites analyzed were Aljarafe (S-B), Bermejales (U-B), Torneo (U-T) in Seville, Asomadilla (S-B) in Cordova, and Ronda del Valle (U-B) in Jaen. This selection allowed for considering various meteorological conditions and local factors that influence the behavior of O3. Meteorological conditions included average temperature (TMP), wind direction (WD), and wind speed (WS). The local factors considered were carbon monoxide (CO), nitrogen dioxide (NO2), and particles of 10 μm or less (𝑃 𝑀10). The O3measurements were obtained using the reference monitoring method [41]. Table 1 summarizes the properties of the datasets employed. 3.5. Experimentation process Due to the differences in the behavior of each zone, we decided to model each one individually. The same methodology was applied in each zone, divided into preprocessing, training, and evaluation. 3.5.1. Preprocessing To prepare the data for the experimentation process, four preprocessing steps are applied to each time series: missing data imputation, splitting, standardization, and windowing. Sensor errors are common and pose challenges in computer science and telecommunications. In this case, the sensors in the monitoring stations may not provide data for specific short periods, which can complicate the experimentation process for some models. A data imputation technique is applied to address this issue. In this case, we employed the MissForest [42] technique, which uses a random forest to fill the data, considering all the non-missing features from the current moment during five iterations. A rigorous data imputation comparison is beyond the scope of this work. For that reason, we do not consider additional methods. Following data imputation, a cross-validation approach is employed to assess the generalization performance of each model. In each crossvalidation iteration, an entire year is used as the test set, the previous year is used as the validation set, and the remaining years are used as the training set. This means that, in each fold, the training size would consist of 16 years (≈140160 instances or ≈90%) while the validation and test size would be one year each one (≈8760 instances or ≈5% each one). Table 2 Hyperparameters ranges for each model employed during the optimization. Model Hyperparameter Range DT Maximum depth [2, 32] L1 Regularization [0, 1] FF & LSTM Layers [1, 3] Units [4, 64] Dropout [0, 0.5] TSL Regularization [0.01, 1e−7] During cross-validation, each data split is standardized by transforming the data to have a zero mean and a standard deviation of one. The mean and standard deviation of the training set are calculated and used to standardize the validation and test sets. Moreover, when dealing with time series problems, it is necessary to incorporate lagged temporal information as input for the models. This means transforming the original time series 𝑋𝑁,𝐷 with 𝑁samples and 𝐷features into 𝑋𝑁,𝐷,𝑊 , where 𝑊represents the number of lagged information (window size) included for each sample and feature. Finally, since the goal is to predict the future, the matrix 𝑌𝑁,𝐻 is constructed for each sample, providing 𝐻future Ozone (O3) values that constitute the prediction. 3.5.2. Training Each selected model has its own set of hyperparameters that greatly influence its effectiveness. A search algorithm is used to determine the best set of hyperparameters, which aims to find the optimal configuration across all the tested folds. The chosen search algorithm is Bayesian search [43], which has shown outstanding performance in recent studies [44,45]. In each iteration, the algorithm generates a set of hyperparameters. The model is then parameterized with these hyperparameters, trained on each training set, and subsequently tested on the corresponding test set within the cross-validation process. During training, neural networks utilize the validation set to stop the training process when no further improvement is observed for the validation efficacy in at least ten epochs. As the optimization algorithm requires establishing a range for the hyperparameters, Table 2 summarizes it for each model. Furthermore, the window size (𝑊) used during preprocessing (see Section 3.5.1) is also optimized to determine the best amount of previous information together with the learning rate and the batch size in the cases of neural networks. In all cases, the hyperparameters have been selected to cover a range which should generate simple models. Finally, a list of ranges is provided in Table 3 with all the hyperparameters relevant to the experimentation. Note that the hyperparameters whose range has a single value are not optimized. Finally, as the Bayesian Search itself incorporates hyperparameters, the default values have been chosen for use in general optimization problems. Specifically, the Upper Confidence Bound acquisition function has been used with a 𝜆value of 2.576, a choice that has been consistently shown to strike a good balance between exploration and exploitation. 5 M.J. Jiménez-Navarro et al. Fig. 2. Results obtained by the Bayesian test applied to TLSTM (left) and L1 (right). Table 3 Other relevant hyperparameters. Hyperparameter Range Window size [24, 72] Horizon 24 Shift 24 Learning rate [0.0001, 0.01] Batch size [16, 128] Epochs 100 Early stopping patience 10 𝜆2.576 3.5.3. Evaluation The primary metric chosen was the Root Mean Squared Error (RMSE) to evaluate the effectiveness of the model as this metric is widely used in air pollution problems. The final step in the workflow is to determine the optimal set of hyperparameters and models. This determination is based on the average efficacy in all years of cross-validation. 4. Results and discussion This section presents the results and discusses the implications and hypotheses on the metrics obtained. Section 4.1 provides an overview of the results for all datasets using grand averages. Section 4.2 presents the RMSE distribution among the different experiments. Section 4.3 provides the best hyperparameters obtained from evaluating our methodology. Section 4.4 analyzes the results on a yearly basis. Section 4.5 shows the error obtained daily and monthly. Section 4.6 presents the best and worst predictions achieved by the best model. Section 4.7 illustrates which features were selected by the best model. Finally, Section 4.8 contains the distribution of the number of features obtained in each experiment. 4.1. General analysis Table 4 displays the average cross-validated metrics obtained after applying our methodology. Although datasets are sampled from sensors with comparable conditions, the errors of the different models vary considerably in some cases. It appears that the behavior of Aljarafe and Asomadilla is quite similar in terms of RMSE, which does not seem to be related to geographic proximity, as other areas studied are closer. Table 4 The average RMSE results obtained for the best models. We conducted a statistical comparison between the L1 and TLSTM models and marked the best results with an asterisk to indicate that their differences are significant. DT L1 FF LSTM TFF TLSTM Aljarafe 15.9 15.2 15.8 15.1 15.3 14.9* Asomadilla 15.9 15.2 15.8 15.2 15.2 14.9* Bermejales 17.9 17.0 17.8 17.2 17.2 16.9 Ronda del Valle 18.1 18.0 18.1 17.7 17.4 17.0* Torneo 14.7 13.9 14.5 14.1 13.9 13.9 Average 16.5 15.9 16.4 15.8 15.8 15.5 On the other hand, the Bermejales and Ronda del Valle datasets appear to be the most difficult to predict, with an RMSE difference of 2 to 3 compared to Aljarafe or Asomadilla. Interestingly, the Torneo dataset achieved a lower RMSE, suggesting that it might be the easiest area to predict. We can generally observe that more straightforward approaches significantly improve the results in almost all datasets. This suggests that the problem is relatively simple and that using more complex approaches hampers generalization. Overall, our TLSTM consistently improves the simple LSTM by simply embedding the feature selection. The TLSTM achieved the best RMSE for all datasets, followed by the TFF, LSTM, and the L1 models with only a difference of 0.3–0.4 in terms of RMSE. On the other hand, DT and FF generally provide the worst results. There is a clear improvement in adding the TSL to the neural networks compared with their original versions (FF and LSTM), especially in the FF neural network. We hypothesize that TSL adds flexibility to the neural network, enabling it to adjust to simpler models. Thanks to our proposal, the neural network can adapt its complexity to match the behavior of the problem. For example, in simpler problems that require fewer features, TSL may eliminate a few, whereas in more complex scenarios, TSL may retain a greater number of features. For example, this flexibility cannot be achieved with the L1 model, as it is only suitable for linear problems. Finally, we compared the results between our best TSL combination and the model with embedded feature selection. Specifically, the L1 and TLSTM models have been compared using a Bayesian test [46] based on the findings in Section 4.4. The best result with a significance level above 75% is indicated with an asterisk. Refer to Fig. 2 for detailed comparison results. 6 M.J. Jiménez-Navarro et al. Fig. 3. RMSE distribution of the different experiments. 4.2. Error distribution analysis Fig. 3 illustrates the RMSE distribution of every experiment performed grouped by model and dataset. This plot may allow us to observe the difficulty in finding optimal hyperparameters. In general, the neural networks are more robust to the selected hyperparameters than the DT and L1 models employed within the optimization space, even when the search space in the neural network is considerably greater. This robustness gives us greater confidence in the model configuration even if it was selected without any hyperparameter optimization process. The FF and TFF models exhibit a more significant difference in all datasets. However, the improvement is less significant between LSTM and TLSTM, which could be explained as the LSTM includes some selection in their forget gates. However, including TSL appears to be beneficial for LSTM on average. Finally, there seems not to be a relationship between the results obtained by our approach and the characteristics of the datasets shown in Table 1. 4.3. Hyperparameters analysis This section presents and discusses the best hyperparameters found by our methodology. The goal is to identify existing patterns among the different datasets. Table 5 summarizes the best hyperparameters for each dataset and model. Note that the hyperparameters are in JSON format using the ‘‘key:value’’ format. In general, there seems to be no recognizable pattern that relates the optimal hyperparameters and the results in all datasets. Regarding the DT model, it is clear that the best depth for the datasets is around five or six, which is considerably less than the maximum available. This supports the idea that the problem requires simpler models, as increasing the depth would increase the probability of overfitting. The L1 model shows more diversity, with the predominant low values ranging from 0.0001 to 0.02. For Aljarafe, Asomadilla, Bermejales, and Torneo, a value of 0.0001 was the best, indicating a considerably low penalization for feature selection. However, the regularization is significantly greater for Ronda del Valle, suggesting that the model should tends to select less features. The FF and TFF models have consistently opted for the same number of layers, typically one or two, without utilizing the maximum available (three) layers in any dataset. As for the number of units, the FF model typically ranged from 12 to 59, whereas the TFF model required considerably fewer units, ranging from 7 to 18. The batch size and learning rate do not show a consistent pattern among the datasets. The FF model seems to have increased its complexity in order to be able to process all features, decreasing its efficacy as a consequence. The dropout rate remained generally low, less than 0.3, except for the Torneo and Bermejales datasets in FF. The LSTM and TLSTM models typically have a similar number of layers across datasets, with differences only seen in Aljarafe, where TLSTM requires fewer layers. However, the situation changes when considering the number of units compared to the FF and TFF models. TLSTM generally requires more units, ranging from 36 to 48, whereas LSTM ranges from 11 to 45. Additionally, TLSTM demands significantly lower batch sizes and learning rates compared to LSTM. Dropout rates in LSTM are consistent across datasets, except for Bermejales. In general, time series forecasting problems typically do not require neural networks with many layers. Thus, it is reasonable that hyperparameter optimization identified a low number of layers and units as optimal. In general, for TSL applied to all neural networks, the regularization term of the layer converged to a close range for all datasets, between 0.002 and 0.01. Consequently, we can assume that a value within this range could be a good initial value for datasets with similar circumstances. Typically, a large window size in all datasets is necessary, often surpassing twice the forecasting horizon. Not all time steps may be used in the DT, L1, TFF, and TLSTM models, as they include a feature selection process. This will be discussed in Section 4.7. 4.4. Yearly analysis This section presents the results for each year obtained after applying cross-validation to the best model. Note that the DT and FF models have been removed for visual clarity, as they seem to be worse in all datasets than the L1, TFF, and TLSTM models. Table 6 displays the average RMSE for each model in all datasets and years evaluated. The results, as discussed in Section 4.1, indicate that the three most promising models which embeds some feature selection are the L1, TFF, 7 M.J. Jiménez-Navarro et al. Table 5 Best hyperparameters found by the experimentation process. Dataset Model Window Hyperparameters Aljarafe DT 40 # Max depth: 7 FF 34 # Layers: 1, # Units: 19, Batch size: 48, Learning rate: 0.0068, Dropout: 0.0650 L1 38 Regularization: 0.0001 LSTM 67 # Layers: 2, # Units: 12, Batch size: 90, Learning rate: 0.0059, Dropout: 0.2574 TFF 51 # Layers: 1, # Units: 16, Batch size: 95, Learning rate: 0.0008, Dropout: 0.2798, Regularization: 0.0096 TLSTM 65 # Layers: 1, # Units: 48, Batch size: 18, Learning rate: 0.0080, Dropout: 0.4003, Regularization: 0.0039 Asomadilla DT 45 # Max depth: 5 FF 67 # Layers: 2, # Units: 12, Batch size: 90, Learning rate: 0.0059, Dropout: 0.2574 L1 38 Regularization: 0.0001 LSTM 67 # Layers: 2, # Units: 12, Batch size: 90, Learning rate: 0.0059, Dropout: 0.2574 TFF 64 # Layers: 1, # Units: 18, Batch size: 44, Learning rate: 0.0059, Dropout: 0.3719, Regularization: 0.0097 TLSTM 71 # Layers: 2, # Units: 48, Batch size: 51, Learning rate: 0.0002, Dropout: 0.3432, Regularization: 0.0075 Bermejales DT 25 # Max depth: 4 FF 68 # Layers: 1, # Units: 45, Batch size: 93, Learning rate: 0.0014, Dropout: 0.4986 L1 38 Regularization: 0.0001 LSTM 68 # Layers: 1, # Units: 45, Batch size: 93, Learning rate: 0.0014, Dropout: 0.4987 TFF 47 # Layers: 2, # Units: 7, Batch size: 30, Learning rate: 0.0022, Dropout: 0.0096, Regularization: 0.0027 TLSTM 65 # Layers: 1, # Units: 36, Batch size: 29, Learning rate: 0.0004, Dropout: 0.0099, Regularization: 0.0024 Rondadelvalle DT 28 # Max depth: 6 FF 32 # Layers: 1, # Units: 56, Batch size: 114, Learning rate: 0.0005, Dropout: 0.4473 L1 56 Regularization: 0.0273 LSTM 26 # Layers: 2, # Units: 11, Batch size: 112, Learning rate: 0.0015, Dropout: 0.3736 TFF 51 # Layers: 1, # Units: 16, Batch size: 95, Learning rate: 0.0008, Dropout: 0.2799, Regularization: 0.0097 TLSTM 71 # Layers: 2, # Units: 48, Batch size: 51, Learning rate: 0.0003, Dropout: 0.3433, Regularization: 0.0075 Torneo DT 45 # Max depth: 5 FF 24 # Layers: 2, # Units: 59, Batch size: 115, Learning rate: 0.0063, Dropout: 0.1786 L1 38 Regularization: 0.0001 LSTM 26 # Layers: 2, # Units: 11, Batch size: 112, Learning rate: 0.0015, Dropout: 0.3736 TFF 51 # Layers: 1, # Units: 16, Batch size: 95, Learning rate: 0.0008, Dropout: 0.2799, Regularization: 0.0097 TLSTM 71 # Layers: 2, # Units: 48, Batch size: 51, Learning rate: 0.0003, Dropout: 0.3433, Regularization: 0.0075 Table 6 Best results obtained for each model and year in terms of RMSE (μg∕m3). Aljarafe Asomadilla Bermejales Ronda del valle Torneo L1 TFF TLSTM L1 TFF TLSTM L1 TFF TLSTM L1 TFF TLSTM L1 TFF TLSTM 2006 16.8 17.1 16.2 17.4 17.7 17.1 18.0 18.9 18.7 20.3 20.3 19.8 13.1 13.1 12.9 2007 15.5 16.4 15.6 16.4 16.3 16.2 17.7 18.3 17.9 19.9 19.6 19.4 13.0 13.1 13.0 2008 16.7 17.7 16.9 16.6 16.9 16.3 18.6 19.2 18.5 20.0 19.3 19.1 13.1 13.2 13.2 2009 16.6 16.9 16.7 16.2 15.8 15.8 18.7 18.6 18.7 18.6 18.4 18.1 13.5 13.7 14.4 2010 16.9 16.6 16.1 16.3 15.9 15.6 18.5 18.9 18.1 20.5 19.8 19.2 14.4 14.4 14.1 2011 15.0 16.0 14.8 15.8 15.9 16.2 17.0 17.6 17.4 18.2 17.8 17.6 13.6 13.9 13.7 2012 16.5 16.7 16.0 15.8 15.7 15.6 17.4 17.8 17.6 16.4 15.5 15.8 14.3 14.4 14.5 2013 14.9 14.4 14.2 14.6 14.4 14.7 17.2 17.2 16.7 17.2 16.7 16.1 13.7 13.7 13.5 2014 15.0 15.0 14.4 15.6 15.1 14.8 16.7 16.8 16.3 17.8 16.9 16.5 13.8 14.0 13.9 2015 17.2 16.1 16.0 14.9 14.9 14.8 16.4 16.6 16.7 18.1 17.1 16.5 14.2 14.0 14.4 2016 15.1 15.1 14.4 14.7 14.5 14.4 16.3 16.5 16.1 17.8 16.8 16.5 14.2 13.6 13.7 2017 15.2 15.1 14.5 14.8 15.1 15.2 17.2 17.8 17.2 19.2 18.1 17.7 14.6 14.2 14.1 2018 14.1 14.3 13.9 15.2 14.7 14.2 17.5 18.1 17.7 18.5 17.8 16.7 14.5 14.0 14.0 2019 13.7 13.8 14.0 15.1 14.6 14.6 16.9 17.3 16.8 17.0 16.8 16.4 14.8 14.5 14.3 2020 13.1 13.0 12.8 13.6 13.6 13.1 16.1 15.8 15.2 15.1 14.8 14.6 13.7 13.5 13.9 2021 12.9 12.9 13.1 14.1 13.7 14.1 16.2 15.7 15.9 16.7 16.1 15.4 14.2 14.4 14.1 2022 13.6 13.7 13.3 13 15.7 12.7 15.5 14.9 14.7 16.3 15.4 15.6 14.3 14.4 14.3 2023 15.2 15.2 14.6 13.6 13.3 13.3 14.0 13.4 14.3 16.3 15.6 15.4 14.0 14.3 14.3 and TLSTM models. Notably, including the TSL significantly improves the performance compared to FF in most cases, affirming the benefits of feature selection and model simplification for this problem. Overall, TLSTM achieves the best results in most datasets, followed by TFF. The L1 model obtains remarkable results, but the reliability is slightly worse than TFF and TLSTM excepting the Bermejales and Torneo datasets. It is crucial to emphasize that our primary interest lies in assessing the model’s efficacy for the last year, as it represents the most ‘‘realistic’’ scenario where only past data is used for training and future data are used for evaluation. In this context, the TLSTM model outperforms the L1 and TFF models in three datasets except for the Bermejales and Torneo datasets where TFF and L1 obtain the best results respectively. 4.5. Error map analysis In this section, we conduct a detailed analysis of the predictions obtained by our proposal for the year 2023, which is paramount in validating the model’s effectiveness. Fig. 4 illustrated the aggregated errors for each hourly prediction (yaxis) for all months (x-axis). In the Asomadilla, Bermejales, and Torneo datasets, the error is significantly lower in December compared to other months, with a similar trend observed for November and January. In Aljarafe and Ronda del Valle, January and December exhibit relatively low errors, but the distribution appears more scattered. Additionally, it is worth noting that errors at the beginning of each month are notably higher than on subsequent days. 8 M.J. Jiménez-Navarro et al. Fig. 4. RMSE aggregated for each month and hour of the day in Aljarafe, Asomadilla, Bermejales, Ronda del Valle, and Torneo. 4.6. Best/worst predictions analysis Fig. 5 displays each dataset’s best and worst 24-hour Ozone predictions. Generally, as discussed earlier, the best predictions are observed at the end and beginning of the year for most datasets, except for Asomadilla and Ronda del Valle, where the best predictions occur in June. Conversely, the worst predictions seems to occur around the middle of the year in Aljarafe and Torneo while the other worst predictions appear at the beginning. It is essential to highlight the potential dangers of make decisions considering the predictions when the predictions have a considerable error, leading to underestimating ozone values and poor recommendations. In conclusion, while the model produces good overall results, it is crucial to exercise caution when interpreting predictions during periods of increased error. In epochs and hours with relatively low error rates, the predictions appear to be accurate enough to be trusted. 4.7. Selected features analysis This section analyzes the features selected by our proposal to gain insight into the behavior learned by the model and the problem. Fig. 6 illustrates the selected features (y-axis) for each time step in the input window (x-axis) across all datasets. Two distinct behaviors are observed among Aljarafe, Asomadilla, Bermejales, and Ronda del Valle, in contrast to Torneo. The model opted for a larger window size in almost all datasets with a sparse feature selection. The selected features were close to 50% in all datasets, which suggests that most features in most time steps were not beneficial for forecasting. Specifically, Asomadilla was the dataset that required fewer features, with only 49.7% of the features selected. Surprisingly, all features showed relevance in at least one-time step, making the complete removal of any feature impractical. Notably, the most frequently selected feature across almost all datasets was Ozone as expected. On the other hand, average temperature (TMP), wind speed (WP), and NO2appeared to be less crucial than pollutants, which held greater importance. 4.8. Number of features distribution analysis Fig. 7 represents the number of features selected by the models along the different experiments performed during the hyperparameter optimization. The FF and LSTM models do not include a feature selection process, which always selects all available features. The differences between the experiments are caused by the window size optimization, which also influences the number of input features. The DT model has difficulty selecting a small number of features, probably when the optimization process is stuck in choosing a large depth for the tree. TLSTM, which performed the best in terms of RMSE, contains considerably more features than the L1 or TFF models. Generally, feature selection for the TSL in LSTM is more challenging than for the TFF model. Nevertheless, the TFF obtained remarkable results, which could also be a good option. Finally, the L1 and TFF models are the models that select fewer features in all datasets, with L1 being remarkably better. However, selecting too many features may cause the poor performance shown in Section 4.1, oversimplifying the model in some cases to near zero features. Even though TSL does not selects fewer features, there seems to be a more robust misconfiguration of the model and also obtained better efficacy than L1. 9 M.J. Jiménez-Navarro et al. Fig. 5. Best and worst predictions for all datasets. Prediction is represented in blue, while ground truth is represented in orange. 5. Conclusions and future works In this work, we present a practical application of TSL by enhancing the previous proposal in more than 1 μg∕m3in most datasets, creating a more reliable model with more accurate predictions of the behavior of ozone pollution. Therefore, a promising methodology has been provided for an efficient and environmentally aware traffic organization or users with respiratory problems. Our approach includes an embedded feature selection method that simplifies and reduces the problem dimensionality, thereby streamlining the model in cases where complexity is unnecessary to achieve good results. As a result, we reduce both computational costs and improve the interpretability of the model. Our approach will allow us to shed some light on deep learning models with a simple yet effective approach that helps us determine the relevant factors that influence the model prediction. TSL could be applied to most time series forecasting problems, not only limited to environmental data. Generally, in time series forecasting problems, less than 50% or 15% of inputs are relevant to obtaining good results. This sparse input information deviates from the typical approach but proves to be effective in our method. Notably, our approach not only yields remarkable results but also exhibits greater robustness to hyperparameters compared to other models in the comparison. Furthermore, TSL enables neural networks to handle a broader range of problems, regardless of their complexity. Notably, we obtained promising results in the Ozone forecasting problem, even though the solution might require a simpler model. Looking ahead, we have several ideas planned further to improve the results and interpretability of the model. Firstly, instead of sharing a single TSL for every output in a multi-step problem, we assign an independent TSL for each specific output. This tailored approach can lead to more accurate predictions, giving more precise information about the relevant features not just for the model but for every specific class that is not widely studied in the literature. The challenge here is to modify the architecture to include this mechanism efficiently without requiring extra layers/weights. Additionally, we consider the locality principle, where the moments near a specific moment should have a greater probability of inclusion in the model. This localized focus may further refine the accuracy of the predictions, reducing the error and producing a smoother instant selection. The main challenge is to determine the optimal number of near instances of every selected instant and its grade of influence over the near instances. Lastly, we plan to experiment with residual connections, as recent studies have shown that they can significantly enhance results. As observed in previous works, incorporating these connections into our model may yield even better efficacy. The challenge in this improvement is to modify the different regularization functions to consider the layers at every level.