A Hybrid LSTM-Attention Model for Multivariate Time Series Imputation: Evaluation on Environmental Datasets

1. Introduction

The rapid expansion of machine learning (ML) applications has dramatically increased the demand for large-scale, high-quality datasets that enable reliable and generalizable model training. The growing accessibility of cost-effective and scalable sensing technologies, such as Internet of Things (IoT) devices, remote sensing nodes, and wireless sensor networks has transformed how data are collected across domains. These systems continuously generate time-stamped data in areas ranging from healthcare monitoring and environmental systems to industrial automation and urban infrastructure, thereby providing a rich source of information for predictive modeling and decision support. However, this exponential increase in data acquisition also amplifies challenges related to data completeness, reliability, and consistency, which directly affect the robustness of ML-driven insights. In particular, data missingness has emerged as one of the most prevalent and detrimental issues, often arising from sensor drift, calibration errors, power failures, communication disruptions, or environmental interferences.

Missing or incomplete data is not confined to any single domain but appears across virtually all data-driven fields. It is particularly common in healthcare systems (Li et al. 2021; Wells et al. 2013; Zuo et al. 2021), geographic information systems (GIS) (Griffith and Liau 2021), transportation and traffic management (Kaur et al. 2022), environmental surveillance networks (Decorte et al. 2024; Okafor and Delaney 2021), and industrial process monitoring (Carbery et al. 2022). In these contexts, data gaps can lead to unreliable trend detection, erroneous forecasts, and reduced interpretability of analytical models. Furthermore, incomplete datasets can introduce statistical bias, distort model parameters, and reduce predictive accuracy when used in training or validation stages (Escobar et al. 2020; Hasan et al. 2021). These challenges are magnified in time-dependent domains, where data continuity and temporal correlations are essential for learning meaningful representations.

To mitigate the effects of missingness, two general strategies have been traditionally adopted, that are downsampling and imputation (Emran 2015). Downsampling eliminates incomplete observations entirely, simplifying analysis but at the cost of discarding potentially valuable information. This approach can be especially harmful in time series analysis, where removing data points disrupts temporal dependencies that are crucial for forecasting and anomaly detection (Okafor and Delaney 2021). In contrast, imputation methods attempt to estimate or reconstruct missing values using statistical inference, probabilistic reasoning, or machine learning-based modeling (Das et al. 2019). The effectiveness of an imputation approach depends not only on the chosen algorithm but also on the mechanism of missingness, which can be MCAR (Missing Completely at Random), MAR (Missing at Random), or MNAR (Missing Not at Random) (Brown and Kros 2003; Jäger et al. 2021; Kang 2013). Understanding this mechanism is vital for selecting appropriate techniques that minimize bias and preserve underlying data distributions.

Traditional imputation techniques are often categorized into Single Imputation (SI) and Multiple Imputation (MI) frameworks. SI approaches, including mean substitution, median filling, and zero imputation, offer simplicity and computational efficiency but tend to underestimate data variability and uncertainty (Kang 2013). MI, introduced through iterative modeling, generates multiple plausible versions of the complete dataset, capturing the inherent uncertainty in the imputation process (Jamshidian and Mata 2007). Classical statistical models such as regression, Expectation-Maximization (EM), and stochastic simulation have been widely used; however, their performance is typically limited when facing nonlinear or high-dimensional data structures (Brown and Kros 2003; Emmanuel et al. 2021). Consequently, recent years have witnessed a paradigm shift toward machine learning and deep learning-based solutions that can capture complex interdependencies among variables and temporal patterns.

A substantial body of research has benchmarked both conventional and learning-based imputation strategies across diverse applications. Kang (Kang 2013) surveyed traditional clinical data imputation methods such as Last Observation Carried Forward (LOCF), regression, and EM, while subsequent studies explored data-driven techniques including K-Nearest Neighbour (KNN), Random Forests (RF), Decision Trees (DT), Support Vector Machines (SVM) and ensemble learners (Emmanuel et al. 2021). Comparative analyses have demonstrated that model performance varies significantly with dataset characteristics and the type of missingness, suggesting there is no universal best method. For instance, Hegde et al. (Hegde et al. 2019) observed that Probabilistic Principal Component Analysis (PPCA) outperformed Multiple Imputation by Chained Equations (MICE) in healthcare datasets with MCAR missingness, while Stekhoven et al. (Stekhoven and Bühlmann 2012) showed that the non-parametric missForest algorithm achieved superior results on mixed-type biological datasets. Early neural network-based methods also demonstrated potential advantages: Gupta et al. (Gupta and Lam 1996) reported better results using backpropagation networks than with mean or regression imputation, and Silva-Ramírez et al. (Silva-Ramírez et al. 2011) validated multilayer perceptrons (MLPs) for robust imputation across both simulated and real-world datasets. More advanced models, such as deep autoencoders and variational autoencoders (VAEs), have since emerged as powerful tools capable of reconstructing complex patterns of missingness in nonlinear, high-dimensional spaces (Beaulieu-Jones and Moore 2017; McCoy et al. 2018).

