Spectral Method for Liming Recommendation in Oxisol Based on the Prediction of Chemical Characteristics Using Interval Partial Least Squares Regression

Abstract

Thousands of chemical analyses are carried out annually with the aim of recommending soil correction; however, these analyses are expensive, destructive, time-consuming, and can be harmful to the environment. As an alternative to conventional analysis methods, diffuse reflectance spectroscopy has been proposed as an option for evaluating the chemical characteristics of soil. The selection of variables has also emerged as an alternative to improve the performance of PLSR (partial least squares regression), as it decreases the root mean square error (RMSE) and increases the accuracy of the models. However, few studies have used a previous selection of variables for the construction of PLSR models to estimate the chemical characteristics of soil. In this context, the hypothesis in this study was that it is possible to calculate the liming recommendation in Oxisol based on the chemical characteristics estimated by PLSR, with a previous selection of variables using iPLS (Interval PLS). The objective was to calculate the need for liming based on chemical characteristics estimated via iPLS selection and PLSR modeling of specific wavelengths of soil reflectance. The experimental area was treated with different application rates of limestone, with and without incorporation, and phosphogypsum was applied in additional treatments. Soil assessments were carried out 5, 12, 24, and 36 months after the application of the treatments, using six layers: 0.00–0.05, 0.05–0.10, 0.10–0.20, 0.20–0.30, 0.30–0.40 and 0.40–0.60 m. Samples were subjected to conventional laboratory analyses, and spectral readings (400–2500 nm) were obtained with a spectroradiometer. The spectral curves were subjected to the iPLS variable selection method to generate PLSR models of the chemical characteristics used to calculate the liming recommendation. The chemical characteristics of the soil, such as Ca²⁺, sum of bases (SB), effective cation exchange capacity (CTCe), cation exchange capacity (CTC), and base saturation (BS), could be estimated, with values of R² ranging from 0.83 to 0.92 in the calibration and validation steps, and from 0.84 to 0.90 for the prediction step (in the fourth assessment). The liming recommendation calculated based on the chemical characteristics predicted from the PLSR models showed a strong correlation (r > 0.86) with the liming recommendation calculated by conventional laboratory techniques. The fourth soil assessment yielded the best correlation coefficient (r = 0.95).

1. Introduction

A large proportion of the world’s agriculture in tropical and subtropical environments is carried out on oxisols. Most Brazilian oxisols impose limitations on the establishment and development of the main crops due to the acidity of these soils, the presence of toxic elements (such as Al³⁺ and Mn²⁺), and low levels of basic cations (such as Ca²⁺ and Mg²⁺), which restrict the development of the plant root system and reduce their ability to absorb water and nutrients, especially in regions or periods of low water availability [1,2,3,4,5]. The adverse effects of soil acidity on plants occur both in the surface layer (0.0–0.2 m) and in the subsurface (depth > 0.2 m) of acidic soils. Thus, to promote greater efficiency of water and nutrient absorption by plants, and hence to enlarge yields, the surface acidity must be reduced and the chemical environment in the subsurface improved.

The application of limestone (liming) aims to correct acidity (increasing the pH), which reduces the toxic effects of high concentrations of Al³⁺ and Mn²⁺ and increases the levels of Ca²⁺ and Mg²⁺ in the soil. Liming also increases the effective CTC and decreases the leaching of bases, and is, therefore, considered one of the practices that most contributes to the increase of fertilizer efficiency, and as one of the pillars for obtaining higher yield and increased agricultural profitability and sustainability [2,6,7,8]. Since liming mainly improves the characteristics of the soil surface, especially in a no-till system, the application of phosphogypsum is an alternative to surface liming, as this improves the root environment in the subsurface of the soil and, under certain conditions, makes limestone incorporation unnecessary [9,10].

In 2020, Brazilian agribusiness used about 45 million tons of limestone, representing an increase in consumption of almost 60% over the previous 20 years and 45% over the previous 10 years [11]. Despite this increase in limestone consumption in Brazil, it is estimated that the ideal consumption would be around 70 million tons per year [11]. This increase in limestone consumption is directly related to soil analysis. The characteristics used to estimate liming recommendations are obtained through conventional laboratory analyses; although these are reliable techniques, they are expensive, destructive, time-consuming, and require reagents which can cause harm to the environment. Furthermore, the sample density is limited by the costs of laboratory tests and operational logistics [12,13,14,15].

As an alternative to conventional soil analysis methods, diffuse reflectance spectroscopy has been used to quickly, non-destructively, and cost-effectively assess soil characteristics [16,17]. The use of this technique to determine the chemical, physical and mineralogical properties of soil has shown favorable results in the regions of the visible (Vis) (350–750 nm), near-infrared (NIR) (750–1100 nm), and shortwave infrared (SWIR) (1100–2500 nm) ranges. The availability of such a large set of spectral information allows researchers to estimate multiple soil characteristics from a single reflectance spectrum [16,18,19,20,21].

Several studies have demonstrated the potential of diffuse reflectance spectroscopy to estimate exchangeable soil ions, such as Ca²⁺, Mg²⁺, K⁺, P, and SO₄²⁻, and pH [22,23,24,25,26,27,28]. This is made possible by the interaction between electromagnetic radiation and soil particles, which causes chemical bonds to vibrate at specific frequencies. This specificity allows for the evaluation of different soil characteristics once the reflectance is known [29,30].

Multivariate statistical techniques can be used to build prediction models for specific soil characteristics based on soil reflectance. However, these techniques need to be robust and accurate, given the large volumes of data. Among the different statistical methods, partial least squares regression (PLSR) is one of the most used [27,31,32,33,34].

Initially, it was supposed that the inclusion of more spectral regions would improve the accuracy of PLSR models, but several studies have found that the use of the full spectrum can cause overfitting. The selection of variables technique has emerged as an alternative to minimize this problem, as it improves the performance of regression, decreases the root mean square error (RMSE), and increases the accuracy of the models [35,36,37,38]. In this context, an interactive algorithm called Interval PLS (iPLS) was developed to generate local PLS models in subintervals of the spectral curve of interest. iPLS selects variables with relevant information in different spectral subdivisions, thus focusing on critical spectral regions and removing interference from other regions [39].

The iPLS algorithm selects a relevant subset of intervals (a subset of wavelengths) for each of the soil characteristics individually, which results in superior models compared to the use of the entire Vis-NIR spectrum, and hence decreases the complexity of the models by removing redundant information from the spectra [40]. When iPLS has split the spectra into different ranges, a PLSR validation model is applied to each of them. The interval that gives the lowest prediction error is modeled again by adding the remaining intervals, one by one, until the error no longer decreases [41].

Although many studies have reported the use of the Vis-NIR-SWIR range to predict soil characteristics, studies that have used a previous selection of variables for constructing PLSR models to estimate the necessary chemical characteristics are still rare in terms of calculating liming recommendations. These could make it possible in the future to develop cheaper sensors (or even cameras using specific wavelengths) which would not need to operate in the full Vis-NIR-SWIR range.

Based on current progress in this area, the hypothesis of this study is that it is possible to calculate liming recommendations for oxisols based on the chemical characteristics estimated by PLSR, with previous application of variable selection by iPLS. Our objective was therefore to calculate the need for liming based on the chemical characteristics estimated by iPLS selection and PLSR modeling of specific wavelengths of soil reflectance.

2. Materials and Methods

2.1. Experimental Area

The experiment (Figure 1) was carried out at the experimental farm of the Cooperativa Agroindustrial Mourãoense (COAMO), located in Campo Mourão, State of Paraná, Brazil (24°05′28″ S, 52°21′31″ W; altitude 605 m) and was evaluated between 2016 and 2019.

According to the Köppen classification, the climate in this region is classified as (Cfa), a humid subtropical climate with hot summers. The annual average temperature is 20–21 °C (in the coldest quarter, the average temperature is 16–17 °C, and in the warmest, 27–28 °C), and the annual precipitation ranges from 1600 to 1800 mm [42]. The weather data (minimum and maximum temperatures and rainfall) for the experimental area between June 2016 and May 2019 are shown in Figure 2. During this period, the accumulated precipitation was approximately 5484 mm.

Before the experiment, the area had been managed with a no-till approach since 2008, with the following crops: oats (2008), corn (2008/2009), oats (2009), soybeans (2009/2010), wheat (2010), soybeans (2010)/2011) and corn (2011/2012). Since 2012, alternating crops of soybeans and wheat have been cultivated.

