Chemometrics is essential for fully extracting the potential of contemporary instruments (22). In numerous investigations, it has been used to estimate different soil parameters using spectral data, such as near- or mid-infrared data, using various types of models, such as Partial Least Square Regression (2). As indicated by (23), the PLSR is the most employed learning modeling method. It has proven to be a good mathematical model for predicting various soil properties, such as soil organic carbon, total nitrogen, cation exchange capacity, clay, sand, pH, and other properties (3, 4). Predictions based on unprocessed soil sample laboratory spectra appear promising for quantitatively estimating soil chemical properties (24). However, the accuracy of this information can be affected by different factors, such as moisture, variations in soil roughness, and soil particle size. Therefore, data preprocessing methods that can enhance spectral information must be implemented (7, 24). Numerous preprocessing techniques, including standard normal variate, Savitzky–Golay smoothing, 1st derivative, and 2nd derivative, have been used in various studies and have demonstrated good prediction capabilities for soil properties (25, 26). The standard normal variate is an easy approach for spectral data normalization that works row by row, aimed at correcting for light scattering using the formula: SNV=xi − x¯iSi, where xi represents the signal of the observation, x¯i is its mean, and Si is its standard deviation (27).

The derivation of Savitzky–Golay pre-processing is one of the most extensively utilized drift-noise reduction strategies. The first and second are two orders, each with a special role. The second derivative is recommended to address the effect of the offset and linear trend of the baseline, while the first is intended only for a constant baseline effect (28). The Savitzky–Golay derivative uses the same procedure as the Savitzky–Golay smoothing described by (29), but with an added step. The derivative of this function is calculated after fitting the polynomial to a window around the moving point. This value was then used as the derivative estimate for the central point (29).

In areas where datasets with many variables are available, variable and feature selection has attracted considerable interest in the research community. The variable selection process has many advantages, such as enhancing data comprehension and visualization, lowering the need for measurement and storage, shortening the time needed for application and training, and reducing the impact of dimensionality to increase prediction accuracy (30). The use of some variable techniques, such as variable selection with interval-PLS (iPLS) and variable importance in projection (VIP), has proven to be more effective and increases the accuracy of the model (31, 32). The selectivity ratio is also beneficial for increasing predictive accuracy by concentrating on a small, highly relevant subset of variables (33).

Interval partial least squares (iPLS) is an approach for finding intervals of variables that are highly important for prediction. The algorithm splits the variables (predictors) into intervals and attempts to find the best combination of intervals that has the highest predictive accuracy. This interval selection can be performed using two methods: forward and backward. The forward method is when the intervals are included successively, whereas the backward method is when the intervals are excluded successively (34). In this study, forward iPLS was applied to the data after selecting the optimal preprocessing on R.

Variable importance in projection-partial least squares (VIP-PLS) is a multivariate statistical methodology employed to evaluate the significance of individual indicators in influencing an aggregate index. The primary objective of VIP-PLS is to establish a hierarchical ranking of indicators based on their respective degrees of importance to an aggregate index. This facilitates a clear understanding of the relative contributions of individuals (35). The VIP score is calculated using the following equation: VIPj=∑f=1FWjf2·SSYf·JSSYtotal·F, where w_jf represents the weight value for j variable and f component, SSY_f is the sum of squares of explained variance for the fth component, and J is the number of X variables. SSY_total is the total sum of squares explained by the dependent variable, and F is the total number of components (36).

The selectivity ratio is the ratio of the explained variance to the residual variance for each variable in a dataset, which is calculated after applying a target projection in a multivariate analysis. This ratio is based on the ability of a variable to discriminate between the different groups of samples. A high selectivity ratio indicates that a variable has a strong discriminative advantage (37). The selectivity ratio is defined as follows: SRi=SSi, explainedSSi, residual. The SS_i,residual represents the portion of the variance of variable i that is explained by the target projection component, that is, the part aligned with the PLS regression direction. In contrast, SS_i,residual denotes the remaining or unexplained variance of that variable after subtracting the target-projected part (the residual variance in the model). Therefore, the ratio quantifies the relevance of each variable’s variance to prediction versus residual noise, with higher SR values indicating greater predictive importance (36).

