Soil Sample Preparation

Modelling complex spectral data (soil) with the resemble package (V)

 Continuing with the vignette Modelling complex spectral data with the resemble package (Leonardo Ramirez-Lopez and Alexandre M.J.-C. Wadoux) 

Continuing from last post: Once we have calculated the PC terms or components (11 in the case of the last PCA analysis using the method OPC), we define planes defined by the combinations of two of those terms (for example: PC1-PC2, PC2-PC3, PC1-PC3,…), and the training spectra is projected on the plane to get the scores of every spectrum vs. each PC component. All those scores are kept in a score matrix “T”. All the projections form a cloud that in the case of just two terms would be a 2D cloud, making easy the interpretations of the distances between every sample and the mean or their neighbors. But in the case of more dimensions it is a multivariate cloud, making the visual inspection more difficult, so we have to check the projections individually in 2D planes or 3D planes.

Algorithms like the Mahalanobis distance to the mean or to the neighbors will help us to check if the sample can be an outlier, it has very close neighbors (so it is represented by samples in theory similar), or if the sample has not closer neighbors and is a good sample to improve the structure of the database and make it more robust.

Let´s see in the case of the previous code one of those score planes, the one formed by the PC1 and PC2 terms:

plot(pca_tr_opc$scores[,1],pca_tr_opc$scores[,2],

    xlim = c(min(pca_tr_opc$scores[,1]),

    max(pca_tr_opc$scores[,1])),
    ylim = c(min(pca_tr_opc$scores[,2]),

    max(pca_tr_opc$scores[,2])),

    xlab="PC1 ", ylab="PC2 ")

We can project the testing data on the same plane, getting the scores of the samples:

pca_projected <- predict(pca_tr_opc, newdata = testing$spc_p)

par(new=TRUE)

plot(pca_projected[,1],pca_projected[,2], col = "red", 

    xlim = c(min(pca_tr_opc$scores[,1]),

    max(pca_tr_opc$scores[,1])),
    ylim = c(min(pca_tr_opc$scores[,2]),

    max(pca_tr_opc$scores[,2])),
    xlab=" ", ylab=" ")

If we had only two PCs, this plane would be enough to show us the cloud, but we have 11  PCs. We can add graphically one more dimension (3D) and we see the cloud more clearly.

library(plotly)

T_training <- as.data.frame(pca_tr_opc$scores)
plot_ly(T_training, x=~T_training[,1], y=~T_training[,2],
             z=~T_training[,3], alpha = 0.7)

Practice yourself and rotate the plot to different angles.

Modelling complex spectral data (soil) with the resemble package (IV)

Continuing with the vignette Modelling complex spectral data with the resemble package (Leonardo Ramirez-Lopez and Alexandre M.J.-C. Wadoux) 

Now we will use the PCA with the method “opc” in order to find the optimal number of components, bases on the its rationale behind that if two spectra are close in the X space (near neighbors), their constituents values will be closer as well on its value, so the optimal number of components will be the one that makes minimum the RMSD (root mean square difference) between them.

For more details you can find more info from the developers of this algorithm : L. Ramirez-Lopez, Behrens, Schmidt, Stevens, et al. (2013) 

optimal_sel <- list(method = "opc", value = 40)

pca_tr_opc <- ortho_projection

              (Xr = training$spc_p,
               Yr = training$Ciso,
               method = "pca",
               pc_selection = optimal_sel)

pca_tr_opc     # to obtain details of the PCA calculations.

 

We specify a maximum value of 40, and the “opc” method estimate that 11 is the best option. If we plot it, we can see graphically the reason:

The vignette shows an interesting code, that if you run it will get the XY plot of the reference Ciso value (for every spectrum) and the reference Ciso value for its closer neighbor, and we can se a high correlation what is really the idea behind the “opc” method.

I found these option very interesting, so we will continue exploring the vignette that sure will help for the purpose of its title.

Modelling complex spectral data (soil) with the resemble package (III)

Continuing with the vignette Modelling complex spectral data with the resemble package (Leonardo Ramirez-Lopez and Alexandre M.J.-C. Wadoux) 