With the rapid evolution of IoT infrastructures and real-time sensing systems, time series imputation has become a central research theme. Unlike cross-sectional data, time series are inherently sequential and temporally correlated, making traditional static imputation methods suboptimal. Decorte et al. (Decorte et al. 2024) investigated the performance of recurrent models such as Multi-Directional RNNs (M-RNN) and Bidirectional Recurrent Imputation for Time Series (BRITS) in environmental monitoring data with simulated MCAR conditions. Maksims et al. (Kazijevs and Samad 2023) benchmarked multiple deep learning architectures, including attention-based autoencoders, LSTM variants, and adversarial RNNs, across five healthcare datasets. Okafor and Delaney (Okafor and Delaney 2021) further demonstrated that variational autoencoder frameworks achieved superior reconstruction accuracy and improved calibration in greenhouse gas sensor networks when compared with traditional algorithms such as MICE, missForest, and KNN.

Recently, researchers have explored hybrid frameworks that integrate the strengths of deep sequential models with complementary algorithms. Examples include coupling Long Short-Term Memory (LSTM) networks with Decision Trees (Haibin and Yongliang 2023), optimizing learning with Genetic Algorithms (Tsoulos et al. 2024), or employing Transfer Learning strategies, for example, TrAdaBoost to handle domain adaptation challenges (Chen et al. 2021). Moreover, advances in attention mechanisms pioneered by transformer architectures in natural language processing (Domhan 2018; Strubell et al. 2018; Tang et al. 2018; Tao et al. 2018; Vaswani et al. 2017; Zhang and Thorburn 2021), have led to a surge of interest in attention-augmented imputation models. These models effectively learn context-aware temporal representations, capturing both short- and long-range dependencies within multivariate sequences.

Several studies and research experiments that have used hybrid methodologies to solve a variety of complex problems, particularly focusing on the innovative use of transformer mechanisms in conjunction with other neural network architectures, have greatly inspired us to explore further for imputing missing data in environmental datasets. A closely associated study conducted by Dagtekin et al. provides evidence that supports the potential efficacy of the combination of these techniques, as they present their findings of better imputation performance and prediction accuracy for water quality data through the integration of LSTM networks alongside self-attention mechanisms (Dagtekin and Dethlefs 2022). Similarly, Shang et al. have investigated the combination of LSTM, attention and AdaBoost algorithms to impute missing data in traffic flow dataset and have reported promising results (Shang et al. 2024).

LSTM and attention hybrids are used extensively in NLP or healthcare but they are underexplored in environmental IoT monitroing. Secondly, many studies use attention mechanism either before or after the LSTM layer(s). In this paper, we present a hybrid LSTM-Attention model for missing data imputation in environmental monitoring datsets. We present a novel "sandwiched" LSTM - Attention - LSTM design. This interleaving enables the first LSTM layer to capture local temporal dependencies, the attention mechanism to highlight long-range relationships, and the second LSTM to re-integrate these enriched features into a coherent sequential representation. This ordering was deliberately chosen to balance sequential memory retention with global context modeling. Many prior works train models on complete data and then test performance on missing data. In our study, we trained the proposed model and other baseline models on 30% missing data to align it with real-world scenarios. We test imputation performance through MAE and RMSE across 10%-90% missingness, which provides a broad robustness analysis on low to moderate missingness. The proposed model demonstrates that it preserves accuracy under extreme missingness scenarios (70-80%), which is beneficial in real-world events of IoT sensor failures.

In the subsequent section of this paper, we briefly introduce the problem and highlight the purpose and significance of the work. The current state of the research in the problem domain is briefly reviewed and this section end with highlighting the contributions of this study. The next section described the proposed model in detail along with sufficient details of the data and mthodology adopted. Section 3 provides concise experimental results and their interpretations. The last section discusses the findings and their implications along with the limitations of the work. Future research directions are also mentioned in this section.

2. Materials and Methods

2.1. Dataset Description

