Handling Missing Air Quality Data Using Bidirectional Recurrent Imputation for Time Series and Random Forest: A Case Study in Mexico City

1. Introduction

Air pollution refers to environmental contamination by any chemical, physical, or biological agent that modifies the natural characteristics of the atmosphere. It is primarily caused by motor vehicles, industrial activity, forest fires, and household combustion devices. Major air pollutants are particulate matter (PM₂.₅ and PM₁₀), carbon monoxide (CO), ozone (O₃), nitrogen dioxide (NO₂), and sulfur dioxide (SO₂). Air pollution causes respiratory and cardiovascular diseases, among others, and is an important source of morbidity and mortality, shortening the average person’s lifespan by 1 year and 8 months (World Health Organization (WHO), 2024 [1]). According to the State of Global Air Report 2024 [2], air pollution accounted for more than 1 in 8 deaths globally in 2021 (8.1 million deaths), becoming the second risk factor for early death, surpassed only by high blood pressure. Over 700,000 deaths in children under five were associated with air pollution that year, representing 15% of all global deaths in this age group. Furthermore, air pollution is a major risk factor for death from heart disease, diabetes, lung cancer, and other noncommunicable diseases (State of Global Air (SoGA) [2]). The World Health Organization estimated that 99% of the global population breathes air exceeding its guideline limits, with the highest exposures in low- and middle-income countries. Since many air pollutants are also major contributors to greenhouse gas emissions, reducing air pollution offers dual benefits, improving public health and supporting climate change mitigation efforts (WHO Air pollution [3]).

Mexico City has a metropolitan area that exceeds 22 million inhabitants and has complex air pollution problems, caused primarily by heavy vehicular traffic, industrial emissions, extensive urbanization, and the burning of fossil fuels (Kim et al., 2021 [4], He et al., 2023 [5]).

This situation is aggravated by topographical factors that restrict atmospheric dispersion, because Mexico City sits over 2000 meters above sea level and is surrounded by mountains on three sides. These conditions result in naturally lower oxygen levels that contribute to incomplete fuel combustion and the accumulation of both primary and secondary pollutants, including fine particulate matter (PM_2.5), ozone (O₃), and nitrogen dioxide (NO₂). Although government initiatives (ProAire [6]) have led to improvements over the past decades, key pollutants continue to exceed guidelines at many monitoring stations across the city (He et al., 2023) [5]).

Effective air quality management relies on accurate environmental monitoring. Mexico City Atmospheric Monitoring System [6,7] operates a network of monitoring stations that collect hourly data on pollutant (CO: Carbon monoxide), NO: Nitric oxide, NO₂: Nitrogen dioxide, NOₓ: Nitrogen oxides, NO + NO₂, SO₂: Sulfur dioxide, O₃: Ozone, PM₁₀, PM_2.5, PMCO: coarse fraction: PM₁₀–PM_2.5) concentrations and meteorological (RH: relative humidity, Temp: Temperature, WDR: Wind Direction, WSP: Wind Speed) parameters to monitor air quality in real time.

However, missing data is a recurring issue, caused by sensor malfunctions, maintenance, power outages, and network transmission failures (Zhang & Zhou, 2024 [8]; Wang et al., 2023 [9]). These issues compromise the completeness and reliability of the dataset, obstructing subsequent analyses such as trend evaluation, forecasting, and epidemiological modeling (Hua et al., 2024 [10]). Imputing missing values in air quality time series is therefore essential to ensure data completeness and analytical robustness (Alkabbani et al., 2022 [11]).

While traditional imputation methods are limited in capturing the complex spatio-temporal dependencies inherent in environmental data, recent advances in machine and deep learning have opened promising avenues (see Table 1). Multiple strategies have been developed to handle missing data in atmospheric monitoring systems. Traditional methods, such as mean/median substitution or linear interpolation, are simple but often fail in complex, multivariate contexts with non-linear dependencies or temporal structure (Wang et al., 2023 [9]). More advanced statistical approaches like Multivariate Imputation by Chained Equations (MICE) and regression-based methods offer improvements but are sensitive to assumptions about data distributions (Hua et al., 2024 [10]).