The results of predicting TN using PLS and different preprocessing methods showed that the Savitzky–Golay 1st derivative was the preprocessing method that yielded the best prediction, with an R²_cv of 0.88 and an RMSE_cv of 0.01%. Savitzky–Golay 2nd derivative, SNV coupled with SG1D, SG1D coupled with SNV, and SNV with SG2D also yielded a good prediction higher than the raw data prediction ranging between 0.85 and 0.88 for the R²_cv and 0.02% and 0.03% for the RMSE_cv, respectively. SNV and SG smoothing yielded good predictions but were lower than the raw data (R²_cv <0.8 and RMSE_cv = 0.03%).

All preprocessing techniques employed did not show any improvement and generated R²_cv below the R squared of the raw data, except Savitzky–Golay smoothing, which yielded similar results with an R²_cv of 0.88 and an RMSE_cv (0.28%) higher than the RMSE_cv of the raw data. The inconsistent performance observed for total carbon using raw MIR spectra likely reflects fundamental spectroscopic limitations rather than deficiencies in model calibration. Total carbon includes both organic and inorganic carbon fractions, which exhibit different MIR spectral behaviors. Organic carbon is mainly associated with broad and chemically informative absorption features, whereas inorganic carbonates show more specific bands that may overlap with organic signals (45, 46).

Consequently, raw spectra may already capture the dominant organic carbon information, whereas common preprocessing techniques can unintentionally suppress or distort these broad features. This likely explains why preprocessing did not improve and, in some cases, reduced the predictive performance compared with the raw spectra.

The combination of Savitzky–Golay’s 1st derivative and SNV yielded good predictions for total organic carbon (R²_cv = 0.81, RMSE_cv = 0.13%), with a better calibration performance of SNV applied to SG1D (SG1D + SNV). Savitzky–Golay’s 2nd derivative and SG2D applied to SNV resulted in the same predictive performance as the raw data models (R²_cv = 0.78, RMSE_cv = 0.14%). Conversely, SNV, SG smoothing, and SG2D + SNV generated acceptable results (0.76 < R²_cv < 0.77, 0.14% < RMSE_cv < 0.15%) but did not show any improvement in predictive performance (lower than the raw data).

The results of predicting clay content using PLS and various preprocessing methods highlighted that SG1D + SNV was the optimal preprocessing method with the highest R²_cv (0.86) and lowest RMSE_cv (7.04%). SNV, SG1D, SG2D, SNV + SG1D, SG2D + SNV, and SNV + SG2D also exhibited good predictive performance, with an R²_cv higher than and an RMSE_cv lower than the raw data. SG smoothing demonstrated an R²_cv similar to raw data (0.81) but higher RMSE_cv (8.25%).

With an R² value of 0.86 and an RMSE_cv of 3.72%, the application of Savitzky–Golay 2nd derivative after applying SNV (SNV + SG2D) produced the highest accuracy in the silt content prediction results utilizing PLS. With R²_cv values of 0.81 and a slight difference in the RMSE_cv (4.27%), Savitzky–Golay smoothing had a predictive performance similar to that of the raw data. The SG1D, SNV + SG1D, SG1D + SNV, and SG2D + SNV techniques also elevated R²_cv to 0.86 and lowered RMSE_cv to 3.8%. In contrast, SNV and SG2D did not improve the predictive accuracy of the model and yielded inferior R²_cv and RMSE_cv values of R²_cv = 0.62, RMSE_cv = 4.07%, and R²_cv = 0.6, RMSE_cv = 3.74%, respectively.

As shown in Table 2, the results of the sand content prediction using PLS and multiple preprocessing methods showed that the SNV + SG2D combination yielded the highest R²_cv and the lowest RMSE_cv, 0.86% and 8.87%, respectively. Savitzky–Golay’s 1st and 2nd derivatives, SNV + SG1D, SG1D + SNV, and SG2D + SNV also yielded good prediction performance and higher than the raw data with an RMSE_cv of 9.78%, 9.48%, 9.36%, 9.4%, and 9.28%, respectively, and an R²_cv of 0.83 for SG1D and 0.84 for the other preprocessing methods. SNV did not improve the prediction, resulting in similar results as the raw data (R²_cv = 0.79, RMSE_cv = 10.78%). In contrast to all the previous preprocessing methods, SG smoothing lowered R²_cv to 0.78 and elevated RMSE_cv to 11.04%, suggesting that smoothing is unsuitable for the prediction of sand in this context.