To evaluate the performance and generalizability of the proposed hybrid model, multiple environmental monitoring datasets were used. These datasets represent varied sensing scenarios, including hydrological, meteorological, and air quality measurements collected from field-deployed sensor networks. Each dataset contains multivariate time series data with varying temporal resolutions, enabling a comprehensive assessment of the model’s imputation capabilities. For each dataset, the target variable, that is the feature to be predicted, is identified alongside a set of correlated auxilliary parameters that serve as input features. The relationships among these variables are visualized using correlation matrices, which help to illustrate the degree of association between the target and auxiliary features. Detailed descriptions of each dataset are provided in the following subsections.

2.1.1. Soil Temperature Dataset

The soil temperature dataset (SoilTemp) is a subset of the dataset originally published by (Decorte et al. 2024). This dataset comprises continuous time-series measurements of soil temperature recorded at three different heights using in-situ soil sensors proposed in (Wild et al. 2019). The sampling frequency for the sensor is 15 minutes in degrees Celsius, offering high temporal frequency and the time series spans 6 months of data observations from 19 April 2021 at 00:00:00 to 26 September 2021 at 23:45:00, providing a clearly defined study period.

The target feature in this dataset is the soil temperature measured at 12 cm above ground level named as $temp\_top$ in the dataset, as it exhibited the greatest variability among the three temperature sensors. The variability of the measurements makes it a suitable candidate for training and evaluating the proposed model, thereby enhancing the generalizability of the obtained results. Input features include soil temperature measurements at ground surface level ($temp\_mid$), subsurface soil temperature ($temp\_bot$) 10 cm deep in the ground and relative humidity. The heatmap in Figure 1 shows the correlations between the target and input features. The heatmap inidcates strong correlation between the three temperature measurements, suggesting they can be useful predictors for the target feature. In comparison to this, a weak correlation between the traget feature and the relative humidity is present in the dataset. However, the target feature shows the strongest correlation with the humidity among the three temperature sensors.

2.1.2. Meteorological Dataset

The meteorological dataset is collected from an automatic weather station operating in the monitoring region under the SuDS+ project in Stanley, Durham, UK. The dataset contains records of nine different weather parameters at an hourly temporal resolution. The time series spans 16 months of data recorded from 11 June 2024 to 14 october 2025 to capture temporal patterns.

The target feature predicted in this dataset is the average air temperature measured in degrees Celsius which is a key climatic parameter. Other predictive features include average humidity, wind direction, wind speed, dew point, wind gust, heat index, wind chill and average rainfall in milimeters. As expected, temperature exhibits negative correlation with humidity and a strong positive correlation with heat index, dew point and wind chill parameters. The correlation matrix in Figure 2 highlights these relationships, which the proposed hybrid model leverages to improve the temporal estimation and imputation of air temperature values.

2.1.3. Air Quality Dataset

The air quality dataset utilized in this study is derived from the publicly available dataset presented by (Feinberg et al. 2018), hosted on the U.S. Environmental Protection Agency (EPA) Environmental Dataset Gateway (EDG EPA 2025). The dataset comprises measurements collected from low-cost gas sensors co-located with Federal Equivalent Method (FEM) reference monitors at an urban regulatory site in Denver, Colorado, USA. The monitoring period spanned approximately six months, from 8 September 2015 to 22 February 2016, with data aggregated at 15-minutes intervals.

Two target features are considered for analysis: ozone (O₃) and combined nitrogen dioxide/ozone (NO₂/O₃) concentrations, both expressed in parts per billion (ppb). For both target features, measurements were obtained for corresponding temperature, relative humidity, and reference measurements from the FEM monitors. The gas sensors demonstrated moderate to strong correlations with other auxilliary features. Specifically, the correlation coefficients between O₃ sensors and their reference exceeded 0.9, while those for NO₂/O₃ sensors ranged between 0.32 and 0.57. As illustrated in the correlation heatmap (Figure 3), the O₃ sensors showed strong inter-sensor consistency and positive associations with temperature, alongside negative correlations with relative humidity. These relationships provide valuable contextual dependencies that the proposed hybrid model can exploit for accurate imputation and prediction of missing air quality measurements.

2.2. Proposed Hybrid Model