Machine learning techniques, including k-nearest neighbors (KNN), support vector regression (SVR), and decision trees, have shown enhanced accuracy, especially when capturing local or spatial patterns (Zhang & Zhou, 2024 [8], Wang et al., 2023 [9], Hua et al., 2024 [10]; Camastra et al. 2022 [12]). Hybrid and ensemble models have also emerged to better exploit spatial-temporal correlations (He et al. (2023) [5], Wang et al. (2023) [9]).

In recent years, deep learning has become central to missing data imputation in air quality research. RNN-based models such as Gated Recurrent Unit-Decay (GRU-D)( Che et al., 2018 [13]) and Bidirectional Recurrent Imputation for Time Series (BRITS) ( Cao et al., 2018 [14]) leverage temporal dependencies and trainable masking strategies, significantly improving performance on time series with irregular missingness. Similarly, generative models have been applied for their ability to more realistically reconstruct missing values. Among them, Generative Adversarial Imputation Nets (GAIN) ( Yoon et al., 2018 [15]) introduces adversarial training into the imputation task, improving robustness across diverse datasets. A recent review by Shahbazian and Greco (2023) [16] presents an extensive survey of several imputation methods based on Generative Adversarial Networks (GAN), covering applications in healthcare, finance, meteorology, and air quality. The review highlights the strengths of these models in preserving statistical distributions and generating more realistic imputations compared to traditional techniques. Other novel architectures have also been proposed, such as Graph Recurrent Imputation Network (GRIN) (Cini et al. 2021 [17]), a graph neural network imputation method, and Self-Attention-Based Imputation for Time Series (SAITS) (Du et al. 2023 [18]), a transformer-based model, which further expands the tools for handling missing data in structured environmental datasets.

Recent studies have explored the use of advanced machine and deep learning models to address missing data in air quality monitoring and to improve forecasting accuracy. A brief literature review of recent studies on air quality data imputation and prediction is presented below. Table 1 summarizes the key contributions in the field, with a focus on datasets, modeling techniques, study targets, and main findings.

Colorado Cifuentes & Flores Tlacuahua (2021) [19] proposed a deep learning approach for short-term air pollution forecasting in Monterrey, Mexico, using hourly data (2012–2017) of eight pollutants and seven meteorological variables from 13 monitoring stations. Seasonal decomposition was applied to impute missing values, and a feedforward Deep Neural Network (DNN) was trained to predict O₃, PM₂.₅, and PM₁₀ concentrations up to 24 hours in advance. A key limitation identified was the reduced accuracy of the model in predicting extreme pollution events.

Kim et al. (2021) [3] developed an interpretable deep learning model based on the Neural Basis Expansion Analysis for interpretable Time Series (N-BEATS) architecture to impute minute-level PM₂.₅ and PM₁₀ concentrations using data from 24 stations in Guro-gu and 42 stations in Dangjin-si, South Korea. Real missingness rates of 7.91% and 16.1% were reported, with an additional 20% of data artificially removed for testing. N-BEATS outperformed baseline techniques (mean, spatial average, MICE) and enabled interpretability through trend, seasonality, and residual components. However, it showed limitations in capturing long-term seasonal trends due to the use of fixed-period Fourier components.

Alahamade & Lake (2021) [20] evaluated six imputation strategies using hourly pollutant data (PM₁₀, PM₂.₅, O₃, NO₂) from the UK’s Automatic Urban and Rural Network (AURN) for the period 2015–2017. Their approach involved three groups of models: (1) clustering-based (CA: cluster average, CA+ENV: average within same station type, CA+REG: average within same region) using multivariate time series (MVTS) to group stations by temporal similarity; (2) spatial proximity-based (1NN: nearest station; 2NN: two nearest stations); and (3) ensemble: median of outputs of all five methods. The ensemble model performed best for O₃, PM₁₀, and PM₂.₅, while CA+ENV outperformed others for NO₂ due to its high spatial variability. Performance depended on pollutant and station type, and extreme values were slightly under- or overestimated.

He et al. (2023) [5] developed a predictive modeling framework for daily NO₂ concentrations in Mexico City by integrating ground-based (RAMA: Air quality monitoring network in Mexico City), satellite observations (OMI: Ozone Monitoring Instrument, and TROPOMI: TROPOspheric Monitoring Instrument), and meteorological data (CAMS: Copernicus Atmosphere Monitoring Service). Missing values were imputed using the MissForest algorithm. Random Forest and XGBoost outperformed the Generalized Additive Model (GAM), demonstrating the value of hybrid data sources for spatial-temporal NO₂ estimation.

