Reviewing some posts (Fatty Acids with R)

Thanks to Christof for his e-mail, and their interest in this blog. I know, like Christof ,other people would like to start the Tutorial about the NIT: Fatty Acids with R. 

I have sent in the past the "tocino6.txt" file by mail to some blog readers, but maybe I had not read (or I lost) some e-mails asking for the file, so if this is the case, please let me know it again and I will send you the file by e-mail or trough a link to a  Dropbox folder, where I will add the data of the tutorials or other information.

I deleted, also in the past some pictures accidentally, so I will review first the tutorials of the "NIT:Ftty Acids with R", and when ready I will write a post to the link of all the posts of these tutorial. Also a new label will be added to be easier to find them.

By the moment, the first three posts of this tutorial are ready with new features scripts and pictures, to everyone who wants to follow it.

NIT: Fatty acids study in R - Part 001
NIT: Fatty acids study in R - Part 002
NIT: Fatty acids study in R - Part 003
NIT: Fatty acids study in R - Part 004
NIT: Fatty acids study in R - Part 005


You would need the tocino6.txt file from post 1:5, after these posts we will validate the model eith the file "tocino6val5.txt" and we combine both in "tocino7.txt" to develop the final regression.
These files are available from Dropbox or I cand send tem by mail, just let me know your interest writing to me a mail: cuesta_joseramon@yahoo.es
 

NIT: Fatty acids study in R - Part 007

Once we have chosen the model, we can continue acquiring spectra of new samples. Spectra is exported to a txt or csv file and we imported in R to be reprocessed.
We use the function “predict” from the PLS package. I have done this with 20 new samples. We need first to apply to them the adequate math treatment (same as the used in the model). I call this sample set for prediction “fatt2ac_val”, after apply the “msc” math treatment.
So, let´s see the predictions:

> predict(C16_0,ncomp=12,newdata=fatty2ac_val)

, , 12 comps

    C16_0

220   22.01807

221   20.44803

222   19.79991

223   21.64232

224   20.29058

225   20.20099

226   21.52053

227   19.83305

228   18.95492

229   21.39239

230   21.11044

231   20.67454

232   19.28662

233   20.97292

234   21.70614

235   20.27464

236   19.70897

237   21.30686

238   20.21069

239   19.21576

I have used the model with 12 terms.
If this data has more than the spectra, (the Lab values) we can also validate and to check the number of terms to use.

> RMSEP(C16_0,newdata=fatty2ac_val,ncomps=12)

(Intercept)      1 comps      2 comps      3 comps      4 comps      5 comps 

     1.6280       1.5987       1.2554       1.1071       1.3447       0.9122 

    6 comps      7 comps      8 comps      9 comps     10 comps     11 comps 

     0.8754       0.8413       0.6597       0.6154       0.5669       0.5791 

   12 comps     13 comps     14 comps     15 comps     16 comps 

     0.5935       0.6261       0.5994       0.6315       0.5787

We can see our RMSE error and compare it with the RMSEP obtained in the PLSR “LOO validation statistics” which was 0,5733.We can see also that we would get even lower values for validation with a lower number of components (10).
It seems that the is working almost as expected, but let´s have a look to the plots:

>predplot(C16_0,ncomp=7:15,newdata=fatty2ac_val,

+ asp=1,line=TRUE)

I can see a Bias, of almost 0,50.

The error corrected by the Bias would be 0.34.

Bias can be due to diferent reasons (temperature, sample presentation,particle size,optical path,state of the instrument,...).

This samples are fat triturated as a paste, and put it in a petri dish, and place it in an instrument in transmitance.

The next step would be to add this data to the data base and recalibrate.

Divide the data into a Calibration and a Validation Set. Be sure that in the validation set there are some samples of this last set.

The idea with all this series of post related to this data set was to work with data I used almost daily in my job, and I wanted to see how to proceed in R, and once you get use to it the results are very good for the understanding of chemometrics for multivariate data. I will continue exporting some of my data sets to work in R. The analysis of Fatty Acids by NIR/NIT can give good results for some of them (C16:0, C18:0, C18:1,..).

NIT: Fatty acids study in R - Part 006

In one of the columns, for constituent C16_0, one sample (57) has a value of “zero” (we could see this in the histogram).The reason for that is that the laboratory did not supply this value. The PLS regression will consider the lab value as cero, so we will get a plot like this:

I observed also that the sample 219 has a high residual for the regressions of all the constituents, so I decided to remove these two samples from the sample set in order to continue, and to develop the models.
I am starting with R, so I will appreciate if you add comments in order to do this task in a simpler way.
I create two sample sets, in order to remove these two samples (219 and 57):

> fattyac1<-fattyac_msc[1:56,]

> fattyac2<-fattyac_msc[58:218,]

and I combined this three sets again:

> fattyac_msc1<-rbind(fattyac1,fattyac2)

Well, I can develop my regression now:

Now we have to take the decision of how many terms to choose. Let´s see the validation plot with 7 and 12 components (terms).

plot(C16_0,ncomp=7,which="validation")

abline(0,1,col="red")

plot(C16_0,ncomp=12,which="validation")

abline(0,1,col="red")

It is clear that the decision to choose one model or the other will have a great influence in the predictions. We need a validation set to make a better decision. But I think that it will work better with 12 terms.
It will be important, if possible to find samples with C16:0 values bellow 18 to add to our database in order to develop a better model.
Another decision could be to keep out this extreme sample until we find more, but we can decide to leave it, in order to extrapolate better in this zone.
It is important not to have unique samples in the model. In this case we have one. We have to consider this.

If you want to follow this tutorial, please send me an e_mail. I´ll send you the “txt” file attached.

NIT: Fatty acids study in R - Part 005 (Calibration)

There are several algorithms to run a PLS regression (I recommend to consult the books: “Introduction to Multivariate Analysis in Chemometrics - Kurt Varmuza & Peter Filzmozer” and “Chemometrics with R – Ron Wehrens”).
We are going to use the PLS package, and we are going to develop, maybe the constituent which looks more promising: Oleic Acic (C18:1).
Of course we are going to use the MSC pretreatment. For vcross validation we are going to use “leave one out”.
I decide a maximum of 16 PLS terms.

[1] C18_1_reg<- plsr(C18_1~NITmsc, ncomp = 16,data =fattyac_msc, 
    + validation = "LOO")

Looking to the Cross Validation statistics, it seems that 12 is the best number of terms to use, anyway lets see the plots.

par(mfrow=c(1,2))
plot(C18_1_reg,"validation",estimate="CV")
abline(v=12,col="red")
plot(C18_1_reg,ncomp=12,"prediction")
abline(0,1,col="red")

We can see how after Term 12 the RMSEP increase.
We can see in the XY plot, how we have a few samples (like the sample 219) with a high residual (probably wrong lab value in some cases). If you prefer to see the residual plot, where you can see sample 219 in the upper right corner:

If we check the residual list (C18_1_reg$residuals),the extreme sample 66, fits well in the model, so we decide to keep it. No we can remove sample 219 from the sample training set removing this row from the data frame.
C18_1_reg<- plsr(C18_1~NITmsc,ncomp = 16,
   + data =fattyac1_msc,validation = "LOO")
We repeat steo [1] again, but thid time we change "fattyas_msc" for "fattyas1_msc", and look to the plots and statistics again. In this case I overplot the plots with the outlier (red colour) and without the outlier (Blue):

  

 We keep this last model and we will test it in the future with new samples (independent test set).


If you want to follow this tutorial, please send me an e_mail. I´ll send you the “txt” file attached.

NIT: Fatty acids study in R - Part 004

We already know that sample 66 has a MD of 11.6. Let´s see the values for the six constituents for this sample:

> fattyac_msc[66,1:6]

   C16_0  C16_1  C18_0  C18_1  C18_2  C18_3

66  15.8     2     6    62.3   10.2    0.6

Let´s compare with the summary

> summary(fattyac_msc)

        C16_0           C16_1                 

 Min.   : 0.00   Min.   :1.500     

 1st Qu.:20.10   1st Qu.:2.000     

 Median :21.00   Median :2.200     

 Mean   :21.34   Mean   :2.267     

 3rd Qu.:22.90   3rd Qu.:2.500     

 Max.   :26.00   Max.   :3.500    

   

  C18_0            C18_1     

Min.   : 5.800   Min.   :43.80 

1st Qu.: 8.600   1st Qu.:51.95 

Median : 9.400   Median :54.50 

Mean   : 9.711   Mean   :53.93 

3rd Qu.:10.500   3rd Qu.:56.15 