The soil within the experimental area was classified as a Latossolo Vermelho distroférrico típico [43], corresponding to a ferralsol [44] and an oxisol [45] of very clayey texture (740 g kg⁻¹ clay). This soil analysis was carried out on 18 May 2016 (

Table 1. Chemical characterization of the soil before treatment in 2016 and limestone application rates used.

Treatment	Layer	Ca2+	Mg2+	K+	P	pH	H+Al	CEC	m	BS	Lime Rate
	m	cmolc dm−3			mg dm−3	CaCl2	cmolc dm−3		%		Mg ha−1
IL BS% 44	0.00–0.20	2.86	0.58	0.39	22.81	4.97	5.35	9.18	2.05	41.72	-
	0.20–0.40	2.04	0.39	0.28	7.91	4.75	5.96	8.67	4.91	31.26	-
IL BS% 60	0.00–0.20	2.73	0.63	0.23	23.55	4.79	5.55	9.14	4.01	39.28	2.3
	0.20–0.40	2.14	0.51	0.20	7.14	4.82	5.15	8.00	4.04	35.63	-
IL BS% 70	0.00–0.20	3.24	0.88	0.29	21.37	5.08	4.61	9.02	0.00	48.89	2.3
	0.20–0.40	2.78	0.70	0.25	10.26	5.23	4.44	8.17	0.00	45.65	-
IL BS% 90	0.00–0.20	3.03	0.99	0.30	20.88	5.07	5.55	9.87	0.00	43.77	5.6
	0.20–0.40	2.30	0.74	0.25	7.80	4.84	4.96	8.25	2.95	39.88	-
SL BS% 44	0.00–0.20	3.52	0.74	0.26	21.38	5.00	5.15	9.67	0.88	46.74	-
	0.20–0.40	2.12	0.41	0.20	5.08	4.74	5.15	7.88	6.51	34.64	-
SL BS% 60	0.00–0.20	3.47	0.9	0.26	30.06	5.16	4.96	9.59	0.00	48.28	1.4
	0.20–0.40	2.14	0.56	0.19	8.02	4.92	4.78	7.67	2.69	37.68	-
SL BS% 70	0.00–0.20	3.11	1.00	0.24	13.73	5.24	4.28	8.63	0.00	50.41	2.1
	0.20–0.40	2.55	0.77	0.22	6.12	5.14	4.44	7.98	0.00	44.36	-
SL BS% 90	0.00–0.20	3.57	1.46	0.21	22.44	5.41	4.28	9.52	0.00	55.04	4.1
	0.20–0.40	2.04	0.73	0.18	8.67	4.80	5.15	8.1	4.84	36.42	-
BS60+G50 *	0.00–0.20	4.11	0.86	0.26	23.18	5.29	4.28	9.51	0.00	54.99	0.6
	0.20–0.40	2.40	0.45	0.20	5.68	4.97	4.61	7.66	0.97	39.82	-
BS70+G50 *	0.00–0.20	3.09	0.70	0.19	26.74	5.00	4.61	8.59	0.75	46.33	2.5
	0.20–0.40	1.71	0.34	0.15	5.39	4.60	5.35	7.55	0.12	29.14	-
BS70+G100 **	0.00–0.20	3.16	0.65	0.19	34.34	5.06	4.96	8.96	0.00	44.64	2.8
	0.20–0.40	2.01	0.31	0.14	7.14	4.75	5.35	7.81	8.21	31.50	-

* represents an application rate of 3.71 and ** represents an application rate of 7.42 Mg ha⁻¹ phosphogypsum (G) applied in 2016. The phosphogypsum rate (GR) was calculated as a function of the clay soil content according to the formula GR = 50 × clay (%), as proposed by Sousa et al. [47]. The highest GR application rate was determined based on twice the recommended application rate [GR = 100 × clay (%)]. IL indicates incorporated liming; SL indicates surface liming; G indicates phosphogypsum; Ca²⁺, Mg²⁺ and Al³⁺ contents were extracted using KCl 1.0 mol L⁻¹; K⁺ and P contents were extracted using Mehlich-1; pH CaCl₂ (0.01 mol L⁻¹) at a soil:solution ratio of 1:2.5; H+Al was extracted using the SMP method; CEC is the cation exchange capacity at pH 7 (CEC = Ca²⁺ + Mg²⁺ + K⁺ + H+Al); m represents the aluminum saturation [m=(Al³⁺/Ca²⁺+Mg²⁺+K⁺+Al³⁺) × 100]; BS is the base saturation [BS = (Ca²⁺ + Mg²⁺ + K⁺/CEC)].

).

2.2. Implementation of the Experiment

The experimental design used was that of complete blocks with treatments distributed at random in a 2 × 4 crossed factorial scheme with three additional treatments and four blocks (giving a total of 11 treatments). These treatments included two types of limestone scale application (superficial and incorporated) and four base saturation levels (V%), corresponding to 44% (without correction), 60%, 70%, and 90%. Three additional treatments were applied that involved the use of phosphogypsum: V% 60 + 3.71 Mg ha⁻¹ of phosphogypsum (60G50), V% 70 + 3.71 Mg ha⁻¹ of phosphogypsum (70G50), and V% 70 + 7.42 Mg ha⁻¹ of phosphogypsum (70G100). Each experimental plot was 12 m in length and 7 m in width, giving a total of 44 experimental plots of 84 m².

The limestone application rates were calculated based on the elevation of base saturation recommended for the state of Paraná [5]. The application rates of phosphogypsum were calculated using the formula proposed by Sousa et al. [46], where the clay content was multiplied by 50 to obtain the recommended application rate and by 100 for twice the recommended application rate. The applied dolomitic limestone had a relative total neutralizing power (RTNP) of 80%. The treatments were applied on 10 June 2016, and the limestone application rates used are given in Table 1.

Limestone was incorporated using a moldboard plow with an effective depth of incorporation of 0.20 m, followed by harrowing with a plow harrow containing 20 28” discs (Ø 0.71 m) and then light harrowing with a harrow containing 42 20” discs (Ø 0.51 m) to level the ground.

2.3. Soil Assessment and Chemical Analysis

Four soil assessments were performed at 5, 12, 24, and 36 months after the treatments were applied on 10 June 2016. The sampled profile was stratified into six layers: 0.00–0.05, 0.05–0.10, and 0.10–0.20 m (with a cutting shovel) and 0.20–0.30, 0.30–0.40, and 0.40–0.60 (with the aid of a Dutch auger). Two subsamples were collected within each experimental plot to obtain a composite sample. On each soil assessment date (5, 12, 24, and 36 months after applying the treatments), 264 samples were collected, giving a total of 1056 soil samples.

The soil samples were oven-dried at 40 °C, and passed through a 2 mm sieve [5]. Exchangeable Ca²⁺, Mg²⁺, and Al³⁺ were extracted using a KCl 1 mol L⁻¹ (1:10 v/v soil/solution). Exchangeable Ca²⁺ and Mg²⁺ contents were determined by atomic absorption spectrophotometry with an air-acetylene flame. Exchangeable Al³⁺ was determined by titration (NaOH 0.0125 mol L ⁻¹). Exchangeable K⁺ was extracted using Mehlich-1 solution (H₂SO₄ 0.0125 mol L⁻¹ and HCl 0.05 mol L⁻¹ at a ratio of 1:10 (v/v) soil/solution) and determined by flame emission spectroscopy. Phosphorus (P) was extracted by Mehlich-1 (H₂SO₄ 0.0125 mol L⁻¹ and HCl 0.05 mol L⁻¹ at a ratio of 1:10 (v/v) soil/solution) and was determined by UV-Vis spectrophotometry. Extraction of S-SO₄²⁻ was performed using calcium phosphate (500 mg L⁻¹ P) in acetic acid 2.0 mol L⁻¹ (1:2.5 v/v), and the concentration was determined using a UV/Vis spectrophotometer. Soil pH was determined in a CaCl₂ 0.01 mol L⁻¹ suspension (at a soil/solution ratio of 1:2.5 v/v). Following the procedures described by Pavan et al. [47], the potential acidity was determined using SMP buffer solution. Except for the potential acidity, the other methodologies were those of Silva et al. [48].