Wang et al. (2023) [9] developed BRITS-ALSTM, a hybrid deep learning model that integrates BRITS with an attention-based Long Short-Term Memory (LSTM) decoder. The model was applied to hourly air quality data (2019–2022) from 16 monitoring stations in the alpine regions of China, covering six pollutants (PM₂.₅, PM₁₀, O₃, NO₂, SO₂, and CO). Both real missingness (5%–22%) and artificial gaps (30%) were considered. BRITS-ALSTM consistently outperformed classical and deep learning baselines, showing strong performance across pollutants and various missing data patterns.

Zhang and Zhou (2024) [8] proposed TMLSTM-AE, a hybrid deep learning model that combines LSTM, autoencoder architecture, and transfer learning to impute missing values in the PM_2.5 time series from Xi’an, China, characterized by single, block, and long-interval missing patterns. Artificial missingness ranging from 10% to 50% was introduced to evaluate performance. TMLSTM-AE was compared with both classical and deep learning approaches. TMLSTM-AE outperformed all baselines, especially in scenarios with long and consecutive missing gaps, demonstrating its capacity to capture spatio-temporal dependencies.

Hua et al. (2024) [10] evaluated several traditional imputation methods (mean, median, k-Nearest Neighbor Imputation (KNNI), and MICE) and deep learning models (SAITS, BRITS, MRNN, and Transformer) using six real-world air quality datasets from Germany, China, Taiwan, and Vietnam. The study introduced artificial missingness ranging from 10% to 80%. Results showed that SAITS and BRITS achieved the lowest imputation errors, while KNN performed well on large datasets with high missing rates (up to 80%). In contrast, MICE was more effective on smaller datasets but incurred higher computational costs. The Transformer-based model underperformed across most scenarios. The authors concluded that no single method is universally optimal; performance depends on dataset size, missingness rate, and application context.

As summarized in this analysis (Table 1), recent studies have increasingly adopted deep learning techniques, such as RNNs, attention-based models, and hybrid architectures, for imputing air quality data, particularly under complex and irregular missingness patterns.

BRITS has demonstrated high performance in reconstructing missing values in multivariate and temporally structured data, using bidirectional recurrent neural networks and a learnable masking mechanism. While some studies in Mexico have explored statistical and machine learning approaches for air quality prediction, the application of deep learning-based imputation remains limited. Conversely, tree-based ensemble methods such as Random Forest have been widely used for imputation in environmental datasets due to their flexibility and ability to handle nonlinear relationships, although they do not explicitly model temporal structure.

The objective of this study is to compare the performance of BRITS, a state-of-the-art deep learning model, with Random Forest, a robust machine learning baseline, for imputing missing hourly observations in air quality time series from Mexico City, considering multiple pollutants (CO, NO, NO₂, NOₓ, SO₂, O₃, PM₁₀, PM₂.₅, PMCO) and meteorological variables (RH, TMP, WDR, WSP). This study aims to identify the strengths and limitations of each approach and to contribute to the development of more accurate and reliable strategies for environmental data imputation.

2. Materials and Methods

2.1. Database

2.1.1. Study Area

Mexico City, the capital and largest metropolitan area in Mexico, is located at over 2,000 meters above sea level and is surrounded by mountains on three sides. In 2020, this city had a population of approximately 9.2 million, while the Valley of Mexico Metropolitan Area had 21.8 million inhabitants, with projections estimating 22.75 million by 2025 [21]. The Valley of Mexico Metropolitan Area includes Mexico City and its adjacent suburban areas. The city covers an area of 1,494.3 km² and faces persistent air quality challenges due to its topography, urban density, and emission sources.

2.1.2. Data Collection and Integration

This study uses hourly data from the Mexico City Atmospheric Monitoring Network (Red Automática de Monitoreo Atmosférico, RAMA [7]), comprising air pollutants (PM₂.₅, PM₁₀, NO, NO₂, NOₓ, SO₂, CO, and O₃) and meteorological variables (temperature, relative humidity, wind speed and wind direction). Data were collected from January 1, 2014, to December 31, 2023, including both dry and rainy seasons.