All the derivative preprocessing and their combination with SNV improved the prediction accuracy of the moisture content, with R²_cv ranging between 0.8 and 0.85 and RMSE_cv between 0.68% and 0.7%. The highest R²_cv (0.85) and lowest RMSE_cv (0.67%) were generated by SG1D + SNV, ranking it as the best preprocessing methods to predict MC in this context. SNV alone could not enhance the model’s predictive accuracy and delivered results similar to raw data (R²_cv = 0.8, RMSE_cv = 0.79%). In contrast, SG smoothing increased the RMSE_cv to 0.82% and decreased the R²_cv to 0.78.

The performance metrics showed that SG smoothing and SG2D + SNV decreased the predictive performance of the model compared to the raw data model, with R²_cv values of 0.61 and 0.58 and RMSE_cv values of 0.39 and 0.41, respectively. The SG1D + SNV also decreased the predictive performance, even though this preprocessing yielded a similar R²_cv (0.66) to the raw data but a higher RMSE_cv (0.37), which led to a higher predictive error. In contrast, SNV + SG2D generated a similar RMSE_cv but a higher R²_cv (0.67). The remaining preprocessing methods (SNV, SG1D, SG2D, and SNV + SG1D) improved the model performance by lowering the RMSE_cv and increasing the R²_cv, suggesting SNV + SG1D as the optimal preprocessing method for predicting pH (R²_cv = 0.73, RMSE_cv = 0.32).

Overall, none of the used preprocessing methods improved or at least maintained the predictive raw data model, except for SNV and SG smoothing, which showed better performance and nearly identical results (R²_cv = 0.62 and RMSE_cv = 24.29 mg/kg for SNV and R²_cv = 0.61 and RMSE_cv = 24.74 mg/kg). The low MIR prediction accuracy obtained for soil phosphorus agrees with previous studies reporting unreliable spectroscopic prediction of phosphorus due to its weak association with MIR-active soil components (47).

In contrast to phosphorus, all preprocessing improved the prediction accuracy of the CEC model, except for SG smoothing, which yielded results similar to those of the raw model’s metrics (R²_cv = 0.46, RMSE_cv = 6.3 meq/100 g). The results also show that the best preprocessing for CEC prediction using PLSR is SG2D + SNV, with the highest R²_cv (0.59) and the lowest RMSE_cv (5.52 meq/100 g). The moderate predictive performance obtained for CEC can be explained by the limitations of MIR spectroscopy. The cation exchange capacity is not spectrally active, and its prediction relies on indirect correlations with clay mineralogy and soil organic matter rather than direct absorption features (48). explicitly highlights this limitation, showing that accurate MIR-based CEC prediction is strongly dependent on the dataset size and modeling strategy.

Overall, the results presented in Tables 2, 3 suggest that the application of Savitzky-Golay (SG) preprocessing, either alone or in combination with standard normal variate (SNV), improves the model performance for most soil properties. However, the order of SNV (applying SNV to the derivative or applying the derivative to the SNV) plays a crucial role in decreasing or increasing the model performance. This combination yielded strong predictive performance for most soil parameters, indicating that it was the most effective preprocessing technique tested, except for TN, TC, and P₂O₅. SNV and Savitzky–Golay derivatives and their combination offered the best results for different soil properties and for different regression techniques (partial least squares, support vector regression, cubist, and convolutional neural networks), as proven in a similar study by (49). In a study by (50), the optimal preprocessing to develop a model for total nitrogen monitoring in situ and in-door for different depths and the optimal preprocessing combinations for the best-performing models (both indoor and in situ) consistently included Standard Normal Variate and derivative preprocessing. Other studies have found other optimal preprocessing combinations to predict soil properties; however, their optimal combination always includes derivative pretreatment, which highlights its role in improving model accuracy (26, 51).

The R²_cv and RMSE_cv of the models generated by the combination of optimal preprocessing and iPLS indicated a high level of prediction ability for most soil parameters, with an R²_cv over 90% for TN, TC, TOC, silt, sand, MC, and pH, and an R²_cv of 77% for CEC and 67% for P₂O₅, as illustrated in Figure 3. The resulting RMSE_cv was lower than that of the raw data and optimal preprocessing. Figure 4 illustrates that the selected intervals varied for different soil properties. Sand required the fewest intervals, utilizing 40 intervals. In contrast, P₂O₅ required the highest number of intervals (150). The remaining properties, TN, TC, TOC, Clay, silt, MC, pH, and CEC, required 45, 55, 110, 88, 60, 50, 66, and 70 intervals, respectively. These results indicate that the combination of iPLS and optimal preprocessing can improve the predictive accuracy of soil properties. Table 2 and Figure 3 show that the predictive performance improved for nearly all properties in this methodological integration. These values represent some of the highest accuracies among the preprocessing methods employed in this study. Demonstrating the robustness of the optimal preprocessing + iPLS method.