The sum of bases (SB) was calculated by adding the Ca²⁺, Mg²⁺, and K⁺ contents. In naturally acidic and intensely leached soils, the values for Na⁺ and NH₄⁺ are very low (practically negligible) and these were therefore excluded from the calculation of BS. The effective cation-exchange capacity (CECe = SB + Al³⁺), cation exchange capacity at pH 7 (CEC = SB + H+Al), base saturation (BS = (SB/CEC) × 100) and aluminum saturation (m% = (Al³⁺/CECe) × 100) were also calculated. The chemical characteristics of the soil were subjected to a descriptive statistical analysis.

2.4. Acquisition and Processing of Spectral Data

Soil samples were placed in Petri dishes 9 cm in diameter, and spectral readings were acquired in the spectral range 350–2500 nm (Vis-NIR) with a Fieldspec 3 Jr spectroradiometer (ASD Inc, Boulder, CO, USA), with a spectral resolution of 3 nm between 350 and 700 nm and 10 nm between 700 nm and 2500 nm.

The geometry used to collect spectral data was adopted from Demattê et al. [49] and Nanni et al. [50]: a white Spectralon plate was used as a reflectance standard with 100% reflectance, as per Labsphere Reflectance Calibration Laboratory (LABSPHERE Inc., North Sutton, NH, USA); and a 650 W halogen lamp was used as a light source, positioned 0.35 m from the sample and with a zenith angle of 30°. The optical sensor was placed in an upright position 0.08 m away from the sample. The soil reflectance was measured at the center of the sample, as shown in Figure 3.

For each soil sample, three readings were performed by rotating the Petri dish through 120° [50,51], to avoid any bias due to the arrangement of particles. The reflectance factor was obtained, which expresses the ratio between the radiant flux of the soil sample and the radiant flux reflected by the reference standard under the same geometry and lighting conditions [52].

The collected spectra were corrected using the “splice correction” routine of the ViewSpecPro 6.0 software package (ASD Inc, Boulder, CO, USA). This routine aims to eliminate steps (or offsets) in the original spectra caused by phase differences between the measurements of the three detectors used with the spectroradiometer: VNIR (350–1000 nm), SWIR 1 (1000–1800 nm), and SWIR 2 (1800–2500 nm). Parts of the spectral curve generated for each sample were removed due to noise, mainly at the beginning and end of the spectral curve. The final spectrum band used for analysis was 400–2400 nm, with a total of 2001 bands.

Several techniques were tested for the preprocessing of the curves, to eliminate or smooth the curves and reduce interference related to the sensor, in order to give better calibration and prediction results. This preprocessing is usually performed in chemometrics with spectral data to reduce optical interference in spectra unrelated to the chemical composition of the sample [22,28,50]. Unscrambler X 10.4.1 software (Camo Software, Oslo, Norway) was used for this step.

Of the preprocessing methods tested at this stage, the Savitzky–Golay approach [53] resulted in the best performance of the model, with a first-order polynomial and a five-point segment about the central point. This technique is a digital polynomial filter in which each value in the signal series is replaced with a new value obtained from a polynomial fit to successive subsets of adjacent data points. The fit is performed using linear least squares (with 2n + 1 neighboring points), where n can be equal to or greater than the order of the polynomial. The more neighbors used in the averaging process, the smoother the signal becomes. Least squares smoothing suppresses noise while maintaining signal information [54].

2.5. Discriminant Analysis

After smoothing the curves, a discriminant analysis was performed on the spectral responses from the soil samples with 11 treatments, from the six sampled layers (0.00–0.05, 0.05–0.10, 0.10–0.20, 0.20–0.30, 0.30–0.40, and 0.40–0.60 m).

SAS software (SAS Institute Inc., Cary, NC, USA) was used for this discriminant analysis, using the STEPDISC and DISCRIM procedures. The STEPDISC procedure was executed by the FORWARD method at 5% probability, and began with no variables in the model. PROC STEPDISC inserted the variable that contributed the most to the model’s discriminatory power at each step, as measured by Wilks’ lambda. The process ended when none of the unselected variables met the entry criteria [55].

The variables that best explained the differences among the treatments and the sampled layers were selected to compose the discriminant model. This procedure generated discriminant models through linear combinations of predictor variables (wavelengths) [55].

Using SAS software, a simulation was performed to test the quality of the models obtained. In this simulation, 70% of the spectral curves were randomly selected and used to construct a new discriminant model, which was later tested using the remaining 30% of the spectral curves. This simulation process was repeated 50 consecutive times [50,56,57].

2.6. Partial Least Squares Regression

After the discriminant analysis, variables were selected based on the iPLS method, using PLS Toolbox 8.9.2 (Eigenvector Research Inc., Wenatchee, WA, USA), a plugin for MatLab software (MathWorks Inc., Natick, MA, USA). PLS Toolbox was operated in “forward” mode, where intervals are successively included in the analysis. The program was set to automatic mode, in which the iPLS algorithm continued to add (or remove) gaps until the RMSE no longer improved. The size of the intervals was defined as “one” (1 = one variable per interval), meaning that the number of selected intervals was equal to the number of selected variables [58].

PLSR was performed after selecting the wavelengths, using Unscrambler X 10.4.1 software (Camo Software, Oslo, Norway) with the NIPALS algorithm [59]. PLSR was conducted in three steps, consisting of calibration, validation, and prediction. The model quality parameters obtained in the validation steps indicated whether the models were suitable for use in predicting the chemical characteristics of the soil.

The database was randomized using Unscrambler X 10.4.1 software (Camo Software, Oslo, Norway), and was divided into two groups: the calibration and validation set (containing 75% of the samples) and the prediction set (containing the remaining 25% of the samples). PLSR models were created to quantify the following chemical characteristics of the soil: Ca²⁺, Mg²⁺, K⁺, P, SO₄²⁻, pH, Al³⁺, H+Al, SB, CECe, CEC, m%, and BS.

In this study, we used the following statistical parameters as indicators of the goodness of fit of the regression models: the coefficient of determination (R²), root mean square error (RMSE), systematic error (BIAS), relative percentage deviation (RPD), and the ratio of performance to interquartile range (RPIQ).

R² values were classified according to Wenjun et al. [60], Williams [61], Viscarra Rossel et al. [34], and Saeys et al. [62], where: R² < 0.50 represented models with poor predictions, which were only able to distinguish between high and low values; R² between 0.50 and 0.65 indicated models with moderate predictions that could be used for evaluation and correlation; R² between 0.65 and 0.80 represented good prediction models that allowed for quantitative predictions; R² between 0.80 and 0.90 indicated very good quantitative prediction models; and R² ≥ 0.90 represented excellent quantitative prediction models.

The RPD values were interpreted following Chang et al. [63] and the results were divided into three categories: highly reliable models with RPD > 2 (which can be used for quantitative analysis and making predictions); reliable models with 1.4 < RPD < 2; and unreliable models for predictions with RPD < 1.4.

Four categories of models were defined based on the RPIQ values: an RPIQ score of between 2.02 and 2.70 indicates a poor model, which can only distinguish between high and low; an RPIQ of between 2.70 and 3.37 represents a regular model, in which quantitative predictions are possible; an RPIQ of between 3.37 and 4.05 represents a good quantitative model; and RPIQ > 4.05 implies an excellent model [64].

2.7. Spectral-Based Method for Liming Recommendation

The chemical characteristics predicted by the PLSR calibration models were used for liming recommendation (LR) using the base saturation method, as follows:

$$LR\left(Mg \right)=\times CEC$$

where BS₂ is the desired base saturation (70%), BS₁ is the current base saturation; RPTN is the relative power of total neutralization (100%); and CEC is the cation exchange capacity at pH 7 (cmol_c dm⁻³), according to SBCS/NEPAR [5].

A correlation analysis was then performed between the estimated liming recommendation using the characteristics obtained by the PLSR models (based only on the prediction dataset), and the calculated liming recommendation, using the characteristics obtained by conventional laboratory techniques.

3. Results and Discussion

3.1. Soil Chemical Characteristics