Annual files for each variable were loaded, concatenated, and transformed into a long format with standardized columns. Date and hour were combined into a single datetime column. All variables were then merged into a unified multivariate time series, aligned by monitoring station and timestamp. Missing values originally coded as -99 were replaced with NaN.

2.1.3. Missing Data Patterns Analysis and Station Filtering

A comprehensive analysis of the percentage of missing data per variable and monitoring station revealed several significant patterns. While some stations exhibit relatively low missingness across all variables, many show critical data quality issues. Several stations, including LAA, COY, and TPN, have complete (100%) missing data for most variables, making them unsuitable for time series modeling.

Additionally, PMCO, PM₂.₅, and PM₁₀ consistently have 100% missing values across many stations. In contrast, pollutants such as CO, NO₂, and O₃ generally exhibit lower percentages of missing data across a wider range of stations, indicating relatively higher data availability. Figure 1 shows a heatmap of the missing data percentages for all monitoring stations included in the dataset. Variables with more than 30% missing data are highlighted in red.

Due to the high proportion of missing values across many stations, a missingness threshold was used to select the stations included in the imputation procedure. Only stations with at least 70% of observed data for each variable were selected. Figure 2 shows the RAMA monitoring network, with the selected stations highlighted in blue. Some variables were excluded from the selected stations when they exceeded the missing data threshold, for example PM₁₀ at NEZ station and PMCO, PM₂.₅, and PM₁₀ at UAX, FAC, and NEX stations, which had complete 100% missing values.

2.1.4. Missing Data Patterns in Selected Monitoring Station

Figure 3 shows the total percentage of missing values per variable for each of the 10 selected monitoring stations (CUA, FAC, MER, MGH, NEZ, PED, SAG, TLA, UAX, VIF). The missingness rates vary across stations and variables ( i.e., some stations have more missing data than others, and within each station, certain variables (e.g., PM₁₀, PM₂.₅, and PMCO) may exhibit particularly high or complete missingness). Although most of these stations have less than 20% missing data overall, certain variables such as PM₁₀, PM₂.₅, and PMCO, show very high or complete absence of data in some stations. For this first stage of the imputation project, only the best-informed stations were selected to ensure the quality of imputed values.

To better understand the temporal structure of missingness, nullity matrices were generated for each station. As illustrated in Figure 4, these visualizations reveal prolonged gaps in data across multiple variables, suggesting that the missing values are not randomly distributed but may be associated with systematic sensor failures.

2.2. Model Training and Evaluation Pipeline

• Missing data identification: Binary masks were created to identify real missing values in the dataset.
• Data splitting: Each station’s dataset was divided into 80% for training and 20% for testing.
• Hyperparameter optimization: Hyperparameter tuning was performed in sequential steps due to computational constraints.
• Artificial missingness for evaluation: To enable controlled performance assessment, 20% of the observed values in the test set were randomly removed under a Missing Completely At Random (MCAR) assumption. This masking allowed direct comparison between imputed and ground-truth values.
• Normalization: For the BRITS method, all variables were standardized using MinMaxScaler normalization, applied separately to the training and test datasets.
• Model Training and Evaluation: After hyperparameters optimization, models were retrained using the full training set and evaluated on the masked test set. Imputation performance was assessed using several evaluation metrics including MAE, RMSE, Wasserstein Distance, and TOST equivalence tests.

2.3. Model Training

To assess imputation quality across different methodological paradigms, two models were implemented: a classical machine learning method (Random Forest) and a deep learning model (BRITS). These models provide contrasting approaches to reconstructing missing environmental data, ranging from ensemble-based predictions to recurrent architectures that can model temporal dependencies.

2.3.1. Random Forest (RF) Hyperparameter Optimization

The Random Forest Regressor offers a robust and interpretable solution for multivariate imputation. We used IterativeImputer from scikit-learn with a Random Forest estimator to iteratively predict missing values using the available features at each timestamp. This approach can model complex nonlinear dependencies and often achieves high accuracy.

Hyperparameter tuning was performed manually using a staged grid search, due to computational constraints. This process used five steps to optimize one hyperparameter at a time:

• n_estimators: Number of trees in the ensemble.
• max_depth: Maximum depth of individual trees.
• max_iter: Number of iterations in the iterative imputation process (specific to IterativeImputer).
• min_samples_split: Minimum number of samples required to split an internal node.
• min_samples_leaf: Minimum number of samples required to be at a leaf node.

The default values and evaluated ranges for these hyperparameters are presented in Table 2, section (a).

2.3.2. BRITS Hyperparameter Optimization

BRITS is a deep learning model designed specifically for imputation data in multivariate time series, without relying on linear dynamics or strong prior assumptions. It is based on a bidirectional recurrent architecture that processes time series data in both forward and backward directions using Recurrent Neural Networks (RNNs), allowing the model to learn temporal dependencies from both past and future contexts (Cao et al. [14]).

Unlike traditional methods such as MICE, where missing values are pre-filled and treated as known inputs, BRITS treats missing values as latent variables embedded within the model’s computation graph. These are dynamically updated through backpropagation during training, enabling more accurate and context-aware imputations. As reported by Hua et al. [10], BRITS has shown strong performance on air quality datasets by effectively capturing nonlinear and correlated missing patterns.

Hyperparameter tuning was carried out manually in three sequential stages using a grid search strategy. In each stage, a subset of parameters was optimized while others were fixed:

• Stage 1: RNN units and subsequence length were tuned, with learning rate = 0.005, batch size = 64, use_regularization = False, and dropout_rate = 0.2 held constant.
• Stage 2: learning rate and batch size were tuned, keeping use_regularization = False and dropout_rate = 0.2 fixed.
• Stage 3: use_regularization and dropout_rate were jointly optimized.

Brief descriptions:

• RNN units: Number of hidden units in the RNN layers. Higher values improve model capacity to learn temporal patterns.
• subsequence length: Input window size used during training.
• learning rate: Controls the optimizer’s step size.
• batch size: Number of samples per training batch.
• use_regularization: Whether to apply dropout/L2 to prevent overfitting.
• dropout_rate: Proportion of neurons randomly deactivated during training.

The default values and evaluated ranges for these hyperparameters are presented in Table 2, section (b).

After hyperparameter tuning, the final model for each station was trained using 80% of the data and the corresponding optimal hyperparameter values. Model performance was then evaluated on the remaining 20% test set. To assess imputation accuracy, a mask was applied to randomly hide 20% of the observed values within the test set, simulating missing data. The model’s imputations were compared against these masked ground-truth values to compute the evaluation metrics.

2.4. Model Evaluation

To evaluate the imputation performance, four metrics were used: Mean Absolute Error (MAE), Root Mean Squared Error (RMSE), Two One-Sided Tests (TOST) AND Wasserstein Distance.

• Mean Absolute Error (MAE) quantifies the average absolute difference between imputed values ŷₜ and observed values yₜ:

$$MAE={\displaystyle \frac{1}{n}}\sum _{t=1}^n |yₜ-ŷₜ|$$

(1)

Where $n$ is the number of data points.

2.

Root Mean Squared Error (RMSE) measures the square root of the average squared differences between imputed and observed values:

$$RMSE=\sqrt{{\displaystyle \frac{1}{n}}\sum _{t=1}^n {(yₜ-ŷₜ)}^2}$$

(2)

3.

Wasserstein Distance, also known as Earth Mover’s Distance (EMD), measures the dissimilarity between two probability distributions by calculating the minimum effort required to transform one distribution into another [23]. In the context of air quality time series, it is useful for comparing distributions of imputed and observed values. Given two cumulative distribution functions, $P\left(x\right)$and $Q\left(x\right)$, the first-order Wasserstein distance is defined as:

$$W(P, Q) = \int \left|P\right(x) -Q(x\left)\right| dx$$

(3)

where the variable $x$ represents all possible values that a specific air quality variable (e.g., PM₂.₅ concentration) can take over time. The integral computes the area between the two cumulative curves, summarizing the total difference in distribution between the imputed and observed values. A smaller Wasserstein Distance means that the imputed values follow the same pattern as the real ones, showing that the imputation method maintains the original data characteristics well.

4.