Max.   :14.000   Max.   :62.30 

 C18_2            C18_3      

 Min.   : 5.500   Min.   :0.3000 

 1st Qu.: 7.600   1st Qu.:0.5000 

 Median : 8.500   Median :0.6000 

 Mean   : 8.503   Mean   :0.6032 

 3rd Qu.: 9.100   3rd Qu.:0.7000 

 Max.   :14.700   Max.   :1.3000  

Sample 66 has the higher value for C18:1 (oleic acid), but it is not isolated in the histogram. For some reasons this sample differs from the others especially from 100 to 1050 nm. We will wait forward to take a decision about this sample.
Until now we have been managing with the X matrix.
Now we start to study the Y matrix. First thing to do is to have a look to the summary, and of course to the histograms.
If you want to follow this tutorial, please send me an e_mail. I´ll send you the “txt” file attached.

> hist(C16_0,col="red")

> hist(C16_1,col="blue")

> hist(C18_0,col="green")

> hist(C18_1,col="brown")

> hist(C18_2,col="violet")

> hist(C18_3,col="orange")

NIT: Fatty acids study in R - Part 003

We are going to use the spectra treated with MSC to continue with this tutorial, and as we saw in the first post, the spectra looks like this:

One way to see the variability in the spectra is to center the spectra, we only have to substract the mean spectra from every spectrum.

One easy way to do it is:

NITmsc_c<-scale(fattyac_msc$NITmsc,
center = TRUE,scale = FALSE)
par(mfrow=c(1,1))
matplot(wavelengths,t(NITmsc_c),lty=1,pch=NULL,
type="l",xlab="wavelengths",ylab="Abs") 

Another way to see the variabiñity in the spectra is the standard deviation spectrum, in this case we calculated it (in red color) and overplotted over the previous plot of centered MSC spectra.
sdNITmsc_c<-as.matrix(apply(NITmsc_c,2,sd))
par(new=T)
matplot(wavelengths,sdNITmsc_c,lty=1,pch=NULL,
  type="l",col="red",lwd=2,xlab="wavelengths",ylab="Abs")

 

Looking to the loadings:

We saw the scores (T matrix) in the previous post, let´s have a look to the loadings (P matrix).

>  fatmsc_nipalsP4pc<-fatmsc_nipals$P[,1:4]

> loading1<-fatmsc_nipalsP4pc[,1]

> loading2<-fatmsc_nipalsP4pc[,2]

> loading3<-fatmsc_nipalsP4pc[,3]
> loading4<-fatmsc_nipalsP4pc[,4]

It is interesting to find interpretation to these plots, for this is important to look carefully to all the plots (scores,loadings, original spectra,…), and also to consult bibliography about the position of the NIT bands for C-H, O-H,.., bonds.
In this case MSC does not remove totally the scater, so scatter continues beeing a high source of variance represented by the first PC.But once is removed we can see better bands in the other loadings.
Lets find the peaks in the loadings:
findPeaks(loading2)   #X878 X932 X972
findPeaks(loading3)   #X904 X944
findPeaks(loading4)   #X898  X916  X966  X998 X1004 X1010 X1028 
Let´s mark the most relevants over the loadings plots

  
  If you want to follow this tutorial, please send me an e_mail. I´ll send you the “txt” file attached.

NIT: Fatty acids study in R - Part 002

> library(chemometrics)

> fatmsc_nipals<-nipals(fat_msc,a=10,it=160)

> CPs<-seq(1,10,by=1)

> matplot(CPs,t(fatmsc_nipals$T),lty=1,pch=21,

  + xlab="PC_number",ylab="Explained_Var")

In the 2D plot, we can see that with 3 or 4 principal components, almost all the variance is explained. We see also how samples are well projected over the first PC, but how one sample seems to be a clear outlier when projected over the second PC.
Also another (or the same) sample seems to be an outlier when projected over the 4th PC.

  A look to the PC planes will show us the distribution of the samples, for all the combinations of the first four Principal Components.
> pairs(fatmsc_nipals$T[,1:4],col="red")

Let´s calculate the Mahalanobis distances to study better our sample population and to find outliers:

> fatmsc_nipals4pc<-fatmsc_nipals$T[,1:4]

> Moutlier(fatmsc_nipals4pc,quantile =0.975, plot=TRUE)