Laboratory analyses of the soil samples from the six different layers (0.00–0.05, 0.05–0.10, 0.10–0.20, 0.20–0.30, 0.30–0.40, and 0.40–0.60 m) indicated a large variability in the results of the four soil assessments (

Table 2. Descriptive analysis of the chemical characteristics of the soil over the four soil assessments.

Statistics *	Ca2+	Mg2+	K+	P	S-SO42−	pH	Al3+	H+Al	SB	CECe	CEC	m	BS
	cmolc dm−3			mg dm−3		CaCl2	cmolc dm−3					%
First assessment (five months)
Minimum	0.23	0.14	0.06	0.64	1.65	3.53	0.00	2.54	0.56	0.94	5.04	0.00	6.93
Maximum	6.80	3.22	0.61	133.31	168.46	6.58	0.91	7.76	8.59	8.61	11.77	58.42	75.65
Mean	1.77	0.73	0.22	28.69	17.64	5.08	0.18	4.75	2.72	2.90	7.47	9.79	34.21
Std. dev.	1.23	0.52	0.12	33.24	28.39	0.49	0.17	0.92	1.72	1.63	1.37	10.72	15.92
CV %	69.38	71.75	57.57	115.87	160.89	9.56	93.89	19.39	63.49	56.34	18.32	109.51	46.54
Second assessment (12 months)
Minimum	0.31	0.17	0.01	0.25	0.00	4.32	0.06	2.54	0.64	1.25	5.18	0.73	8.36
Maximum	7.93	3.40	0.87	101.40	112.80	6.09	1.25	7.15	10.00	10.10	12.74	66.22	78.80
Mean	2.31	0.84	0.23	15.08	21.48	4.99	0.28	4.77	3.37	3.65	8.14	11.69	38.19
Std. dev.	1.56	0.71	0.20	17.89	21.04	0.39	0.22	0.99	2.30	2.18	1.76	11.22	17.62
CV %	67.49	84.97	90.26	118.68	97.96	7.86	78.44	20.79	68.22	59.87	21.64	95.91	46.15
Third assessment (24 months)
Minimum	0.36	0.13	0.05	0.36	1.53	4.22	0.02	2.86	0.64	1.40	5.33	0.23	8.35
Maximum	8.59	3.39	0.71	107.42	258.18	6.21	1.56	10.07	11.49	11.63	13.52	70.86	78.20
Mean	2.29	0.85	0.23	16.64	19.67	4.97	0.33	5.10	3.37	3.70	8.47	13.68	36.54
Std. dev.	1.63	0.75	0.15	17.16	18.79	0.43	0.27	1.11	2.38	2.25	2.21	13.03	17.00
CV %	71.17	87.80	65.16	103.11	95.54	8.60	81.87	21.85	70.60	60.67	26.07	95.23	46.52
Fourth assessment (36 months)
Minimum	0.19	0.12	0.02	0.76	3.70	4.32	0.04	2.60	0.57	1.27	4.46	0.51	9.58
Maximum	5.46	3.54	0.77	77.90	38.53	6.10	1.22	7.15	8.56	8.64	12.20	68.15	75.08
Mean	1.81	0.82	0.22	13.83	10.22	4.94	0.30	4.08	2.86	3.16	6.95	13.52	38.21
Std. dev.	1.14	0.65	0.17	15.22	5.93	0.37	0.22	0.69	1.81	1.69	1.65	12.38	15.38
CV %	62.82	78.23	76.93	110.04	58.02	7.44	75.00	16.83	63.22	53.46	23.69	91.58	40.24

* The data follow the normal Kolmogorov–Smirnov distribution (p-value ≤ 0.01). Ca²⁺, Mg²⁺ and Al³⁺ contents were extracted using KCl 1.0 mol L⁻¹; K⁺ and P contents were extracted using Mehlich-1; pH CaCl₂ (0.01 mol L⁻¹) with a soil:solution ratio of 1:2.5; H+Al was extracted using the SMP method; CECe is the effective cation exchange capacity; CEC is the cation exchange capacity at pH 7; m is the aluminum saturation; BS is the base saturation; CV is the coefficient of variance; std. dev. is the standard deviation.

). This was as expected, due to the large variability of the chemical characteristics in the soil profile, particularly when compared to the surface layers (which had higher SB, CECe, and BS content) and subsurface layers (with lower SB, CECe, and BS content). The smallest coefficients of variance were obtained for H+Al and CEC, as this oxisol showed only slight mineralogical and texture variations along the profile. The highest values of the coefficient of variance (CV) were observed for P (118.68% in the second soil assessment), S-SO₄²⁻ (160.89% in the first soil assessment), and Al³⁺ (93.89% in the first soil assessment). The lowest CVs were observed for pH (7.44% in the fourth soil assessment), H+Al (16.83% in the fourth soil assessment), and CEC (18.32% in the first soil assessment).

The minimum and maximum levels of Ca²⁺ increased after the application of the treatments (5 to 24 months after application). In the fourth soil assessment (36 months after application), the Ca²⁺ content was found to have decreased. The Mg²⁺ content did not vary significantly among the different soil assessments. The SB increased after the treatments were applied, but the values in the fourth soil assessment were similar to those in the first soil assessment. This can be explained by the extraction of bases by the cultures and their leaching. Similar behavior was observed for CECe, CEC, and BS (Table 2).

3.2. Spectral Modeling: Discriminant Analysis

The selection of variables for the discriminant analysis performed among treatments, using PROC STEPDISC at 5% probability, did not select any variable. Thus, the discriminant analysis between treatments was interrupted. The discriminant analysis performed to discriminate the different soil layers sampled was completed. The accuracy of the discriminant analysis of the four soil assessments is illustrated in Figure 4.

It was observed that the 0.00–0.05 m layer had the highest accuracy in terms of the classification of samples among soil assessments, ranging from 96.27–100% in the calibration step and 92.46–100% in the validation phase. The discriminant model of the 0.00–0.05 m layer in the fourth soil assessment reached 100% accuracy in both the calibration and validation steps. This separation of the 0.00–0.05 m layer was probably due to the higher content of organic matter and SB, which are usually higher in the first centimeters of the soil profile, thus distinguishing it from the other layers.

The 0.40–0.60 m layer also obtained a high accuracy percentage, ranging from 84.58–90.63% in the calibration step and 74.06–86.49% in the validation step for the different soil assessments. This difference was probably due to the lower content of organic matter and SB in the 0.40–0.60 m layers.

The intermediate soil layers between 0.05 and 0.40 m had an accuracy ranging from 65.59–95.22% in the calibration step and 51.74–85.14% in the validation step. The lower accuracy of these discriminant models is due to the larger variation observed among the spectral curves of these layers.

It was observed that the models had a high percentage accuracy, indicating that this statistical technique is valid and sensitive in terms of the differentiation of the layers. In the fourth soil assessment, even in the intermediate layers, the percentage of success was higher than 89.1% in the calibration step and 68.4% in the validation step (Figure 4).

The results from the discriminant analysis indicate that the differences among the sampled soil layers had a larger influence on the PLSR models than the differences between the applied treatments.

3.3. Spectral Modeling: iPLS

An example of the selection of iPLS variables from the reflectance curves (400–2400 nm) for CEC prediction in the second soil assessment can be seen in Figure 5. A validation PLSR model was applied to each of the intervals of the reflectance curves, and after variable selection, a decrease was observed in the RMSE from 0.608 (red line, using all variables) to 0.34 cmol_c dm⁻³ (blue line, using only the selected variables).

The wavelengths selected by the iPLS algorithm for each soil characteristic are shown in Figure 6. Note that for each chemical characteristic, different wavelengths were selected in the four assessments due to the variability in these characteristics (Table 2) over the assessment period. For characteristics with lower variability, such as CEC, the selected wavelengths were more similar among the soil assessments. Furthermore, for the Ca²⁺ and Mg²⁺ content, the selected wavelengths of close to 1.690 and 2.300 nm are related to the applied dolomitic limestone (CaMg (CO₃)₂) [65].

3.4. Spectral Modeling: PLSR

The statistical parameters of the PLSR models developed in the first soil assessment are shown in

Table 3. Statistical parameters of PLSR models developed from the first soil assessment (with iPLS selection).

PLSR		Ca2+	Mg2+	K+	P	S-SO42−	pH	Al3+	H+Al	SB	CECe	CEC	m	BS
		cmolc dm−3			mg dm−3			cmolc dm−3					%
	Variables	109	138	153	67	98	232	25	126	96	132	71	136	63
	Factors	5	7	3	3	1	4	3	3	5	5	3	5	4
Calibration	R2	0.76	0.58	0.76	0.74	0.00	0.59	0.10	0.18	0.81	0.86	0.84	0.26	0.70
	RMSE	0.55	0.35	0.05	15.16	28.24	0.30	0.16	0.82	0.70	0.63	0.54	9.19	8.40
	BIAS	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	RPD	2.04	1.54	2.04	1.96	1.00	1.56	1.05	1.10	2.29	2.71	2.50	1.16	1.83
	RPIQ	2.58	1.68	1.92	2.76	0.35	2.09	1.36	1.30	2.97	3.19	4.04	1.44	2.52