The hybrid model presented in this paper is a sequential neural network designed for time-series imputation tasks, using a combination of LSTM layers and a Multi-Head Attention mechanism to capture temporal dependencies and enhance feature representation. Figure 4 shows the overall architecture of the proposed model. The proposed architecture sandwiches the attention layer between the two LSTM layers instead of cascading the Attention and LSTM layers on top of each other. The model maps an input sequence ${\mathbf{X}}_{1:T}$ to a scalar prediction $\widehat{y}$ via a first LSTM encoder producing hidden states $\mathbf{H}$. A multi-head self-attention module is then applied to $\mathbf{H}$, returning attention outputs $\mathbf{A}$. The concatenation $\mathbf{C}$ of outputs from first LSTM and Attention layers is passed through a second LSTM, which re-encodes $\mathbf{C}$ to produce the final hidden state ${\mathbf{s}}_T$. The output is obtained by a dense linear layer, and the network is trained by minimizing the mean squared error (MSE) as loss measure.

Before feeding the data into the model, missing values in the time series are reconstructed using linear interpolation to ensure temporal continuity. Given two observed data points $(t_1,x_{t_1})$ and $(t_2,x_{t_2})$, the interpolated value at time t, where $t_1<t<t_2$, is computed as

$$x_t=x_{t_1}+(x_{t_2}-x_{t_1})·{\displaystyle \frac{t-t_1}{t_2-t_1}}.$$

(1)

This method provides a simple yet effective way to fill small gaps in the data while preserving the temporal trend.

The interpolated data are subsequently normalized using Min-Max scaling to restrict all feature values within the range $[0,1]$, which improves numerical stability and accelerates convergence during training. Each feature value $x_t$ is scaled as

$$x_t^{\prime }={\displaystyle \frac{x_t-x_{min}}{x_{max}-x_{min}}},$$

(2)

where $x_{min}$ and $x_{max}$ represent the minimum and maximum values of the feature, respectively.

To capture sequential dependencies, feature engineering is applied using a sliding window technique with a window size of 30, denoted as T and represented by the parameter $window\_size$ in this paper. Each input sample at time step t is expressed as

$$X=[x_1,x_2,\dots ,x_T],$$

(3)

where $x_i\in {\mathbb{R}}^d$, with $d=1$ for univariate data and $d>1$ for multivariate cases. This preprocessing pipeline ensures that the input sequences provided to the model are both complete and properly scaled, enabling effective learning of temporal relationships across multiple correlated features.

2.2.1. First LSTM Layer (LSTM₁)

The preprocessed input is then passed to the first LSTM layer to learn long-term dependencies in sequential data (Hochreiter and Schmidhuber 1997). It utilizes 64 hidden units to learn local temporal dependencies and transforms the input sequence into a sequence of hidden states:

$$h_t^{\left(1\right)}={\mathrm{LSTM}}_1(x_t,h_{t-1}^{\left(1\right)},C_{t-1}^{\left(1\right)}).$$

(4)

In our implementation, the standard hyperbolic tangent activation functions used in the candidate and output computations of the LSTM were replaced by ReLU to enhance gradient propagation and model stability. The overall LSTM formulation and gating mechanism remains consistent with the original definition by (Hochreiter and Schmidhuber 1997). Accordingly, the operations for the forget gate, input gate ,candidate state update and output gate are defined as the following equations.

$$\begin{array}{cc}\hfill f_t& =\sigma (W_f[h_{t-1},x_t]+b_f)\hfill \end{array}$$

(5)

where $f_t$ is the forget gate activation, $W_f$ and $b_f$ are the weight matrix and bias, $h_{t-1}$ is the previous hidden state, and $x_t$ is the input at any cuttrent time.

$$\begin{array}{cc}\hfill i_t& =\sigma (W_i[h_{t-1},x_t]+b_i),{\tilde{C}}_t=\mathrm{ReLU}(W_C[h_{t-1},x_t]+b_C),\hfill \end{array}$$

(6)

where $i_t$ is the input activation function result, and ${\tilde{C}}_t$ represents the candidate values to add to the cell state.

$$\begin{array}{cc}\hfill C_t& =f_t\odot C_{t-1}+i_t\odot {\tilde{C}}_t\hfill \end{array}$$

(7)

where ⊙ represents element-wise multiplication, $f_t$ is the output of forget gate, $C_{t-1}$ is the previous cell state, $i_t$ is the new input and $C_t$ is the new cell state.

$$\begin{array}{c}\hfill o_t=\sigma (W_o[h_{t-1},x_t]+b_o),h_t=o_t\odot \mathrm{ReLU}\left(C_t\right)\end{array}$$

(8)

where $o_t$ is the activation of the output gate, and $h_t$ is the new hidden state.

Hidden states produced at each time step are the output of this layer and are collected as

