Evaluation of Leymus chinensis quality using near-infrared reflectance spectroscopy with three different statistical analyses

Sample collection and pretreatment

The sampling sites were from the west to east across the grassland in Heilongjiang Province of Northeast China. A total of 203 samples of sheepgrass hay were randomly collected from sheepgrass fields of the hay factories in 2013, with latitudes ranging from 44.475°N to 51.728°N, longitudes from 123.209°E to 132.944°E. The sites were chosen to be representative of sheepgrass production fields and contain a range of soils and climate are described in Appendix S1. These locations produce approximately 117.19 million tons of sheepgrass hay per year. All samples of sheepgrass hay were sampled at blooming stages, identified and collected before being clipping and packaging. The collected samples (not involve endangered or protected species), a total sample size of 203 from 21 sampling sites distributing over 8 regions (not privately owned or protected in any way) in Heilongjiang Province (see Fig. 1 and Appendix S1), were then forwarded to the lab and stored at 4 °C for further analyses.

To be representative, each sample consisted of one quarter square meter clipped at 4 cm, transported to the lab, oven dried (65 °C, 48 h), ground (1 mm sieve), and mixed. Meanwhile, the mixed samples were divided by the quartile method; half for duplicate chemical analyses and half for the near infrared reflectance spectra.

In fact, all samples (203) were used to evaluate and develop NIRs models in this study. To ensure the adaptability of the calibration models, some samples (51) were used for prediction set and the rest of samples (152) for calibration set. All tasks including spectral measure and chemical analyses were finished on the same day or the next day.

Chemical properties analyses

The first experiment was designed to develop a database to evaluate relationships between the quality indices of sheepgrass hay and NIRs measurements. In order to accomplish this goal, three characteristics were measured at blooming stages of growth. These included CP, ADF, and NDF, which were usually regarded as the principal forage quality paremeters, serving as the primary nutrition source in the diet of dairy cattle (Yan et al., 2011; Mertens, 1992). All chemical analyses were conducted in duplicate using analytical grade chemicals are described in Appendix S2. The CP was quantified using the Kjeldahl procedure for nitrogen (N) determination (%CP =%N × 6.25) by the methods of Association of Official Analytical Chemists-AOAC (AOAC, 1990). The ADF and NDF were conducted by using the method described by Van Soest, Robertson & Lewis (1991).

Spectroscopic measurement

By using an near infrared spectrometer (DA7200, Perten Corporation, Hägersten, Sweden) to collect for the spectra of samples, the NIR spectra were recorded at 5 nm intervals from 950 to 1,650 nm in Appendix S2. Approximately 50 g of dried sample was scanned in a 7.5 cm diameter sample cell with a quartz window in the room temperature maintained at 25 °C. With the working parameters of DA7200 software, samples were scanned from 950 to 1,650 nm in reflectance mode (R mode). Data form was converted to the absorbance (the logarithms of the reflectance reciprocal). In order to eliminate an error caused by loading the sample, each sample were repeated to scan three times and each time with scanning background light three times. The average spectral scanning were as for the final spectra of samples deposited in the computer for calculating using the Grams32 software (Perten Corporation, Hägersten, Sweden).

Statistical analysis

To obtain linear correlations of the NIRs with the chemical values, spectra were expressed as the absorbance A (A=log1/R). The NIRs were conducted by Grams32 software and Microsoft Office Excel 2003. In performing measurements with the NIR System DA7200 spectrometer, a number of data were generated for each sample and it is obvious that some data reduction method is needed to facilitate data interpretation. Thus, there is a need for data processing methods that transform the measured spectral data into the sample properties of interest. For producing such models that determine the equations describing the relationship between spectral data and chemical values, three different methods (MLR, PCR and PLS) were used. MLR uses a multiterm linear polynomial to describe this relationship using only several spectral data measured as “characteristic” wavelengths (by the way of correlation spectrum). With PCR, all spectral data are used, reducing the dimensionality of the data set by looking for orthogonal directions in spectral data space along which the variance of the data set is maximised. These directions are called principal components (PC). Thus, the first principal component (PC 1) is determined as the direction in spectral data space that corresponds to the largest variation in the data set. The second principal component (PC 2) is calculated as the direction perpendicular to PC 1 along which the remaining variation in the data set is the largest. This procedure is repeated until no variation is left in the data set. PLS uses not only the spectral information in the data set but also incorporates information about sample properties, e.g., concentrations to determine the most useful orthogonal directions in spectral space.

An ideal equation would have a coefficient of determination (R²) of 1.0 and root mean squared error (RMSE). Approximately 75% of the data were used in creating the equation, while the remaining 25% were randomly removed for validation. Validation procedures for the created equations were the same as described above, using the remaining 25% of the data. In this study, the calibration set (152 samples) and the validation set (51 samples) were randomly divided by the method applied by several publications (Miriam et al., 2011; Peltre et al., 2011; Hayes, 2012). To restrain invariable background signals and to improve the visual resolution, the Savitsky-Golay 2nd order derivations were also applied in this study. The selection of the models developed was largely dependent upon determination coefficient for the calibration set (R²) and validation set (r²), root mean squares error of calibration (RMSEC), the RPD (ratio of standard error of performance to standard deviation) and root mean squares error of prediction (RMSEP) used to evaluate the practicability of NIRs for determining to the quality indicators of sheepgrass hay.

The formula for RMSEC and RMSEP are described briefly below: ${\displaystyle \mathrm{RMSEC}=\sqrt{\frac{1}{n_{\mathrm{c}}}\sum _{\mathrm{i}=1}^{n_{\mathrm{c}}}{\left(\overline{y_{\mathrm{i}}}-y_{\mathrm{i}}\right)}^2}}$ ${\displaystyle \mathrm{RMSEP}=\sqrt{\frac{1}{n_p}\sum _{\mathrm{i}=1}^{n_{\mathrm{p}}}{\left(\overline{y_{\mathrm{i}}}-y_{\mathrm{i}}\right)}^2}}$ where $\overline{y_{\mathrm{i}}}$ is the predicted value of the ith observation; y_i is the measured value of the i-th observation; n_c is the number of observations in calibration set; n_p is the number of observations in validation set.

