Practicing with "NIR hyperspectral image analysis using R" - part 3 .(A.A.Gowen - NIR News)

Today I have been practicing with  part 3 of these nice tutorials about NIR hyperspectral image. The author (@eefieg) leaves some parts without code in order that the reader practice and get the same or similar results. For example compare two classification images obteined with two different methods, so if we substract one from the other, in the best of the cases, we would not see any frames, just a uniform color (the color which represents the zero value).

Finally I got (I think) the same residual image as Aoife, in the article.

 

r_PC1<-Im1_PC1[,]<= -2
r_PC1<-1*r_PC1         

g_PC1<-(Im1_PC1[,]> -2)&(Im1_PC1[,]<=0.7)
g_PC1<-3*g_PC1       

b_PC1<-(Im1_PC1[,])>0.7
b_PC1<-2*b_PC1       

col_PC1<-r_PC1+g_PC1+b_PC1
dev.new();
image(col_PC1)

diff<-col_PC1 - Im1_Class
dev.new();
image(diff)

An introduction to the "prospectr" package

prospectr: An interesting package that we will use to apply to our tutorials and data. I like that it has an algorithm (SELECT) used in Win ISI and the name of the function is "shenkWest" (Shenk and Westerhaus, develop the SELECT algorithm for Win ISI), so we have more tools to use R in spectroscopy.

Don´t forget that we have also other packages like:

ChemoSpec

hyperSpec

cluster

mvoutlier

pls

signal

soil.spec

caret

...
You will find  PDF tutorials for these packages on the Web, and for others like ChemometricsWithR and chemometrics, there are even some books.

 

Comparing regressions with different scatter math-treatments (Shootout 2002 Tutorial)

In a previous post, we have developed the calibration without any math-treatment with the Training Set from Instrument 1, without any treatment , knowing that it was not the best choice, and we look to the LOO (leave one out) cross validation errors to check the performance.

Now we develop the regressions with some anti-scatter math-treatments to compare the cross validation errors and to decide which of them performs better. Anyway the shootout supplies also a test set, so validating with this test set will give us a better idea about how the calibration is performing with independent data.

 

First thing to do is to convert the “X1. Training” and “X1. Test” matrix to the math treatment we want to use: SNV (Standard Normal Variate), Detrend , SNV + Detrend and MSC.

>nir.train1_snv<-data.frame(X= I(X1_snv),Y=I(Y))

>nir.train1_detrend<-data.frame(X= I(t(X1_detrend)),Y=I(Y))

>nir.train1_snvdt<-data.frame(X= I(t(X1_snvdt)),Y=I(Y))

>nir.train1_msc<-data.frame(X= I(X1_msc),Y=I(Y))

Now we can develop the PLS regressions:

 

>mod1_snv<-plsr(Y~X,data=nir.train1_snv,+

 ncomp=10,validation="LOO")

>mod1_detrend<-plsr (Y~X,data=nir.train1_detrend,+

 ncomp=10,validation="LOO")

>mod1_snvdt<-plsr(Y~X,data=nir.train1_snvdt,+

 ncomp=10,validation="LOO")

>mod1_msc<-plsr(Y~X,data=nir.train1_msc,+

 ncomp=10,validation="LOO")

We can plot the RMSEP values versus the number of components (or terms), to have a better idea of the performance of the models (black line is the model without mat-treatments or raw spectra, green line is with just Detrend, the rest (SNV, SNV+DT and MSC) are almost overlaped.

 

But if we want to see it with more details, we have to see the numbers provided by the summary of the models. I mark in yellow the smallest values.

 

 

But the question can be: Do we have to choose the number of components which gives the small RMSEP?.

We will reply to this question soon.

In a next post we will do the same with derivative mixed with scatter corrections to see if we get better values for RMSEP and we will check it with an external validation (don´t forget that these RMSEP are for Cross Validation).

 

 

 

 

 

 

SNV and Detrend (Shootout 2002 NIR with R)

I want to say “thanks” to Aoife Gowen (@eefieg) for the references, in her last article to this blog. The article is the part 4 of a series called “NIR hyperspectral image analysis using R” and the title of this article is: "Pretreatments and partial least square discriminant analysis". (NIR News Vol.25 Nº 5 August 2014).
 

Lately I am playing with the Shootout 2002 data so we can summarize on this data things that we have seen before in the blog.
 

In the last post, I was using the raw spectra and that is not (in some cases) the best way to develop a calibration. So we treat the raw spectra with math treatments to improve the correlation of the spectra with the reference values, and to make more transferable calibrations. If we use a correct math-treatment the regressions will use less terms, and will improve the statistics, and the transferability.

Now I write the script used to treat the matrix X1 (training spectra matrix in Instrument 1):

>par(mfrow=c(2,2),ps=14)
>matplot(indices,t(nir.training1$X),type="l",col=1,

 + main="Raw Spectra")
#Applying "SNV" to the Shootout 2012 Data>X1_snv<-scale(t(nir.training1$X),center=TRUE,scale=TRUE)
>matplot(wavelength2,X1_snv,type="l",

 + xlab="Wavelength (nm)",ylab="1/R",lty=1,

 + col=2,main="SNV math_treat")
#Applying "Detrend" to the Shootout 2012 Data

>library(pracma)
>X1_detrend<-detrend(t(nir.training1$X))  
 +matplot(wavelength2,X1_detrend,type="l",

 +lty=1,xlab="Wavelength(nm)",ylab="1/R",col=3,

 +main="Detrend math_treat")
#Combining SNV with Detrend

>X1_snvdt<-detrend(X1_snv,tt="linear")

>matplot(wavelength2,X1_snvdt,type="l",lty=1,

 +xlab="Wavelength(nm)",ylab="1/R",col=4,

 +main=" SNV+Detrend math-treat")

 

As we are writing this script, we see on the Active Graphics Window the spectra , so we can compare them.

These are some of the most familiar anti-scatter math-treatments, apart from others like MSC.

 

We will continue with this type of exercises in the next posts.

XY Plot to compare predictions (same model / diferent instruments)

As we mention several times, plots help us to understand better what is happening. In the previous posts we are working with the shootout 2002 data in order to work with as many tools as possible we have in R to understand this data set and to make a model we can use to predict as better as possible a validation set.
Plotting the predictions of the Test set acquired in Instrument 2 with the model developed with the Training set of samples acquired in Instrument 1, versus the predictions of the test set acquired in instrument 1, can give us an idea that there is something more than a simple bias, and that we have to take into account a group of samples that behaves differently.
Lets see that plot:
1) We make the predictions of the Test1 set with Model 1 and stored it in:
  anl_tst1_mod1<-predict(mod1,ncomp=3,
+ newdata=data.frame(Y=I(Y.test),X=I(X1.test)))
2) We make the predictions of the Test2 set with Model 1 and stored it in:
  anl_tst2_mod1<-predict(mod1,ncomp=3,
+ newdata=data.frame(Y=I(Y.test),X=I(X2.test)))
3) Now we plot the predictions versus reference for anl_tst1_mod1
  plot(Y.test,anl_tst1_mod1,col="blue",ylim=c(150,240))
4) Now we write the script:
   par(new=TRUE)
5) Now we plot the predictions versus reference for anl_tst2_mod1
  plot(Y.test,anl_tst2_mod1,col="red",ylim=c(150,240))
6) Now we get this plot:

We can see that most of the samples behave similarly with a bias

 

If interested to follow this tutorial see also:

Working with the Shootout 2002 (NIR) in R

Plotting the predictions (Working with the Shootout 2012 NIR data in R)

Comparing Spectra and Coefficients (Shootout 2002 NIR with R)