$$S=[h_1^{\left(1\right)};\dots ;h_T^{\left(1\right)}]\in {\mathbb{R}}^{T\times d_h}.$$

(9)

To mitigate overfitting, a Dropout layer is applied immediately after, with a dropout rate p selected as 0.3 to produce the following output.

$$\tilde{S}=\mathrm{Dropout}(S;p)$$

(10)

2.2.2. Multihead Attention Layer

An Attention mechanism follows the initial LSTM layer to enhance feature representation by focusing on the most relevant parts of the sequence. The attention mechanism is a key element of transformer-based neural networks, enabling effective modeling of contextual dependencies. The multi-head attention variant extends this concept by projecting the input into multiple subspaces, one per head, to compute parallel attention scores that capture diverse contextual relationships within the sequence Vaswani et al. (2017). In the proposed model, six attention heads are used, with key vector dimensions equal to the LSTM output size, allowing the model to learn multiple attention patterns in parallel. For each attention head $i\in \{1,\dots ,N\}$, queries, keys, and values are computed as:

$$Q_i=\tilde{S}{W}_i^Q,K_i=\tilde{S}{W}_i^K,V_i=\tilde{S}{W}_i^V,$$

(11)

where $W_i^Q,W_i^K,W_i^V\in {\mathbb{R}}^{d_{\mathrm{model}}\times d_k}$ are learnable projection matrices and $d_k=d_v=d_{\mathrm{model}}/h$. The attention output for $hea{d}_i$ is computed through:

$${\mathrm{head}}_i=\mathrm{softmax}\left({\displaystyle \frac{Q_i{K}_i^{\top }}{\sqrt{d_k}}}\right)V_i,$$

(12)

and outputs of all heads are concatenated and linearly projected to form the output $\mathrm{A}$.

$$A=\mathrm{Concat}({\mathrm{head}}_1,\dots ,{\mathrm{head}}_N)W^O.$$

(13)

To integrate global contextual information with local temporal patterns, the attention outputs and outputs from the first LSTM’s hidden states are concatenated at each time step:

$$c_t=[h_t^{\left(1\right)};a_t].$$

(14)

2.2.3. Second LSTM Layer (LSTM₂)

The concatenated sequence $C=[c_1;\dots ;c_T]$ is then passed through a second LSTM layer with the same configurations as the ${\mathrm{LSTM}}_1$:

$$s_t={\mathrm{LSTM}}_2(c_t,s_{t-1},{\tilde{C}}_{t-1}).$$

(15)

This layer returns only the final hidden state $s^{*}=s_T$, capturing high-level temporal representations. A Dropout layer similar to the first LSTM layer is applied at this stage too for additional regularization.

$${\tilde{s}}^{*}=\mathrm{Dropout}(s^{*};p)$$

(16)

2.2.4. Output Layer

A final dense layer maps the final feature representation ${\tilde{s}}^{*}$ to the single scalar output $\widehat{y}$ which is then inverse-transformed to the original scale.

$$\widehat{y}=\sigma (W_d{\tilde{s}}^{*}+b_d),$$

(17)

where $W_d$ and $b_d$ are trainable parameters, and $\sigma (·)$ denotes the linear activation function to predict the final value.

2.3. Experimental Design

This section outlines the methodological framework employed to evaluate the proposed hybrid deep learning model for missing data imputation. The experiments extend our previously published univariate results (Laeeq et al. 2026) by investigating the model’s performance on multivariate sensor datasets. The complete workflow includes missing data simulation, model training with cross-validation, and comparative evaluations against multiple baseline and ablation variants.

To systematically assess model performance, missing data were simulated on the target variable of each dataset using Bernoulli random sampling at a rate of 30%. The missing indices were generated randomly across the time steps to mimic realistic stochastic data loss in sensor networks. The auxiliary features remained complete and were used to support temporal and cross-variable correlations during imputation.

Each dataset was partitioned into training and test sets with a ratio of 80:20 without shuffling to preserve the temporal sequence. While the training data was further partitioned into training and validation sets to train the proposed hybrid model using 5-fold cross-validation, to ensure robustness and generalizability. Mean Squarred Error (MSE) was used as the loss measure during the training process. The best-performing model, as determined by validation loss, was then evaluated on unseen test data under varying levels of missingness from 10% to 90%. Performance was quantified using the Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE), defined as:

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{({\widehat{y}}_i-y_i)}^2}$$

(18)

$$\mathrm{MAE}={\displaystyle \frac{1}{N}}\sum _{i=1}^N|{\widehat{y}}_i-y_i|$$