Validation	R2	0.74	0.51	0.75	0.73	0.00	0.56	0.07	0.14	0.79	0.84	0.83	0.19	0.68
	RMSE	0.57	0.38	0.06	15.62	28.72	0.32	0.17	0.85	0.73	0.68	0.56	9.67	8.71
	BIAS	0.06	0.00	0.00	−0.01	−0.07	0.00	0.00	0.00	0.00	0.00	0.00	0.01	−0.02
	RPD	1.96	1.43	2.00	1.92	1.00	1.51	1.04	1.08	2.18	2.51	2.43	1.11	1.77
	RPIQ	2.49	1.54	1.60	2.68	0.35	1.96	1.28	1.26	2.85	2.96	3.89	1.37	2.43
Prediction	R2	0.75	0.00	0.81	0.82	0.00	0.53	0.01	0.27	0.72	0.81	0.73	0.30	0.63
	RMSE	0.71	0.74	0.06	14.01	28.52	0.30	0.18	0.81	1.04	0.61	0.66	8.97	10.37
	BIAS	0.36	0.58	0.01	2.02	−3.45	−0.02	0.06	0.02	0.23	−0.03	0.45	−1.11	5.84
	RPD	2.00	1.00	2.29	2.36	1.00	1.46	1.01	1.17	1.89	2.80	1.92	1.20	1.64
	RPIQ	2.00	0.79	1.60	2.98	0.35	2.09	1.21	1.32	2.00	3.30	3.30	1.47	2.04

Ca²⁺, Mg²⁺, and Al³⁺ contents were extracted using KCl 1.0 mol L⁻¹; K⁺ and P contents were extracted using Mehlich-1; pH CaCl₂ (0.01 mol L⁻¹) at a soil:solution ratio of 1:2.5; H+Al was extracted using the SMP method; CECe is the effective cation exchange capacity; CEC is the cation exchange capacity at pH 7; m indicates the aluminum saturation; BS is the base saturation; R² is the coefficient of determination; RMSE is the root mean square error; BIAS is the systematic error; RPD represents the relative percentage deviation; RPIQ is the ratio of performance to interquartile range.

. The characteristics of Ca²⁺ and K⁺ content, P, SB, CECe, CEC, and BS had values for R² ranging from 0.70 to 0.86 in the calibration and validation steps, and were therefore classified as good or very good prediction models. In the prediction step, the values of R² ranged from 0.63 to 0.82, giving models with moderate to very good quantitative prediction capabilities. Our values were similar to those found by Rodrigues et al. [28] for Ca²⁺, K⁺, and P (with R² above 0.76 in the calibration and validation steps), and were higher than those found by Veum et al. [66] for Ca²⁺, K⁺, and P.

Of the PLSR models developed in the first soil assessment (Table 3), those for K⁺, P, and CECe stood out, with values of R² of above 0.81. For P, the prediction R² was 0.82, while for K⁺, the value was 0.81, both of which were higher than those reported by Zornoza et al. [22]. With R² values of above 0.80, these models provide a very good estimate of the characteristics. For Mg²⁺, it was not possible to perform the prediction step, possibly due to the low value of R² in both the calibration and the validation steps. For S-SO₄²⁻, Al³⁺, and m%, it was not possible to develop good models, as they could not explain the variability of the data.

For the RMSE value in the validation step (Table 3), there was a slight variation compared to the calibration step, indicating models with good predictive capacity. These low RMSE values agree with the coefficients of determination, which were classified as good and very good. The same results were found in the other evaluations (

Table 4. Statistical parameters of PLSR models developed from the second soil assessment (with iPLS selection).

PLSR		Ca2+	Mg2+	K+	P	S-SO42−	pH	Al3+	H+Al	SB	CECe	CEC	m	BS
		cmolc dm−3			mg dm−3			cmolc dm−3					%
	Variables	136	181	172	168	36	220	237	175	124	130	93	82	170
	Factors	9	5	6	5	1	9	4	9	6	6	3	1	7
Calibration	R2	0.84	0.65	0.74	0.76	0.07	0.62	0.11	0.48	0.86	0.85	0.90	0.23	0.76
	RMSE	0.59	0.42	0.10	8.14	18.76	0.23	0.20	0.69	0.83	0.77	0.46	9.77	8.57
	BIAS	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	RPD	2.50	1.69	1.96	2.04	1.04	1.62	1.06	1.39	2.67	2.58	3.16	1.14	2.04
	RPIQ	3.09	1.67	1.99	2.74	0.69	2.22	1.29	1.88	3.29	3.55	5.42	1.53	2.73
Validation	R2	0.80	0.62	0.73	0.74	0.06	0.54	0.05	0.38	0.84	0.84	0.90	0.16	0.72
	RMSE	0.65	0.43	0.11	8.59	18.96	0.26	0.21	0.76	0.87	0.82	0.47	10.23	9.15
	BIAS	0.00	0.00	0.00	0.00	0.01	0.00	0.00	0.00	0.00	0.00	0.00	−0.03	0.02
	RPD	2.24	1.62	1.92	1.96	1.03	1.47	1.03	1.27	2.50	2.50	3.16	1.09	1.89
	RPIQ	2.80	1.63	1.81	2.59	0.68	1.96	1.23	1.71	3.14	3.34	5.30	1.46	2.55
Prediction	R2	0.81	0.47	0.70	0.67	0.09	0.20	0.00	0.00	0.85	0.86	0.91	0.03	0.72
	RMSE	0.68	0.53	0.10	8.91	23.86	0.36	0.23	1.26	0.90	0.81	0.51	11.19	9.45
	BIAS	0.12	−0.03	−0.03	1.51	−1.49	−0.12	−0.07	−0.89	−0.17	0.22	0.28	3.07	−0.01
	RPD	2.29	1.37	1.83	1.74	1.05	1.12	1.00	1.00	2.58	2.67	3.33	1.02	1.89
	RPIQ	2.68	1.33	1.99	2.50	0.54	1.42	1.12	1.03	3.03	3.38	4.88	1.33	2.47

Ca²⁺, Mg²⁺, and Al³⁺ contents were extracted using KCl 1.0 mol L⁻¹; K⁺ and P contents were extracted using Mehlich-1; pH CaCl₂ (0.01 mol L⁻¹) at a soil:solution ratio of 1:2.5; H+Al was extracted using the SMP method; CECe is the effective cation exchange capacity; CEC is the cation exchange capacity at pH 7; m represents the aluminum saturation; BS is the base saturation; R² is the coefficient of determination; RMSE is the root mean square error; BIAS is the systematic error; RPD is the relative percentage deviation; and RPIQ is the ratio of performance to interquartile range.

,

Table 5. Statistical parameters of PLSR models developed from the third soil assessment (with iPLS selection).

PLSR		Ca2+	Mg2+	K+	P	S-SO42−	pH	Al3+	H+Al	SB	CECe	CEC	m	BS
		cmolc dm−3			mg dm−3			cmolc dm−3					%
	Variables	81	61	88	24	18	32	41	2	44	82	58	65	35
	Factors	8	10	5	4	7	10	8	2	7	8	4	8	4
Calibration	R2	0.82	0.73	0.78	0.77	0.18	0.74	0.36	0.05	0.86	0.86	0.87	0.39	0.67
	RMSE	0.65	0.40	0.07	8.01	18.81	0.23	0.21	1.13	0.89	0.84	0.73	10.09	9.64
	BIAS	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	RPD	2.36	1.92	2.13	2.09	1.10	1.96	1.25	1.03	2.67	2.67	2.77	1.28	1.74
	RPIQ	3.14	2.02	2.85	2.94	0.48	2.89	1.61	1.03	3.67	3.63	4.24	1.60	2.71
Validation	R2	0.81	0.69	0.75	0.77	0.09	0.70	0.29	0.03	0.84	0.82	0.86	0.32	0.65
	RMSE	0.68	0.43	0.07	8.35	19.90	0.25	0.22	1.15	0.93	0.95	0.75	10.74	9.95
	BIAS	0.00	0.00	0.00	0.05	0.00	0.00	0.00	0.00	0.01	0.01	0.00	−0.03	−0.02
	RPD	2.29	1.80	2.00	2.09	1.05	1.83	1.19	1.02	2.50	2.36	2.67	1.21	1.69
	RPIQ	3.00	1.87	2.85	2.82	0.46	2.66	1.53	1.01	3.51	3.21	4.13	1.50	2.62
Prediction	R2	0.71	0.00	0.76	0.28	0.00	0.00	0.00	0.00	0.77	0.87	0.80	0.18	0.66
	RMSE	0.98	0.85	0.07	11.11	108.20	4.99	0.77	0.96	1.08	0.76	0.85	11.80	10.06
	BIAS	0.53	−0.66	0.01	−7.50	104.70	4.96	−0.73	0.12	0.02	−0.03	0.18	−0.09	1.16
	RPD	1.86	1.00	2.04	1.18	1.00	1.00	1.00	1.00	2.09	2.77	2.24	1.10	1.71
	RPIQ	2.08	0.95	2.85	2.12	0.08	0.13	0.44	1.21	3.02	4.01	3.65	1.36	2.59

Ca²⁺, Mg²⁺, and Al³⁺ contents were extracted using KCl 1.0 mol L⁻¹; K⁺ and P contents were extracted using Mehlich-1; pH CaCl₂ (0.01 mol L⁻¹) at a soil:solution ratio of 1:2.5; H+Al was extracted using the SMP method; CECe is the effective cation exchange capacity; CEC is the cation exchange capacity at pH 7; m represents the aluminum saturation; BS is the base saturation; R² is the coefficient of determination; RMSE is the root mean square error; BIAS is the systematic error; RPD is the relative percentage deviation; and RPIQ is the ratio of performance to interquartile range.

and

Table 6. Statistical parameters of the PLSR models developed from the fourth soil assessment (with iPLS selection).

PLSR		Ca2+	Mg2+	K+	P	S-SO42−	pH	Al3+	H+Al	SB	CECe	CEC	m	BS
		cmolc dm−3			mg dm−3			cmolc dm−3					%
	Variables	134	81	95	64	81	137	126	101	90	117	231	146	47
	Factors	10	5	3	4	3	10	9	10	7	7	5	6	10
Calibration	R2	0.92	0.63	0.84	0.84	0.16	0.75	0.49	0.67	0.88	0.89	0.91	0.47	0.87
	RMSE	0.31	0.42	0.07	6.51	5.11	0.18	0.16	0.37	0.62	0.56	0.50	9.27	5.47
	BIAS	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00	0.00
	RPD	3.54	1.64	2.50	2.50	1.09	2.00	1.40	1.74	2.89	3.02	3.33	1.37	2.77
	RPIQ	4.05	2.14	2.85	3.47	0.90	2.50	2.11	2.27	4.28	4.23	4.96	1.75	4.34
Validation	R2	0.89	0.59	0.83	0.82	0.12	0.68	0.38	0.59	0.87	0.87	0.89	0.37	0.83
	RMSE	0.36	0.44	0.07	6.82	5.28	0.21	0.17	0.42	0.65	0.61	0.56	9.01	6.38
	BIAS	0.00	0.00	0.00	0.00	−0.01	0.00	0.00	0.00	0.01	−0.01	−0.05	−0.18	0.04
	RPD	3.02	1.56	2.43	2.36	1.07	1.77	1.27	1.56	2.77	2.77	3.02	1.26	2.43
	RPIQ	3.49	2.05	2.85	3.31	0.88	2.14	1.99	2.00	4.08	3.88	4.43	1.80	3.72
Prediction	R2	0.84	0.01	0.75	0.58	0.19	0.00	0.00	0.00	0.86	0.89	0.90	0.00	0.90
	RMSE	0.54	0.45	0.09	7.90	6.18	3.15	0.36	1.66	0.64	0.84	0.48	15.42	4.79
	BIAS	0.84	0.30	0.01	3.96	1.06	−3.14	−0.08	1.48	0.08	−0.49	0.17	−12.64	−0.12
	RPD	2.50	1.01	2.00	1.54	1.11	1.00	1.00	1.00	2.67	3.02	3.16	1.00	3.16
	RPIQ	2.33	2.0	2.22	2.86	0.75	0.14	0.94	0.51	4.15	2.82	5.17	1.05	4.96

Ca²⁺, Mg²⁺, and Al³⁺ contents were extracted using KCl 1.0 mol L⁻¹; K⁺ and P contents were extracted using Mehlich-1; pH CaCl₂ (0.01 mol L⁻¹) at a soil:solution ratio of 1:2.5; H+Al was extracted using the SMP method; CECe is the effective cation exchange capacity; CEC is the cation exchange capacity at pH 7; m represents the aluminum saturation; BS is the base saturation; R² is the coefficient of determination; RMSE is the root mean square error; BIAS is the systematic error; RPD is the relative percentage deviation; and RPIQ is the ratio of performance to interquartile range.

).