(from previous post) - These 825 samples are divided in two sets, one for training (value equal 1 in the variable train) and another for testing (value = 0). 

 NIRsoil %>%
    count(train)
    train       n
     0         207
     1         618

Let´s create these two dataframes:

training<- NIRsoil[NIRsoil$train == 1, ]
testing<- NIRsoil[NIRsoil$train == 0, ]

Now the vignette explain how to proceed with the dimensionality reduction and the methods we can use for that purpose: pca, pca_nipals and pls.

You can practice each one (as in the vignette code) using the default values for the explained  variance  for every component  (by default > 1%) or you can fit the cumulative variance you want.

PCA method use the Singular Value Decomposition algorithm, PCA_NIPALS, the NIPALS algorithm and PLS use (apart of the predictor variables (absorbances)) the response variables , maximizing the covariance of latent variables from  predictors with the response (we have to specify which response variable to use).

Just follow the vignette examples for these three methods, and check how many components are recommended on each one.

For a more advanced use of the methods we can configure the "pc_selection" method, generating a list with the method ("var", "cumvar", "manual" or "opc"), and the value of the variance each component  must explain ("var" case). In the case of "cumvar" we set the value to the total cumulative variance we want to explain, and we will see in more detail the case of the "opc" in another post.

Take into account that we will create the ortho-projections with the training set, and after, we can project new data (testing set) on the planes created with the training set, to get the scores of this new testing data.

Modelling complex spectral data (soil) with the resemble package (II)

Continuing with the vignette Modelling complex spectral data with the resemble package (Leonardo Ramirez-Lopez and Alexandre M.J.-C. Wadoux), now it is time to see the predictor variables which are reflectance values of the soil samples acquired in a NIR (Near Infrared Reflectance) instrument in the range from 1100 to 2498 nm in two nm steps, so we have 700 data points. We prepare a vector with the wavelengths and we call it "wav" (same as the vignette).

wavs<-NIRsoil$spc %>% colnames() %>% as.numeric()

Now we can se to the raw spectra (spectra without any treatment):

matplot(x = wavs, y = t(NIRsoil$spc),  

        xlab = "Wavelengths, nm",
        ylab = "Absorbance", type = "l", 

        lty = 1, col = "#5177A133")

Now it is time to treat the spectra. In the vignette, the spectra is reduced in the number of data points and treated with a first derivative (Savitzky Gollay) using a first polynomial order and a window of 5. This reduction and signal improvement techniques are very useful to prepare the spectra "X" matrix for new reduction techniques after, saving this way computation time.

We just have to use the code of the vignette, but of course we are free to use other mathematical treatments to reduce the scatter effects  or improve the resolution.

NIRsoil$spc_p <- NIRsoil$spc %>% 

  #we make a reduction of the number of data points 

resample(wav = wavs, new.wav = seq(min(wavs), 

         max(wavs), by = 5)) %>% 

  #and apply to the spectra the Savitzky Golay function

  #polynomial order =1

  #window = 5

  #first derivative

  savitzkyGolay(p = 1, w = 5, m = 1)

Let´s create a new vector considering the wavelength reduction:

new_wavs <- as.matrix(as.numeric(colnames(NIRsoil$spc_p)))

and plot the spectra to see their appearance:

matplot(x = new_wavs, y = t(NIRsoil$spc_p),
        xlab = "Wavelengths, nm",
        ylab = "1st derivative",
        type = "l", lty = 1, col = "#5177A133")

Now in the data frame "NIRsoil" we have two spectra matrices, the raw spectra (spc) and the spectra reduced and math treated with the SG first derivative (spc_p).

We can check the dimensions of these matrices:

names(NIRsoil)
    "Nt" "Ciso" "CEC" "train" "spc" "spc_p"
dim(NIRsoil$spc)
    825 700
dim(NIRsoil$spc_p)
    825 276

In the next post we will continue the preprocessing process and  preparation of the data as the vignette suggest, trying to understand the different procedures to model, as better as possible, the soil spectral data.