(19)

where ${\widehat{y}}_i$ and $y_i$ denote the imputed and true values, respectively. These metrics were computed across all missingness levels to evaluate the model’s robustness under increasing data degradation.

The hybrid model is implemented using Python with open source libraries of Pandas, Numpy and TensorFlow. Table 1 summarizes the configuration and training parameters for the proposed model, used in all experiments.

The proposed model was benchmarked against two widely used imputation baselines: K-Nearest Neighbour (KNN) imputation, which estimates missing values based on the mean of the k most similar instances. For our experiments, we set k to be five.

BRITS (Bidirectional Recurrent Imputation for Time Series), a deep learning-based approach using bidirectional recurrent dynamics to model temporal dependencies. We used the $sequence\_length$ for this model identical to $window\_size$ for our proposed hybrid model. The $learning\_rate$ was set to be 0.00001 and number of epochs used for each training fold was 50 for this model. These values are chosen to avoid zero training loss.

All comparative models were trained on the same datasets and evaluated using identical experimental conditions of 5-fold cross validation and performance metrics for fair comparison.

To investigate the contribution of each module within the hybrid model, three ablation variants were designed, as detailed in Table 2. Each variant omits or alters specific components to assess their impact on model performance. Each ablation model follows the same training and evaluation procedure as the full hybrid model. The resulting comparisons provide insights into the role of temporal sequencing and combining the outputs of LSTM and attention layer to b eused as input for the second LSTM layer, emphasizing the use of sandwich approach in the proposed model architecture.

3. Results

3.1. Multivariate Evaluation

Table 3 presents the MAE and RMSE values for four environmental datasets: SoilTemp, Meteorological, Air Quality NO₂/O₃, and Air Quality O₃ across missingness levels ranging from 10% to 90%. Across all datasets, the model demonstrates strong robustness at lower levels of missingness, with relatively small changes in MAE and RMSE up to around 70%. Beyond this point, both metrics begin to rise more noticeably, and the degradation becomes substantial.

For the Soil Temperature dataset, the model exhibits the most stable behaviour. As shown in Table 3 and Figure 5(a), both MAE and RMSE remain nearly constant up to 70% missingness. A substantial increase is observed only beyond 80%, indicating that soil temperature values are easier to reconstruct even when a large portion of data is missing. The Meteorological dataset displays similarly consistent performance at moderate missingness levels. Interestingly, MAE slightly improves up to 40%, however, Figure 5(b) show a clear increase in both metrics once missingness exceeds 60%, with a sharp rise at 80-90%. In the Air Quality NO₂/O₃ dataset, the model maintains stability up to approximately 40% missingness, but errors increase steadily after that. As illustrated in Figure 5(c), MAE nearly doubles between 10% and 90%, reflecting the higher variability and nonlinearity typically associated with air pollutant dynamics. Finally, the Air Quality O₃ dataset shows small absolute error values due to its lower numerical range, but the upward trend with increasing missingness is still evident (Figure 5(d)). While errors remain minimal up to 60%, both MAE and RMSE rise consistently beyond 60%, confirming that even low-magnitude datasets eventually exhibit degradation as missingness becomes severe.

Overall, the results show that the model performs reliably under moderate to high missingness conditions for all environmental datasets, but predictable degradation occurs at extreme levels of missing data.

3.2. Comparative Study

Table 4 reports the Mean Absolute Error (MAE) comparison and Table 5 shows the corresponding RMSE values for the Hybrid model, KNN and BRITS across four environmental datasets and missingness levels from 10% to 90%. Visual comparison is available in Figure 6 and Figure 7.

For the SoilTemp dataset, the proposed hybrid model consistently achieves the lowest MAE and RMSE across most missingness levels, demonstrating strong robustness to data loss. For example, at 30% missingness the Hybrid MAE is 0.0962 while KNN and BRITS record 0.4045 and 0.5733, respectively; the RMSE trend is also similar. These differences are clearly visible in Figure 6(a) and Figure 7(a).

The Meteorological dataset exhibits a different pattern, KNN outperforms the Hybrid model across the reported missingness levels. However, the RMSE results (Table 5 and Figure 7(b)) provides a more interesting picture. The Hybrid model attains slightly lower RMSE than KNN at low to moderate missingness, indicating that although KNN gives smaller average absolute errors, the Hybrid produces fewer very large errors in these missingness levels. At higher missingness (80-90%) KNN’s RMSE becomes lower than Hybrid’s, reflecting KNN’s relative resilience in reducing large residuals for this dataset at extreme data loss. These contrasting MAE/RMSE behaviours suggest that for Meteorological data KNN yields smaller average deviations while the Hybrid reduces large errors at low-to-moderate missingness.