The Two One-Sided Test (TOST [22]) is a statistical procedure used to assess the equivalence between the mean of imputed and observed values. Unlike conventional statistical tests that seek to detect significant differences, TOST is specifically designed to demonstrate equivalence, that is, to confirm that the difference is small enough to be considered negligible within a predefined tolerance margin. Let μₒ and μᵢ be the means of observed and imputed values, respectively. According to Santamaría-Bonfil et al. [2], TOST is implemented as two one-sided t-tests:

$$H₀: (\mu ₒ-\mu ᵢ) \le -\Delta and H₁: (\mu ₒ-\mu ᵢ) \ge \Delta$$

(4)

where Δ is a defined equivalence margin. Equivalence is declared only if both null hypotheses are rejected, confirming that the mean difference lies entirely within the interval ($-\Delta ,+\Delta$). This criterion ensures that the imputed values do not statistically differ from the observed values beyond an acceptable tolerance. In this study, we used a significance level α = 0.05, and equivalence margins based on ±10% of the observed mean.

3. Results

3.1. Results of Hyperparameter Optimization

Hyperparameter tuning was performed in sequential steps as described in section 2.3 “Model Training”. For Random Forest, hyperparameter tuning revealed notable variability across stations (Table 3). The number of estimators (n_estimators) ranged from 10 to 50, while the maximum tree depth (max_depth) was commonly set to 10 or 20, though some cases used no depth limit. The number of imputation iterations (max_iter) varied from 5 to 20, and the minimum samples required to split a node (min_samples_split) and to remain at a leaf node (min_samples_leaf) showed diverse configurations, indicating sensitivity to local data structure. Overall, smaller leaf sizes (1 to 6) and moderate (2 to 15) split thresholds were preferred.

The optimized hyperparameters suggest that stations with more complex patterns (e.g., PED, SAG) required deeper trees (max_depth = 20) and more estimators (50), while others (e.g., CUA, VIF, MGH) performed better with shallower models (max_depth = 10) and fewer estimators (10–20). Variations in min_samples_split and min_samples_leaf also suggest differences in noise levels and local dynamics across stations.

For BRITS, optimal hyperparameter settings were more consistent (Table 4). A large hidden dimension was generally preferred, with 512 RNN units for most stations. Subsequence length values commonly ranged from 32 to 64, and the learning rate was predominantly fixed at 0.005. The batch size was mostly set to 32, indicating a trade-off between stability and speed. Regularization was selectively applied, and dropout was commonly set to 0.2, except for a few cases where 0.1 or 0.3 yielded better results.

BRITS showed consistent optimal settings across stations, favoring large hidden layers and moderate subsequence lengths, whereas Random Forest required highly station-specific configurations, reflecting its adaptability to local variability but also its sensitivity to irregular missingness patterns.

3.2. Results of Imputation Models

After hyperparameter optimization, the final models were retrained on the full 80% training set using the selected optimal hyperparameter configurations. Evaluation was performed on the 20% test set. Additionally, to simulate realistic missingness and allow controlled performance evaluation, 20% of the observed values in the test set were randomly masked. The imputed values were compared to these hidden ground-truth values using four metrics: Mean Absolute Error (MAE), Root Mean Squared Error (RMSE) on masked subsets, and Wasserstein Distance and TOST for equivalence testing on complete distributions.

3.2.1. Performance Evaluation on Masked Subset

MAE and RMSE metrics were computed using 20% of observed values randomly masked from the validation set.

MAE Results

Figure 5, Figure 6 and Figure 7 show the distribution of imputation errors, expressed as MAE and RMSE, for all variables across the ten stations, comparing RF and BRITS models. Each subplot corresponds to a specific variable with its original scale, because these metrics are reported in the original units of each variable and must be compared only within the same context.

Figure 5 presents MAE results for RF and BRITS across all stations within each variable. RF consistently achieved lower MAE values than BRITS for every variable. The largest differences occurred in reactive gases such as NO, NO₂, and NOx, where RF errors were less than half those of BRITS. For example, in station MER, the MAE for NOx was 8.66 with RF versus 26.62 with BRITS, while in TLA, PM₁₀ showed 4.85 for RF compared to 18.14 for BRITS (Table A1).

RMSE Results