Based on the RPD results (Table 3), the PLSR models in the calibration phase were classified as highly reliable for Ca²⁺, K⁺, SB, CECe and CEC, reliable for Mg²⁺, P, pH, and BS, and unreliable for S-SO₄²⁻, Al³⁺, H+Al, and m%. Our RPD values were higher than those reported by Veum et al. [66] for Ca²⁺, Mg²⁺, K⁺, P, and pH.

Based on the RPIQ values from the first evaluation (Table 3), the PLSR models were classified as poor for Ca²⁺, pH, and BS (in which only high and low values are distinguishable), regular for P, SB, and CECe (in which quantitative predictions are possible), and good for CEC (representing a good quantitative model). These results show that the RPIQ parameter proved to be more rigorous than the other parameters for the classification of the developed models; even so, it was still possible to obtain regular and good models. In particular, the PLSR model for CEC was classified as good based on the RPIQ. Both the R² and the RPD demonstrated that this was a good and highly reliable model.

In general, an improvement was seen in most PLSR models with iPLS selection (Table 3) compared to PLSR models without iPLS selection (

Table A1. Statistical parameters of the PLSR models developed from the four soil assessments (without iPLS selection).

PLSR		Ca2+	Mg2+	K+	P	S-SO42−	pH	Al3+	H+Al	SB	CECe	CEC	m	BS
		cmolc dm−3			mg dm−3			cmolc dm−3					%
First Assessment
Factors		5	5	5	2	1	3	6	7	5	5	5	6	5
Cal.	R2	0.75	0.49	0.76	0.62	0.01	0.42	0.23	0.31	0.77	0.79	0.76	0.31	0.67
	RMSE	0.59	0.37	0.06	20.15	27.69	0.36	0.15	0.76	0.81	0.74	0.67	8.54	8.98
Val.	R2	0.70	0.43	0.74	0.61	0.00	0.40	0.14	0.18	0.72	0.75	0.74	0.21	0.62
	RMSE	0.66	0.40	0.06	20.54	28.04	0.37	0.16	0.82	0.89	0.80	0.71	9.17	9.68
Pred.	R2	0.73	0.39	0.69	0.60	0.02	0.44	0.14	0.30	0.75	0.77	0.72	0.22	0.65
	RMSE	0.70	0.39	0.07	21.39	29.61	0.37	0.16	0.80	0.89	0.81	0.71	10.43	9.88
Second Assessment
Factors		8	5	5	9	1	9	2	9	5	5	6	2	5
Cal.	R2	0.78	0.51	0.61	0.72	0.04	0.53	0.05	0.40	0.75	0.78	0.90	0.14	0.62
	RMSE	0.74	0.50	0.13	9.76	19.10	0.27	0.22	0.80	1.19	1.05	0.56	10.91	11.05
Val.	R2	0.74	0.46	0.59	0.65	0.03	0.43	0.03	0.26	0.71	0.75	0.88	0.13	0.59
	RMSE	0.82	0.53	0.14	10.96	19.32	0.30	0.22	0.89	1.27	1.11	0.61	11.01	11.58
Pred.	R2	0.69	0.32	0.57	0.61	0.05	0.32	0.00	0.16	0.65	0.70	0.86	0.00	0.47
	RMSE	0.78	0.57	0.12	10.04	24.17	0.31	0.20	0.78	1.24	1.10	0.62	9.46	11.91
Third Assessment
Factors		6	10	6	3	1	10	7	8	6	7	5	9	7
Cal.	R2	0.76	0.68	0.76	0.72	0.00	0.53	0.27	0.27	0.79	0.83	0.80	0.39	0.67
	RMSE	0.77	0.43	0.07	9.14	20.67	0.29	0.24	0.92	1.08	0.91	1.01	10.27	9.64
Val.	R2	0.74	0.58	0.73	0.70	0.00	0.41	0.16	0.12	0.76	0.80	0.78	0.26	0.61
	RMSE	0.82	0.49	0.08	9.52	20.87	0.33	0.25	1.01	1.17	1.01	1.07	11.34	10.51
Pred.	R2	0.71	0.48	0.64	0.72	0.00	0.36	0.13	0.20	0.80	0.84	0.81	0.10	0.64
	RMSE	0.92	0.47	0.09	8.81	11.19	0.34	0.24	1.08	1.04	0.88	0.86	11.89	10.34
Fourth Assessment
Factors		6	5	3	4	2	10	10	9	6	5	4	10	6
Cal.	R2	0.86	0.60	0.80	0.82	0.17	0.66	0.53	0.54	0.87	0.88	0.88	0.60	0.79
	RMSE	0.43	0.41	0.08	6.34	4.86	0.21	0.16	0.45	0.66	0.59	0.57	8.15	6.97
Val.	R2	0.84	0.55	0.79	0.81	0.15	0.56	0.38	0.40	0.85	0.86	0.87	0.47	0.76
	RMSE	0.47	0.43	0.08	6.67	4.94	0.25	0.18	0.51	0.70	0.64	0.59	9.54	7.48
Pred.	R2	0.82	0.57	0.65	0.73	0.09	0.43	0.10	0.30	0.84	0.86	0.81	0.27	0.79
	RMSE	0.46	0.43	0.10	8.08	7.06	0.28	0.18	0.63	0.70	0.63	0.72	8.78	7.09

Ca²⁺, Mg²⁺, and Al³⁺ contents were extracted using KCl 1.0 mol L⁻¹; K⁺ and P contents were extracted using Mehlich-1; pH CaCl₂ (0.01 mol L⁻¹) at a soil:solution rate of 1:2.5; H+Al was extracted using the SMP method; CECe: effective cation exchange capacity; CEC: cation exchange capacity at pH 7; m: aluminum saturation; BS: base saturation; R²: coefficient of determination; RMSE: root mean square error.

), although no improvement was observed for S-SO4²⁻, pH, Al⁺³, H+Al, and m%. The R² values in the calibration step were higher by about 16% for Mg²⁺ and P, 9% for CECe and CEC, and 4.5% for SB and BS (Table A1). For Ca²⁺ and K⁺, no increase in R² was observed. There was an average reduction of 16% in the calibration RMSE with iPLS for Ca²⁺, Mg²⁺, K⁺, P, SB, CECe, CEC, and BS compared to the PLSR models without iPLS selection; in particular, the RMSE values for K+, P and CEC showed reductions of 20%, 33% and 24%, respectively. In the validation phase, there was an average increase of 10% in the R² values for Ca²⁺, Mg²⁺, K⁺, P, SB, CECe, CEC, and BS, where the values for Mg²⁺, P, CECe, and CEC showed increases of 16%, 16%, 11% and 11%, respectively. In the same way as for the calibration step, there was an average reduction of 16% in the RMSE for the validation step with iPLS, and in particular, the RMSE values for P, SB, CECe, and CEC showed reductions of 31%, 22%, 18%, and 27%, respectively. K+ and P stood out in the prediction phase, as they showed increases in R² of 15% and 27% and reductions in the RMSE of 17% and 53%, respectively. The prediction RMSE for CEC was also noteworthy, with a reduction of 33%.