For the NO₂/O₃ target, KNN initially achieves lower MAE than the Hybrid at low missingness; as presented in Table 4 and Figure 6(c). However, this advantage diminishes as missingness increases, and the Hybrid overtakes KNN from roughly 50% missingness onward. The RMSE series shows that Hybrid already attains lower RMSE than KNN from as early as 20% missingness, indicating that the Hybrid is better at limiting larger errors even when KNN has lower MAE at the very lowest missingness levels. In short, KNN is initially better on MAE for NO₂/O₃, but Hybrid becomes superior as missingness grows; RMSE shows the Hybrid reduces large residuals earlier than MAE suggests. For O₃, absolute errors are very small for all methods; nevertheless, the Hybrid maintains a small but consistent advantage in MAE across most missingness levels and RMSE values follow the same pattern.

3.3. Ablation Study

Ablation study was conducted to evaluate the performance of three ablation variants (A, B, C) as defined in Table 2 against the proposed hybrid baseline model across four datasets (SoilTemp, Meteorological, Air Quality NO₂O₃, and Air Quality O₃) using MAE and RMSE metrics under varying missing data percentages (10%-90%). The results of this study for all four datasets are summarised in Table 6 and Table 7 and a visual representation is created in Figure 8 and Figure 9.

For SoilTemp dataset, across all missingness levels, the Baseline achieves the lowest MAE and RMSE, confirming that the full hybrid structure is essential for accurate reconstruction. Ablation B performs the worst, with very large errors at all missingness levels, indicating that the LSTM component is critical for capturing temporal dependencies. Ablation A and Ablation C perform moderately but consistently worse than the Baseline. RMSE results in Table 7 and Figure 9(a) follow the same trend, reinforcing the dominance of the Baseline model.

A similar pattern is observed in the Meteorological dataset. Ablation B again yields extremely high error values across all missingness levels (Table 6 and Figure 8(b)), showing that temporal modeling is indispensable for meteorological dataset as well. Ablation A and Ablation C remain competitive with each other but remain consistently inferior to the Baseline. RMSE comparisons reflect the same ordering. The Baseline retains a clear advantage across all missingness levels.

For the NO₂/O₃ dataset, the Baseline again outperforms all ablated variants across both metrics. Ablation B remains the poorest performer (MAE increasing to 3.3611 at 90%), while Ablation A and Ablation C track more closely. Notably, Ablation C performs relatively closer to the Baseline at very low missingness (10%) but diverges at higher missingness. RMSE differences in Table 7 and Figure 9(c) show a sharper degradation for all ablations, with Ablation B exhibiting the steepest rise.

The O₃ dataset exhibits very small absolute error magnitudes, but the relative behaviour remains consistent with the other datasets. The Baseline consistently achieves the best MAE and RMSE values. At low missingness levels, Ablation A and Ablation C perform comparably to the Baseline, but Ablation B again shows noticeably higher errors. The RMSE curves in Figure 9(d) show similarly small but consistent gaps.

4. Discussion

The results of the multivariate evaluation of the proposed hybrid model demonstrate that it provides reliable imputation performance across diverse environmental datasets, maintaining stable MAE and RMSE values up to approximately 40-70% missingness. This robustness indicates that the architecture effectively captures temporal dependencies, allowing it to reconstruct missing patterns even when substantial portions of the sequence are absent. However, the sharp increase in error beyond 70-80% missingness highlights the limitations of the model when the available contextual information becomes insufficient. Air quality datasets, particularly NO₂/O₃, show more rapid deterioration due to their higher variability, whereas Soil Temperature remains stable for longer. These findings underscore that while the model is well-suited for real-world environmental monitoring scenarios with moderate to high data loss, extreme missingness still presents a significant challenge that future work should address through enhanced temporal reasoning or external information integration.

The comparative evaluation highlights the strong performance and robustness of the proposed Hybrid model relative to KNN and BRITS across the four environmental datasets. In the SoilTemp and O₃ datasets, the Hybrid model consistently yields the lowest errors, even under extreme missingness. These datasets exhibit smoother temporal variations, enabling the Hybrid model to effectively leverage both local feature extraction and sequence learning to reconstruct missing values with high stability.