Chemical value

The obtained data (chemical value) were calculated and analyzed by the min, max, average, and standard deviation (SD), while the results of the calibration sets and validation sets were summarized in Table 1. These data sets indicated a wide variability in the chemical indices. For all the total sheepgrass sampled, the CP of 203 samples ranged from 6.2% to 14.33% (SD = 1.24), while the ADF ranged from 35.13% to 42.34% (SD = 2.33) and NDF from 50.71% to 71.08% (SD = 2.67). The datum of CP and NDF were nearly normally distributed around the average (10.54% and 60.89%) but the distribution of ADF measurements was skewed to the right, with the right tail stretching further than the left tail (mean = 38.74%).

This wide heteromorphosis in the quality indices of the sheepgrass hay was beneficial to successfully establish a relationship between the NIRs and the quality indices (Yu et al., 2010). Outliers in the collected sample were considered and thus excluded during the process of calibration. The elimination of outliers was based on the criterion that if mahalanobis distance for a sample was 3 SD or more (Yang, Han & Fan, 2006). For the value of chemical indices in this study, 2 outliers were removed for CP, 8 for ADF, and 1 for NDF (Table 1).

NIR spectra

Some peaks and valleys were obviously shown in the spectra (Fig. 2), which represent the characteristics of sheepgrass hay including hidden information of different components and its quantities. All the near infrared reflectance spectra of the collected sample were separated the two groups of spectra with different slopes, one group displayed an obvious absorbance peak at wavelengths of nearly 1, 100 nm, another group showed significant spectral peak was at approximately 1, 450 nm (Fig. 2). From the NIRs of the total sheepgrass samples (N = 203), the absorbance peaks were very overlapped, primary reason that the spectrum includes combinations and overtones of vibration such as stretching and bending of hydrogen-bearing functional groups such as –CH, –OH, and –NH (Saeys, Mouazen & Ramon, 2005). Meanwhile, other interference information affects the near infrared reflectance spectra (Meissl et al., 2008). Therefore, the quantitative evaluation are of difficulty through NIRs alone. Multivariate methods might be necessary to analyse the response of quality of sheepgrass from spectral characteristics with the support of chemometric methods, e.g., PLS, PCR and MLR analysis (Meissl et al., 2008; Saeys, Mouazen & Ramon, 2005).

Multivariate calibration analysis

For quantitative evaluation methods multiple linear regression (MLR), principal component regression (PCR) and partial least squares regression (PLS) were used to for data modelling and for predicting the investigated chemical constituents. For model validation, full cross-validation was used; each case was predicted by the model derived from all other remaining spectra. The results of full cross-validation is the average of the standard error of prediction values produced during each cross-validation step. The optimum number of Orthogonal factors (PC) for the models was obtained by using the leave-one-out cross validation technique for the calibration set. The recorded log (1/R) spectra were smoothed and transformed to second derivative before the analysis. The optimum models were achieved in Table 2 for the spectral processing method of CP, ADF, and NDF. PLS and PCR methods are based on the regression of the full spectra while MLR is based on discrete parts of the spectra (2 wavelengths) in this study. Table 2 showed that the comparison of the accuracy of the calibration results achieved by near infrared spectroscopy using PLS, PCR, and MLR methods, respectively, for evaluation the chemical indices of samples selected.

The results of the calibration procedures are presented in Table 3. It can be seen in the table that the smallest value for prediction error (RMSEC) with the highest coefficient of determination for the calibration set (R²) was provided by the PLS method. The results were less accurate for CP, NDF and ADF provided by the PCR method. The MLR had the worst prediction. Therefore, predictions were moderately successful for CP, NDF, and ADF was provided by the PLS method with second derivative of log(1/R) spectra.

The results of the cross validation are presented in Table 4. It can be seen in the table that the smallest value for prediction error (RMSEP value for CP, NDF and ADF is 1.24, 1.41 and 1.37, respectively) with the highest the coefficient of determination for the validation set (r² value for CP, NDF and ADF is 0.91, 0.90 and 0.90, respectively) was provided by the PLS method. The results were less accurate for CP (RMSEP = 2.63, r² = 0.86), NDF (RMSEP = 4.50, r² = 0.82) and ADF (RMSEP = 5.86, r² = 0.82) provided by the PCR method. The MLR had the worst prediction for CP (RMSEP = 6.74, r² = 0.74), NDF (RMSEP = 8.59, r² = 0.62) and ADF (RMSEP = 9.91, r² = 0.60). Therefore, predictions were moderately successful for CP, NDF, and ADF was provided by the PLS method. The application of PCR using the whole wavelength region requiring scanning type spectrometers as with PLS, resulted in an almost twice as high prediction error. The same order of magnitude could be obtained for RMSEP value by using MLR. Compared to PCR and PLS, the MLR model uses only two characteristic wavelengths for the calculation, creating the conditions for constructing cheaper single-purpose filter instruments.

Similarly as RPD determined by the above three methods in Table 4, the accuracy was considered good for RPD > 2, acceptable for 1.4 < RPD < 2, and unreliable for RPD < 1.4 according to Albrecht et al. (2009) and Chang et al. (2001). Therefore, the highest RPD was provided by the PLS method for CP(RPD = 3.2), NDF(RPD = 3.0), and ADF(RPD = 3.1) while the worst RPD provided by the PCR and MLR method are also shown in Table 4.