In the second soil assessment (Table 4), the PLSR models of Ca²⁺, K⁺, P, SB, CECe, CEC, and BS achieved the highest values of R², ranging from 0.73 to 0.90 in the calibration and validation steps; in the prediction step, the values of R² ranged from 0.67 to 0.91. This range of R² values indicates that the models have good to excellent quantitative prediction capabilities. These values were higher than those reported by Pinheiro et al. [26] and Veum et al. [66]. The results from the second soil assessment were also superior to those of the first assessment (Table 3).

In addition, the PLSR models developed for the prediction of SB, CECe, and CEC achieved R² values of above 0.85; the value for CEC reached 0.91 (Table 4), and the model had lower numbers of variables (93) and factors (three) compared to those of SB and CECe. For S-SO₄²⁻, Al³⁺, H+Al, and m%, it was not possible to develop good models in the second soil assessment (Table 4), as they could not explain the data variability.

Based on the RPD values (Table 4), the PLSR models of the calibration phase were classified as highly reliable for Ca²⁺, P, SB, CECe, CEC and BS, reliable for Mg²⁺, K⁺ and pH, and unreliable for S-SO₄²⁻, Al³⁺, H+Al, and m%. Note that there was an improvement in the RPD for BS compared to the first evaluation, which was accompanied by an increase in R² values in all phases.

Based on the RPIQ values for the second evaluation (Table 4), the PLSR models were classified as poor for pH (meaning that only high and low values were distinguishable), regular for Ca²⁺, P, SB, and BS (meaning that quantitative predictions were possible), good for CECe, and excellent for CEC. These RPIQ values show an improvement in the predictive capacity of the models compared to the first evaluation, mainly for Ca²⁺, CECe, CEC, and BS. The PLSR model for CEC, in particular, was classified as excellent based on the RPIQ, and was shown to be a very good and highly reliable model based on both the R² and the RPD.

The better accuracy in the second evaluation can be explained by the longer reaction time for the limestone and phosphogypsum, which allowed for higher values of SB, CECe, CEC, and BS (Table 2). As a result of this greater amplitude, more robust models were obtained with a better fit.

In general, there was an improvement in most PLSR models with iPLS selection (Table 4) compared to models without iPLS selection (Table A1). However, again, no improvements were observed for S-SO₄²⁻, pH, Al⁺³, or H+Al, and these results were repeated in the other evaluations. After iPLS selection, the R² values in the calibration step were increased by 22% for Mg²⁺, 18% for K⁺, and 13% for SB and BS. Regarding the calibration RMSE, there was an average reduction of 28% for Ca²⁺, Mg²⁺, K⁺, P, SB, CECe, CEC, and BS compared to PLSR models without iPLS selection, where the results for K⁺, SB, CEC, and BS showed reductions of 30%, 43%, 36% and 29%, respectively. In the validation phase, there was an average increase of 14% in the R² values for Ca²⁺, Mg²⁺, K⁺, P, SB, CECe, CEC, and BS; the R² values for Mg²⁺, K⁺, SB, and BS showed increases of 26%, 19%, 15% and 18%, respectively. In the validation stage, there was an average reduction of 30% of the RMSE with iPLS, with the RMSE for SB, CECe, and CEC showing reductions of 46%, 35%, and 30%, respectively. In the prediction phase, the results for SB, CECe, and BS showed increases of 24%, 19%, and 35% in the values of R², and reductions in the RMSE of 38%, 36%, and 26%, respectively. Another noteworthy result was the prediction RMSE for CEC, with a reduction of 22%.

PLSR models with similar adjustments to the second soil assessment were also found for the third (Table 5) and fourth soil assessments (Table 6), two and three years after the treatment application (2016), respectively. In the fourth soil assessment, the PLSR models of Ca²⁺, K⁺, SB, CECe, CEC, and BS achieved the highest values of R², ranging from 0.83 to 0.92 in the calibration and validation steps and from 0.70 to 0.90 for the prediction step. This range of R² values indicates that these models have a good or very good capacity for quantitative prediction. These values were higher than those reported by Pinheiro et al. [26] and Terra et al. [25], but lower than those found by Rodrigues et al. [28] for Ca²⁺ and K⁺, for all three steps (calibration, validation, and prediction). However, it was not possible to develop good models for S-SO₄²⁻, pH, Al³⁺, H+Al, and m% for the calibration, validation, and prediction steps of the third and fourth soil assessments (Table 5 and Table 6). The prediction models for CEC and SB for the fourth soil assessment (Table 6) were notable, with R² values of 0.89 and 0.90, respectively, indicating a very good or excellent capacity for quantitative prediction.

For the RPD values, the results for the third assessment were similar to those of the second assessment (Table 5). The PLSR models for the calibration phase were classified as highly reliable for Ca²⁺, K⁺, P, SB, CECe, and CEC, reliable for Mg²⁺, pH and BS, and unreliable for S-SO₄²⁻, Al³⁺, H+Al, and m%.

