PCR vs PLS (part 1)

The calculation of the Principal Component Regression coefficients can be obtained by the Singular Value Decomposition of the X math treated matrix, where we obtain the scores of the samples, multiplying the "u" and "d" matrices (u %*% d is what we call often the T matrix), and after, we make a regression of the constituent (Y matrix) vs. the scores (used as independent variables) to obtain the regression coefficients of the scores.

 

After we multiply the regression coefficients of the scores by the loadings "v" (P transposed matrix) to obtain the regression coefficients of the principal components regression.

 

This is the long way to do a PCR, but the quicker option is to use the "pcr" function from the R "pls" package.

## The short way:
library(pls)
Xodd_pcr2<-pcr(Prot[odd,] ~ X_msc[odd,],ncomp=5)
matplot(wavelengths,Xodd_pcr2$coefficients[,,5],
             lty=1,pch=NULL,type="l",col="red",
             xlab="Wavelengths",
             ylab="Regression Coefficients")

We can look to the summary:

summary(Xodd_pcr2)

Data:  X dimension: 327 100 
 Y dimension: 327 1
Fit method: svdpc
Number of components considered: 5
TRAINING: % variance explained
             1 comps  2 comps  3 comps  4 comps  5 comps
X             99.516    99.87    99.97    99.99    99.99
Prot[odd, ]    3.934    15.34    19.31    65.45    68.17

We have develop the PCR models with the half of the total database which we has split in two (odd and even) in the post "Splitting spectral data into training and test sets".

As we can see with just 1 PC component we explain more than 99,5% of the variance in X, but just a few of the Y variable (Protein in this case). We need to add more than just the first term (at least 4) to succeed in the calibration, but the cost is that we can add noise and over fit the model.

 

The external validation (even samples) set will help us to see the performance of the model.

 

 

 

Splitting spectral data into training and test sets

It is common to split the spectral data into a validation or test set and a calibration or training set. This can be done in different ways (random, structurally,...), selecting different percentages for each.

 

This is a simple case and we are going to select 50% of the samples for the calibration or training and the rest (the other 50%) for validation or test. This way we have a Training Set and a Validation test.

 

One simple way to proceed is to sort the samples randomly or structurally (sort by constituent value, date of acquisition, type of product,....), and select the odd samples for the Training Set and the even samples for the Test Set.

## DIVIDE DATA SETS INTO CALIBRATION AND VALIDATION SETS
##We create a sequence with the odd samples
odd<-seq(1,nrow(X_msc),by=2)
##We create a sequence with the even samples
even<-seq(2,nrow(X_msc),by=2)
#We take the odd samples for the training set
X_msc_tr<-X_msc[odd,]
Prot_tr<-Prot[odd,]
#We take the even samples for the validation set
X_msc_val<-X_msc[even,]
Prot_val<-Prot[even,]
matplot(wavelengths,t(X_msc_tr),type="l",xlab="wavelengths",
        ylab="Absorbance",col="blue")
par(new=TRUE)
matplot(wavelengths,t(X_msc_val),lty=1,
        pch=NULL,axes=FALSE,
        type="l",col="red",xlab="",ylab="")

We can see in the plot the training spectra in blue and the test spectra in red. We will continue practicing with these sets in the next posts.

 

 

Building a neural network from scratch in R

 We are use to work with lineal models, but sometimes we see that our models are not so linear, so we can apply logs or other algorithms as LOCAL to model the data. But another option is to apply a Neural Network Model.
Working with NIR data  from many different sources is not easy and we have to check the data, make some groupings and to create different products in order to take certain ranges, eliminate bias and curvature, check for different math treatments to see which one makes the model more linear, etc.

Neural Networks Models deal with all these but is more complex and unknown. It requires also for a lot of data.

But how it works?: Read this link to know how it works and how can you build these models with R.

 

 

 

Using R to prevent food poisoning in Chicago

Amazing to see the applications of R, in this case for food safety, see the next video to see how it is applied:

 

You can see more details in the post from Revolutions blog:

Using R to prevent food poisoning in Chicago

 

This is a talk of the author of the R package used for this application

Help of category variables to understand the spectral population

When we work with a data set, it is important to know all the information we can compile about it, so we can create category variables in order to understand better the sample population. In the soy meal data set, I know that some of them come from Brazil, others from USA and the rest... I don´t know where they can from (probably from the same countries, but the samples were not labeled when acquired). So I created a category variable called "Origin" with a group 1 (from Brazil), a group 2 (from USA) and a group 3 (origin unknown).

 

In "R", 1 correspond to black color, 2 to red and 3 to green. We can plot the Mahalanobis distance ellipse in PC1 and PC2 and to see how the samples are grouped.

 

drawMahal(T_msc,center=apply(T_msc,2,mean),
          covariance=cov(T_msc),quantile=0.975,
          col=soy_ift_conv$Origin,
          xlab="PC1",ylab="PC2")
legend("topleft",legend=c("Brazil", "USA","Unknown"),
       col=c("1","2","3"),pch=1, cex=0.8,
       title="Origin")

It´s a pity not to have all the information about the green samples. Anyway consider also another types of categories variables, for example order the samples by the constituent value and add a category for high, medium and low protein.