Figure 6 shows a similar trend for RMSE. RF had smaller imputation errors and lower dispersion across stations compared to BRITS. RMSE differences were particularly notable for pollutants with high variability, such as NO, NOx, and PM10. For example, in station MER, RF achieved an RMSE of 20.36 for NOx, while BRITS reached 38.50; in TLA, PM₁₀ showed 10.77 versus 24.25, respectively. Because RMSE penalizes large deviations quadratically, these higher values for BRITS suggest greater sensitivity to extreme imputation errors compared to RF.

3.2.2. Distributional Similarity Assessment

Wasserstein Distance and TOST were applied to the complete distribution of each variable, comparing the original dataset with the imputed dataset. These complete distributions showed a different trend than the masked subset.

Results of Wasserstein Distance

BRITS achieved lower Wasserstein Distance values than RF for most variables, indicating superior preservation of the original data statistical distribution (Figure 7). For example, in station MER, NOx had 0.80 for RF compared to 0.40 for BRITS, and in station TLA, PM_2.5 showed 0.95 for RF versus 0.45 for BRITS. Similarly, in station SAG, NO₂ presented 1.35 with RF versus 0.64 with BRITS, and in CUA, NO₂ was 1.22 for RF against 0.59 for BRITS. BRITS also exhibited lower variability across stations, reflecting more consistent alignment with the true distributions. RF, despite achieving the best local accuracy (MAE and RMSE), displayed higher Wasserstein Distance values in many cases, underscoring a trade-off between decreasing individual errors and preserving the complete statistical distribution.

Results of TOST

Figure 8 presents a comparison between original means, imputed means, and TOST equivalence limits for each variable and station. In all cases the TOST equivalence criterion was satisfied, indicating that the imputed means were statistically equivalent to the original means within the predefined tolerance. Each subplot corresponds to a specific variable, while each row within a subplot represents one station. Black open circles indicate the observed mean for the complete dataset, while blue and red open circles represent the means obtained from the RF and BRITS imputation models, respectively. Horizontal dashed lines indicate the TOST equivalence bounds, defined as ±5% of the original mean for each variable. These bounds represent the range within which the imputed mean would be considered statistically equivalent to the original mean. This visualization facilitates a clear assessment of the proximity of the imputed means to the observed mean across variables and stations.

The comparison of differences between observed and imputed means across stations reveals several patterns in the performance of RF and BRITS. Overall, most differences remained small (< ±1 unit), indicating that both models achieved high accuracy for most variables. However, BRITS showed a slight overestimation tendency for several pollutants, particularly for NO, NOx, PM₁₀, and PM_2.5. While RF generally was stable, it exhibited higher variability in some cases, with both negative and positive deviations, especially in variables with greater heterogeneity, such as wind direction (WDR).

Visualization of Kernel Density Estimation (KDE)

KDE is a technique used to estimate the probability distribution of a dataset in a smooth and continuous manner. It is widely applied in exploratory analysis to examine the shape of the distribution (e.g., symmetric, skewed, multimodal). In visualization, a KDE plot is a curve that represents the estimated density derived from the data.

Figure 9 shows KDE-based estimated distributions for the MER station, comparing original data (shaded area) with imputations from RF (black dashed line) and BRITS (red dashed line). Overall, both techniques adequately captured the global shape of the original distributions. However, notable differences are observed for variables with high skewness or multimodality. For pollutants such as CO, NO, NO₂, and NOx (Figure 9 a–d), the imputed distributions closely followed the right-skewed pattern of these variables, although RF tended to generate sharper peaks around central tendency, whereas BRITS provided smoother densities, more closely resembling the original distribution. For SO₂ and O₃ (Figure 9 e, f), both methods exhibited acceptable alignment. For particulate matter (PM₁₀, PM₂.₅, and PMCO; subplots g–i), notable differences emerged: RF concentrates density around central values, reducing distribution spread, while BRITS better reproduced distribution. For meteorological variables such as RH and TMP (Figure 9 j, k), BRITS more accurately captured RH bimodality and TMP variability, whereas RF tended to overfit central regions. Finally, for WDR and WSP (Figure 9 l, m), both methods followed the general pattern, though marked deviations occurred in WDR due to its multimodal distribution, where RF exhibited oscillations absent in the original data. Overall, Figure 9 highlights that BRITS preserved distributional characteristics more effectively in variables with high variability, while RF, despite achieving reasonable central estimates, tended to distort tails and reduced spread in highly skewed variables.

