Applying SG to all our X matrix (Pracma Package)

"R" is without any doubt a great and wonderful community,  and it is nice to see how the package developers and maintainers help you in case you have any doubts.

It was the case some time ago when was writing some posts about the ChemoSpec package and Bryan Hanson helps me with some doubts. After the last post, I wrote a mail to Hans Werner (Pracma Package) , and he replied quickly, telling me the reason the "savgol" function use a vector indeed a matrix, and giving to me some ideas, to convert all the spectra matrix to Savitzky Golay.

Of course one of the ways is to use the apply function. When applying the SG filters there is  a reduction in the number of data-points at both sides of the wavelengths, depending of the window size.

 So I tried this way, to see all the spectra together:

> library(pracma)
# This script is for first derivative
> X1_sg1d_pracma<-apply(nir.training1$X,1,savgol,11,4,1)
> matplot(wavelength2[11:281],(X1_sg_pracma
+[11:281,]),type="l",xlab="Wavelength (nm)",
+ ylab="1/R (SG 1st derivative)",lty=1,col=1,main="SG-1D1104")

 

# This script is for second derivative
> X1_sg2d_pracma<-apply(nir.training1$X,1,savgol,11,4,2)
> matplot(wavelength2[11:281],(X1_sg_pracma
+[11:281,]),type="l",xlab="Wavelength (nm)",
+ ylab="1/R (SG 2nd derivative)",lty=1,col=1,main="SG-2D1104")

 

Savitzky Golay filters with Pracma Package

Pracma package has the function "savgol", where we can apply Savitzky Golay filter to a vector (in our case a spectrum).
The function is:
  savgol(T, fl, forder, dorder)

And the Arguments are:
 T... Vector of signals to be filtered.
fl... Filter length (for instance fl = 51..151), has to be odd.
forder... Filter order (2 = quadratic filter, 4 = quartic).
dorder... Derivative order (0 = smoothing, 1 = first der, etc.).

As you know I´m using in my last post the shoot-out 2002 data to develop a tutorial, and I read an article from the winner of this shoot-out where he use the shoot-out spectra with Savitzky Golay, a filter of 11, quartic, and second derivative using the Unscrambler software.

So I try thess values in the arguments of the Pracma SG filter and the results of the bands look exactly the same that the ones in the article, so this option looks good to work with this data. Anyway I will try also with the other functions from other packages.

In the case of the Pracma package we have to use a vector (a single spectrum), so some work has to be done to convert all the matrix of spectra, but the results looks great.

X1_sg_pracma<-as.matrix(savgol(nir.training1$X[1,],11,4,2))

PCAs with three diferent methods and projections (Test and Val Set)

The shoot-out 2012 is composed with 155 samples for the training set (Blue color), 460 for the Test Set (Red color), and other few samples for the Validation set (Green color).

I have developed 3 different ways of Principal Components Analysis and I would like to show you the score plot of PC1 vs PC2 developed with the Training Set and the projections on that space of the Test and Validation Set.

This first plot is in the case  of using PRCOMP:

Second case is using NIPALS for the calculation of the PCAs:

and third using SVD for the calculation of the PCAs

As you can see no differences, and we can have the conclussions that the Training Set cover the variability for the samples of the Test Set and Validation Set, so we don´t have to extrapolate outside the calibration space.

I will writte a post with the code (quite long) if interested. Let me know.

 

Adding Category Variables to a Data frame in R

Normally I used data frames to manage NIR data, the data frames are composed normally in my case by a X  or Spectra matrix (dataframe$X), and a Y or constituent matrix (dataframe$Y). But when we want to manage and understand plots, like score plots, it is interesting to classify the samples with some category variables.

This category variables can be: "location", "type", "customer", "product",.....

In the case of the shoot-out data the samples can be classify by their content of the main parameter, and can be classified as:

"Low"             (if the sample has less than 160 mg)
"Medium"          (between 160 y 221 mg)
"High"            (more than 221 mg)

Let´s create the variable in the data frame of the training set for instrument 1
nir.training1$type[Y <=160] <- "Low"
nir.training1$type[Y>160 & Y<221] <- "Medium"
nir.training1$type[Y>=221] <- "High"
Now we have a new variable in the data frame called "type"
Check it with:

names(dataframe)

and appart from X and Y we have Type.

We proceed the same way for the other dataframes.

Another thing is that we can create a big data frame with all the spectra from different instruments and sets and create a category variable for the instrument ( A and B), and another for the Set (Training, Test and Validation).

An introduction to the "´Resamble" package

Some time ago I wrotte a post : An introduction to the "prospectr" package . There I gave a list of some  Chemometric Packages for R. Some days ago in a NIRS forum gave the name of another Chemometric package, so the options to use Chemometrics with R are growing day by day.

This package is "Resamble". This is the link to the reference Manual.

 

More details at:

http://cran.r-project.org/web/packages/resemble/

http://l-ramirez-lopez.github.io/resemble/

The author explains that with this package algorithms such as LOCAL and locally weighted PLS regression can be easily reproduced.
I really want to test it as soon as I can.