Improved predictive accuracy for the studied soil attributes was also achieved by applying the variable importance in projection variable selection method to the optimal preprocessing, except for clay, moisture content, and phosphorus, which did not show any improvement using this technique. As shown in Table 2, VIP showed excellent predictive accuracy for TC with an R²_cv of 0.91 and an RMSE_cv of 0.24%, and good predictive accuracy for TN, TOC, silt, and sand, with an R²_cv ranging between 0.82 and 0.89. In comparison, pH and CEC exhibited acceptable accuracy with an R²_cv of 0.75 and an RMSE_cv of 0.31 for pH and 4.32% for CEC.

Compared to VIP and iPLS, SR was not the best technique for improving the predictive accuracy of soil properties, as it did not show any improvement in total nitrogen, moisture content, clay, and pH. SR showed improved predictive accuracy of total organic carbon and silt with nearly similar R²_cv and RMSE_cv of VIP. The same case was observed for total carbon, but the R²_cv and RMSE_cv values were nearly similar to those of iPLS. Sand, P₂O₅, and CEC were improved using SR, with R²_cv values of 0.88, 0.63, and 0.73 and an RMSE_cv of 8.28%, 24.01%, and 4.43%, respectively.

Overall, the combination of optimal preprocessing with variable selection techniques has been shown to be effective in improving the accuracy of the developed models by increasing R²_cv and decreasing the models’ errors. However, when comparing the variable selection techniques, iPLS outperformed both VIP and SR. The iPLS method improved the metrics of all the properties. In contrast, VIP and SR could not improve certain properties, particularly the moisture content. In addition, the improvement rate of iPLS was higher than that of the other techniques (lowest RMSE_cv and highest R²_cv). In terms of model complexity, as illustrated in Figure 5, SR and VIP selected larger wavelength intervals, whereas iPLS selected shorter intervals, which means a less complex model. iPLS selected a lower number of intervals for some properties, such as total carbon, but for others, such as pH, a high number of intervals was needed to improve accuracy. Table 4 lists the intervals selected by iPLS and the number of selected wavenumbers for each property. Spectral interval selection using iPLS was conducted on the preprocessed data. As preprocessing methods modify the physical meaning of individual wavelengths, the selected intervals were used to optimize predictive performance and reduce model complexity rather than for direct chemical interpretation, which should be approached with caution when feature selection is applied (52). The proportion of selected MIR wavenumbers differed substantially among soil properties, highlighting the contrasting levels of spectral redundancy and information density. MC required the largest fraction of the spectrum (68.74%), indicating that predictive information for this property is distributed across the broad MIR regions. In contrast, properties such as P₂O₅ and clay were characterized by very limited spectral selections (3.15% and 5.50%, respectively), suggesting that their predictive signals are concentrated within narrow intervals. Intermediate selection levels were observed for organic-related properties (TN, TC, TOC), pH, and CEC, indicating an equilibrium between spectral coverage and reduced model complexity.

To assess whether the improvement in model performance after variable selection was statistically significant, a two-sided Wilcoxon signed-rank test was conducted on the cross-validated R² values before and after applying the iPLS (α = 0.05). The test revealed a significant increase in predictive accuracy following iPLS variable selection (V = 55, p = 0.0029), confirming that the observed improvement in the R²_cv values was statistically significant. This supports the effectiveness of iPLS in enhancing the PLSR model performance.

The application of iPLS in the present study has proven to boost model accuracy and reduce model complexity. In parallel research by (34), the iPLS performed comparably to other variable selection methods (principal variables, forward stepwise selection, and recursively weighted regression), but the iPLS demonstrated its power in optimizing the prediction model compared to the full spectrum Partial Least Squares; this technique offers a beneficial advantage which is the graphical output that enhances data interpretation. The results of this study demonstrated the combination of optimal preprocessing with variable selection improves the predictive accuracy of the model, but when choosing the optimal preprocessing and optimal variable selection techniques, a better predictive performance is achieved, lowering the error of the model. The studied techniques (iPLS, VIP, and SR) are not commonly used in soil spectroscopy, which opens the window for further research and application and can guide future research in developing robustness and testing other preprocessing and variable selection techniques.