In the Meteorological dataset, however, KNN demonstrates competitive behaviour at low missingness levels, marginally outperforming the Hybrid model in both MAE and RMSE. This result suggests that meteorological variables contain short-range correlations that KNN can exploit effectively through distance-based interpolation. Figure 10 shows the distribution comparison of $temp\_top$ and $tempAvg$, target features in SoilTemp and Meteorological daatsets respectively. $temp\_top$ may show skewness, reflecting localized fluctuations, while, $tempAvg$ distribution is smoother and less skewed, indicating stability from averaging. $temp\_top$ exhibits a wider range, suggesting more variability in sensor readings. On the other hand, $tempAvg$ is more concentratedand hence neigbouring values are favouring imputation. As missingness increases, KNN’s neighbourhood becomes sparse, and its performance deteriorates, allowing the Hybrid model to surpass it. This behaviour illustrates the advantage of Hybrid’s architecture in maintaining reconstruction accuracy under progressively degraded information availability.

A similar pattern emerges in the NO₂/O₃ dataset, where KNN initially achieves lower error metrics than the Hybrid model at 10-30% missingness. Pollutant concentrations often exhibit strong local covariation, explaining KNN’s effectiveness in early stages. However, as the missingness grows, the Hybrid model becomes the most accurate method, consistently outperforming both KNN and BRITS. The rapid performance decline of BRITS across all datasets and missingness levels indicates that, despite being designed for irregular time series, it struggles to maintain stability under the types of data gaps and nonlinear dynamics present in environmental monitoring.

Overall, the comparative study shows that the Hybrid model generalises more effectively than KNN and BRITS across diverse environmental datasets, particularly when missingness becomes substantial. Its ability to integrate local, sequential, and contextual information enables strong performance.

The ablation study confirms that each component of the Hybrid model; multihead attention and the two LSTM temporal modules contribute positively to reconstruction accuracy. The Baseline consistently achieves the lowest MAE and RMSE values and shows the most stable behaviour across all missingness levels. These results validate the necessity of the full hybrid design and highlight the complementary strengths of its constituent components. Across datasets and metrics, Ablation B (removal of the second LSTM) is the most detrimental modification. The consistently poor performance of Ablation B across SoilTemp, Meteorological and NO₂/O₃ datasets demonstrates that the post-attention sequence module (LSTM2) is critical for producing accurate reconstructions when attention outputs are to be re-integrated into a sequential representation.

Ablation A (removal of the first / pre-attention LSTM) and Ablation C (removal of the attention mechanism) both produce moderate degradations relative to the Baseline but do not cause the extreme failure observed for Ablation B. For NO₂/O₃, Ablation C is relatively close to the Baseline at low missingness (MAE 1.7326 vs. Baseline 1.6996 at 10%), but its error grows faster with increasing missingness (see Figure 8(c)), indicating that attention contributes increasingly with larger data gaps. The O₃ target (small-magnitude values) reflects the same qualitative pattern but with much smaller absolute errors: the Baseline maintains the best MAE/RMSE across missingness levels and Ablation B again produces the largest deterioration.

In summary, Removing the post-attention LSTM (Ablation B) causes the largest and most consistent performance drop across datasets and metrics, highlighting the importance of re-integrating attention-enhanced features into the sequential model. Omitting the pre-attention LSTM (Ablation A) or the attention module (Ablation C) leads to moderate but measurable degradations; attention becomes progressively more important as missingness increases. The full Baseline configuration consistently attains the lowest MAE and RMSE values, confirming that the combined sandwiched pipeline (LSTM1 + Attention + LSTM2) provides the most reliable reconstruction performance under varying degrees of data loss and varying types of data distributions.

Future research directions include hyperparameter optimisation and imputing a sequence of missing values as an output of this model. Understanding the inference of the model by analysing the attention weights to see what feature are important in making imputations and what time steps are more significant thatn others can be another future direction. Testing the model on datasets from other domains can be a useful way forward in evaluating the generalisability and robustness of the proposed model.

5. Conclusions

In this study, we presented a hybrid deep learning framework for robust missing data imputation in environmental monitoring datasets. By integrating LSTM with attention mechanism, in sandwiched interleaved manner, our model effectively captures both temporal and spatial dependencies, outperforming baseline and ablation variants in terms of MAE and RMSE. The results demonstrate that the proposed approach not only enhances imputation accuracy under varying levels of missingness but also maintains stability across different environmental datasets. These findings highlight the potential of hybrid architectures for improving data quality in environmental monitoring systems, which is critical for reliable decision-making and downstream tasks.