Figure A1, Figure A2, Figure A3, Figure A4, Figure A5, Figure A6, Figure A7, Figure A8 and Figure A9 present KDE-based estimated distributions for the CUA, FAC, MER, MGH, NEZ, PED, SAG, TLA, UAX, and VIF stations, comparing the original dataset with imputations from RF and BRITS models. KDE analysis showed that BRITS generally preserved distributional variability more effectively than RF, which tended to concentrate density near central values; both models struggled with WDR due to its multimodal nature.

Figure 10 and Figure 11 show the complete time series from 2014 to 2023 for the MER station, which exhibited the low missing rates, ranging from 4% (WDR and WSP) to 16% (PM₁₀, PM₂.₅, and PMCO). Figure 10 shows imputations obtained for BRITS, where imputed values (yellow color) closely follow the pattern of observed values (blue color), even across large gaps, for example, TMP in 2021 and NO, NO₂, and NOx during in late 2014 display imputations that align well with the original trends. In contrast, RF imputations (Figure 11) were less consistent, particularly for particulate matter variables (PM₁₀, PM₂.₅, and PMCO), where the model often generated unrealistic extremes at peak points.

4. Discussion

Our findings demonstrate that BRITS generally outperforms RF in maintaining distributional similarity and minimizing Wasserstein distance, particularly for variables with high variability, while both methods achieved TOST equivalence. These results are consistent with recent studies by Wang et al. (2023) and Hua et al. (2024), which highlighted the robustness of BRITS and other deep learning models under irregular missingness in multivariate time series. RF performed adequately for low-variability variables (e.g., CO and SO₂), confirming its usefulness in certain contexts as noted by He et al. (2023). However, both methods showed limitations when imputing meteorological variables with higher variability, such as wind direction (WDR).

Neither RF nor BRITS successfully reconstructed long periods of missing data. For example, there was a extreme gap in PMCO, PM₁₀, and PM₂.₅ data at the PED station from 2016 to 2019. During this period, BRITS consistently underestimated the actual values, while RF tended to overestimate and produced outliers.

An important aspect to consider is the contrasting hyperparameter behavior of both models, which reflects their underlying adaptability. BRITS exhibited relatively consistent optimal settings across stations, favoring large hidden layers (512 units), moderate subsequence lengths (32–64), and a fixed learning rate (0.005). This indicates its reliance on generalized temporal patterns. In contrast, RF required highly station-specific configurations, such as deeper trees (max_depth = 20) and more estimators (50) for stations with high variability (e.g., PED and SAG), while simpler structures sufficed for others. This heterogeneity suggests that RF adapts its complexity to local dynamics, whereas BRITS applies a more uniform strategy, which may explain its superior stability for short to medium gaps but inferior performance under extreme missingness.

Additional preliminary tests were conducted using two alternative deep learning models, specifically SAITS and GAIN. However, these models produced unsatisfactory results, because they exhibited very limited variability and a strong tendency to reproduce mean values, rather than capturing the original distribution or local dynamics. This smoothing effect makes such models unsuitable for time series in which maintaining natural variability is critical, such as environmental monitoring and subsequent forecasting.

Future research on this topic should focus on developing more robust imputation frameworks capable of handling stations with extreme missingness (greater than 30% and even 100%), which were excluded from this study. This could involve hybrid deep learning and machine learning methods that integrate spatial-temporal dependencies to preserve variability and dynamic patterns in long missing sequences.

5. Conclusions

This study evaluated the performance of RF and BRITS for imputing missing values in multivariate air quality time series from Mexico City monitoring stations. The results evaluation considered two perspectives: local performance on masked subsets using MAE and RMSE metrics, and distributional similarity on complete datasets through TOST and Wasserstein Distance. Results indicate that RF consistently achieved lower MAE and RMSE values than BRITS on masked subsets. In contrast, BRITS generally preserved the statistical distribution of the original data better than RF, for most variables. However, both models showed limitations for variables with high variability, such as wind direction (WDR) and particulate matter (PM₁₀ and PM₂.₅), where errors were larger and equivalence intervals were wider.