$md

  [1]  3.1025142  1.6787234  1.1155225  1.7130001  2.6460755  1.3341133

  [7]  1.6411590  1.8790763  0.8541428  1.1332246  2.3236420  2.8866269

 [13]  1.3146715  1.2759619  1.4875943  1.0603605  1.4880568  0.9665029

 [19]  1.8239722  1.6562767  1.9359163  1.6708538  2.2936930  1.4845431

 [25]  2.6051294  1.1436901  0.3571686  2.0533289  1.4248877  0.7107637

 [31]  1.2504387  2.0814050  1.8502729  3.1433570  2.2836840  0.3547710

 [37]  2.9424866  2.6208234  0.6914464  2.9649615  1.7914432  2.2730964

 [43]  3.0530848  1.1569603  1.5458923  1.3590878  1.9744677  1.8299434

 [49]  1.6564926  2.7850876  2.9147344  2.9858931  2.0337672  2.6220121

 [55]  2.4169714  0.9873321  2.7614810  3.4578931  4.3510717  2.1840045

 [61]  3.4219424  2.6471133  2.5050841  1.6500068  2.2036638 11.6174968

 [67]  3.1637271  3.0938694  2.3489142  2.5605777  1.7651892  0.7602064

 [73]  0.8169349  1.1276683  1.0530317  0.9008947  1.5501520  0.8291586

 [79]  1.8831524  0.5590048  2.3312774  2.0025709  0.7148548  1.9298735

 [85]  1.7581300  1.9388953  1.4556749  2.0408671  1.7715642  2.5011261

 [91]  4.4534119  2.8088303  4.2640203  0.9677583  6.1127505  1.2239764

 [97]  5.9621142  0.9987361  1.5365592  0.8917701  0.6152401  0.8996054

[103]  1.8370282  1.3580873  0.7873400  0.9220825  1.8619488  1.9298884

[109]  1.4912294  0.9832971  0.9842641  1.2018128  0.7935046  0.8925428

[115]  1.2003102  1.4462257  1.2691323  1.8269249  1.2838734  0.9981628

[121]  1.9145605  1.7954542  1.5230153  1.3347716  1.1156095  1.5871748

[127]  1.4889242  1.1780966  1.4165463  1.0057897  1.6742841  1.7999796

[133]  1.2231126  1.3167038  1.7676869  1.7475316  1.5718934  0.7844088

[139]  0.7250911  0.8394164  0.9434329  1.3583476  0.9143295  1.5666855

[145]  0.8250539  0.5027369  1.6273106  1.8940848  0.8493707  1.4611669

[151]  0.3644340  0.7813530  1.6332761  1.0557438  1.2848675  1.0695355

[157]  1.7891441  0.6474083  0.8387371  0.9655893  1.6508979  1.4765710

[163]  2.6846350  1.9820580  2.0689903  1.5834826  1.2542036  0.8494160

[169]  1.3529783  0.8451586  1.6718654  2.5892144  1.3678979  1.4070544

[175]  1.3870741  1.2010282  1.3446915  1.4648297  1.4599712  1.5161282

[181]  1.2140609  2.1280737  1.1751724  1.5939065  0.8337121  1.0548981

[187]  1.2061079  1.1519596  1.4011917  1.1339365  1.3009569  1.1758361

[193]  0.9313623  0.9973675  1.3783733  1.3145118  1.4065661  2.2898204

[199]  1.3149368  1.6195627  1.3458978  1.1028901  1.5325457  1.4918670

[205]  1.6747645  1.0730898  1.3003462  2.2767521  1.2188084  1.4188156

[211]  1.2551781  1.1094945  1.7552917  1.6537534  1.0851287  1.1067528

[217]  1.4062079  1.6325028  2.0682626

$cutoff

[1] 3.338156

We can see that sample 66 has a MD>11, so we can suspect which sample is looking to the plots, but you can selected and overplotted to be sure:
matplot(wavelengths,fat_msc[66,],lty=1,type="l",
   col="red",lwd=3,pch=NULL,xlab="nm",ylab="abs")
par(new=TRUE)
matplot(wavelengths,t(fat_msc),lty=2,type="l",
    pch=NULL,xlab="nm",ylab="abs")

If you want to follow this tutorial, please send me an e_mail. I´ll send you the “txt” file attached, or a dropbox link.