Based on the RPIQ values of the third assessment (Table 5), the PLSR models were classified as poor for Mg²⁺ (meaning that only high and low values were distinguishable), regular for Ca²⁺, K⁺, P, pH, and BS (meaning that quantitative predictions were possible), good for SB and CECe (meaning that these were good quantitative models), and excellent for CEC. There was an improvement in the predictive capacity of the SB model compared with the second evaluation, and this was classified as a good quantitative model. The classification of the PLSR model for CEC remained excellent, based on the RPIQ.

A comparison between PLSR models with iPLS selection (Table 5) and PLSR models without iPLS selection (Table A1) indicates that the R² values for the calibration stage increased by 7% for Ca²⁺ and Mg²⁺, and 8% for SB and CEC. In terms of the calibration RMSE, there was an average reduction of 14% for Ca²⁺, Mg²⁺, K⁺, P, SB, CECe, CEC, and BS compared to the PLSR models without iPLS selection; in particular, the SB and CEC models showed reductions of 21% and 38%, respectively. In the validation phase, there was an average increase of 8% in the R² values for Ca²⁺, Mg²⁺, K⁺, P, SB, CECe, CEC, and BS. The R² values for Mg²⁺ and SB were notable, with increases of 16% and 10%, respectively. Regarding the validation RMSE, there was an average reduction of 18% in the RMSE with iPLS, and in particular, the RMSE values for Ca²⁺, SB, and CEC showed reductions of 21%, 26%, and 43%, respectively. In the prediction phase, the results for K⁺ were notable, and showed an increase of 16% in R² and a reduction in RMSE of 29%. However, in the prediction phase, the PLSR models with iPLS developed in the third evaluation were not improved for most features, compared to the PLSR models without iPLS.

The RPD values found in the fourth evaluation (Table 6) indicate that the PLSR models of the calibration phase can be classified as highly reliable for Ca²⁺, K⁺, P, SB, CECe, CEC, and BS, reliable for Mg²⁺ and pH, and unreliable for S-SO₄²⁻, Al³⁺, H+Al, and m. The classification for the BS model was elevated to highly reliable, demonstrating the evolution of the BS parameters throughout the evaluations.

Based on the RPIQ values for the fourth evaluation (Table 6), the PLSR models were classified as poor for Mg²⁺, pH, Al³⁺, H+Al (meaning that only high and low values are distinguishable), regular for K⁺ (meaning that quantitative predictions are possible), good for Ca²⁺ and P, and excellent for SB, CECe, CEC, and BS. The predictive capacity of the SB, CECe, and BS models improved compared to the other assessments, and these were classified as excellent quantitative models.

The PLSR model for CEC remained excellent based on the RPIQ values over the last three evaluations, a result that can be explained by the association of CEC with mineralogy and soil granulometry [2]. Numerous studies have been performed using Vis-NIR-SWIR reflectance spectroscopy to predict the physical characteristics of soil, such as granulometry [50,67]. The relationship of Vis-NIR-SWIR spectroscopy with the characteristics and physics of soil is already well established; however, works that include exchangeable ions, SB, CECe, CEC, and BS are still scarce, and mainly involve variable selection.

Since both the CEC and BS are necessary to calculate liming recommendations using the BS method [5], the results for R², RPD, and RPIQ show that it is possible to calculate liming recommendations based on the values estimated by the CEC and BS PLSR models.

Once again, an improvement was noted in most PLSR models with iPLS selection (Table 6) compared to models without iPLS selection (Table A1). After iPLS selection, the R² values at the calibration stage were on average 4% higher, and in particular, the results for Ca²⁺ and BS showed increases of 7% and 9%, respectively. There was an average reduction of 18% in the calibration RMSE for models with iPLS for Ca²⁺, K⁺, SB, CECe, CEC, and BS compared to the PLSR models without iPLS selection; notable results were seen for Ca²⁺ and BS, which showed reductions of 39% and 27%, respectively. In the validation phase, there was an average increase of 4% in the R² values for Ca²⁺, Mg²⁺, K⁺, P, SB, CECe, CEC, and BS, and the values of R² for Mg²⁺ and BS increased by 7% and 8%, respectively. Again, in the validation phase, there was an average reduction of 13% in the RMSE for models with iPLS, where the values for Ca²⁺, K⁺, and BS showed reductions of 31%, 14%, and 17%, respectively. The results for K+ and BS were notable in the prediction phase, as these showed increases of 13% and 12% in R² and reductions of 11% and 48% in RMSE, respectively. The prediction results for the RMSE for CEC were also significant, with a reduction of 50%. For Ca²⁺ and Mg^2+, there were increases in RMSE of 15% and 4%, respectively, in the models with iPLS. These results indicate a greater reduction in the RMSE in the models with iPLS, particularly for the CEC and BS, which are characteristics used for liming recommendation in the BS method [5].

In most cases, the reduction in the RMSE values obtained by the PLSR models with iPLS (in the four evaluations) improved the predictive capacity of these models, since the RMSE measures the differences between the estimated values and the real values, and hence quantifies the accuracy of the model. Reductions in the RMSE values therefore gave rise to better accuracy for these models.

Furthermore, for all the characteristics studied for the four soil samples (Table 3, Table 4, Table 5 and Table 6), the systematic error values (BIAS) in the calibration and validation phases were close to zero; these results demonstrate the very low bias of the model in regard to the estimated characteristics [68]. Willmott [69] states that PLSR models are considered to be good when the BIAS approaches zero, similarly to the results in the prediction phase for characteristics with R² above 0.6 (Table 3, Table 4, Table 5 and Table 6).

Thus, there was a gradual evolution of the statistical parameters from the first to the fourth assessment; this may have been because as the treatments (limestone and phosphogypsum) reacted with the soil over the years, greater amplitudes in the chemical characteristics of the soil were generated, allowing us to find better models.

3.5. Spectral-Based Method for Liming Recommendation

Figure 7 shows the results of a correlation analysis between the estimated liming recommendation using the characteristics obtained by the PLSR models, and the calculated liming recommendation using the characteristics obtained by conventional laboratory techniques. The first soil assessment showed the lowest correlation (r = 0.86), and the correlation values increased up to the fourth soil assessment, for which the correlation coefficient was 0.95. These results therefore show that over the four soil assessments (5, 12, 24, and 36 months after the application of treatments), the liming recommendations based on the chemical characteristics obtained by the PLSR models became progressively more correlated with the liming recommendations calculated based on the chemical characteristics obtained by conventional laboratory techniques. This was probably due to the limestone reaction time. Several studies of Brazilian soils have demonstrated that three years is sufficient for all of the applied limestone to undergo reaction [6,70,71]. Thus, for our very clayey Oxisol, liming recommendations based on multivariate statistical modeling of soil reflectance should preferably be made three years after the last liming, as in this case, the prediction models for CEC and BS had values of R² for prediction of 0.90 (Table 6).

4. Conclusions

Our models were shown to discriminate between the layers sampled from the soil profile with high accuracy, although models discriminating among treatments could not be developed.

The use of iPLS selection gave rise to an increase in R² and a reduction in the RMSE values for most PLSR models developed for all four soil assessments, which improved the accuracy of these models.

We were able to estimate the soil chemical characteristics by selecting wavelengths using the iPLS method and PLSR modeling, and, in particular, for Ca²⁺, SB, CTCe, CTC, and BS. The statistical parameters of these models gradually improved until the fourth assessment. The values of R² ranged from 0.83 to 0.92 in the calibration and validation steps, and from 0.84 to 0.90 in the prediction step (in the fourth assessment).

The PLSR models of CEC and BS made it possible to calculate a liming recommendation based on the chemical characteristics estimated by the PLSR models. This estimate was positively correlated (r > 0.86) with the liming recommendation calculated based on the chemical characteristics obtained by conventional laboratory techniques. The fourth soil assessment yielded the best correlation coefficient (r = 0.95).

Our results demonstrate that the use of Vis-NIR-SWIR spectroscopy in conjunction with the selection of iPLS variables technique has strong potential in terms of predicting soil chemical characteristics, and can complement conventional analyses of these characteristics. This approach has several advantages, as it has a low cost, does not use toxic reagents, and can be used when many soil samples are required, for example in precision agriculture. In addition, the possibility of using a smaller number of wavelengths to estimate soil chemical characteristics may enable the use of cheaper sensors and automated recommendations for liming tropical soils. However, more work on wavelength selection for PLSR modeling needs to be done.