Tidyverse & Chemometrics (part 4)

It is important to know the correlation between the constituents of a product (in this case the fish meal), so we merge the "dry matter",  "protein", "fat" and "ash" into a Y matrix and calculate the correlation.

y<-as.matrix(cbind(fish1$DM,fish1$PROTEIN,
                   fish1$FAT,fish1$ASH))

Now let´s give a name to the columns of the matrix:

colnames(y)<-c("DM","PROTEIN","FAT","ASH")

The best way to check the correlation is in a matrix form so let´s load the package "coorplot":

library(corrplot)
corrplot(cor_y,method="number")

and we get this plot:

We can see that there is a high negative correlation between "ash" and "fat", and a very low (almost cero) correlation between "protein" and "dry matter".
Let´s check this in a XY plot:

library(ggpubr)
ggscatter(fish1_df, x = "fat", y = "ash",
          add = "reg.line", conf.int = TRUE,
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "FAT", ylab = "ASH")

ggscatter(fish1_df, x = "protein", y = "dm",
          add = "reg.line", conf.int = TRUE,
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "PROTEIN", ylab = "DRY MATTER")

To know this correlations it is important to interpret better the loadings, regression coefficients ,...

Tidyverse and Chemometrics (part 3)

Once that we have defined the math treatment that we consider a good option for our spectra, we can check visually if we find some outliers, and that was not the case in our second derivative spectra for fish meal, so we continue checking for the samples in the principal component space.

For this calculation I recommend to load the library "chemometrics" and to run the "nipals" algorithm, where we select the X matrix (in our case treated with the second derivative) and select also the number of PCs we want to check.

library(chemometrics)
nir2d_pc<-nipals(nir_2d,10)
#Two matrices are generated:
# T = Scores Matrix
# P = Loadings Matrix

Now we are going to work with the T matrix and we will see how it describes ellipses in the PC space which is a curious characteristic of the NIR spectra. We will define the limits for those elipses based on the Mahalanobis distance.

Chemometrics packages has the function "drawMahal", where we use the T matrix calculated with Nipals to draw the ellipses. There is an ellipse per two PCs, which i what we called PC scores map, but in general an ellipsoid is drawn in all PC space.

Now I write the code to check the maps PC1 vs PC2, PC1 vs PC3 and PC2 vs PC3:

drawMahal(nir_pc$T[,c(1,2)],
          center=apply(nir_pc$T[,c(1,2)],2,mean),
          covariance=cov(nir_pc$T[,c(1,2)]),
          quantile=0.975,xlab="PC1",ylab="PC2",
          col="blue")

drawMahal(nir_pc$T[,c(1,3)],
          center=apply(nir_pc$T[,c(1,3)],2,mean),
          covariance=cov(nir_pc$T[,c(1,3)]),
          quantile=0.975,xlab="PC1",ylab="PC3",
          col="blue")

 

drawMahal(nir_pc$T[,c(2,3)],
          center=apply(nir_pc$T[,c(2,3)],2,mean),
          covariance=cov(nir_pc$T[,c(2,3)]),
          quantile=0.975,xlab="PC2",ylab="PC3",
          col="blue")



As we can see with 40 samples we have a wide distribution, being one of them an outlier in the PC1, but we want to fill the space with more samples spending as less as possible in laboratory analysis and this will be the idea along this tutorial.
What is the reason for the outlier?. Checking the raw spectra and specially the math treated spectra we can realize that one of the samples moves apart from the rest, being the water band, at 1940 nm aprox. lower than the rest, this mean that this sample is more dry than the others and that is the reason the sample is an outlier in the first PC component.

 

Tidyverse & Chemometrics (part 2)

In the previous post we have seen the histograms for the 40 fish meal samples, now it is time to look to the spectra.

Spectra can be is a matrix called "nir", as far as we have 40 samples and the spectrum of each sample goes from 1100 to 2498 every 2 nm, we have a matrix with 40 rows and 700 columns.

To use "matplot" is a quick option lo look to the raw spectra:

wavelength<-seq(1100,2498,2)
matplot(wavelength,t(nir),type="l",col="blue",
        xlab="nm",ylab = "log1/R",
        main="Fish1 spectra")

And we get this plot:

To remove the scatter and to improve the bands resolution we combine SNV and Detrend with the first derivative with a gap of 1 and a segment of 4.

library(prospectr)
nir_snvdt<-detrend(nir,wavelength)
nir_1d<-gapDer(nir_snv,m=1,w=1,s=4)
colnames(nir_1d)
wavelength_1d<-seq(1108,2490,2)
matplot(wavelength_1d,t(nir_1d),type="l",col="blue",
        xlab="nm",ylab = "log1/R",
        main="Fish1 1D spectra")

We can see how the scatter is reduced and how the band`s resolution increase.

We can try with the second derivative, with a gap of 1 and a segment of 10.

nir_2d<-gapDer(nir_snv,m=2,w=1,s=10)
colnames(nir_2d)
wavelength_2d<-seq(1130,2466,2)
matplot(wavelength_2d,t(nir_2d),type="l",col="blue",
        xlab="nm",ylab = "log1/R",
        main="Fish1 2D spectra")

We are going to use these math treated spectra to develop the Principal Component Analysis and PLS regressions in the coming posts.

Tidyverse & Chemometrics (Part 1)

Along this year, I am going to write a number of posts to work with R and Near Infrared spectra using the "tidyverse" package (and also others), so this will be a way I will learn how to use it.

 

So let´s start with some fish meal data (40 samples) acquired during a certain fishing season , so this will be our starter database and we will updated and validated using spectra of new samples from  other seasons.

 

The constituents are Moisture, Protein, Fat and Ash, and they are in a data frame with other columns as the Sample ID, Month, Year and of course the absorbances at a certain range (from 1100 to 2448 nm).

 

The name of the data frame is "fish1_df", and the first thing I do is to get the histograms to check the distributions for the parameters:

library(tidyverse)
#for the Dry Matter:
fish1_df %>%
ggplot(aes(dm))+
  geom_histogram(fill="blue",col="black")+
  geom_density(alpha=.2, fill="#FF6666") +
  geom_vline(aes(xintercept=mean(dm, na.rm=T)),
               color="red", linetype="dashed", size=1)

 

fish1_df %>%
ggplot(aes(protein))+
  geom_histogram(fill="red",col="black")+
  geom_density(alpha=.2, fill="#FF6666") +
  geom_vline(aes(xintercept=mean(protein, na.rm=T)),
               color="blue", linetype="dashed", size=1)

 

ggplot(aes(fat))+
  geom_histogram(fill="green",col="black")+
  geom_density(alpha=.2, fill="#FF6666") +
  geom_vline(aes(xintercept=mean(fat, na.rm=T)),
               color="blue", linetype="dashed", size=1)

fish1_df %>%
ggplot(aes(ash))+
  geom_histogram(fill="brown",col="black")+
  geom_density(alpha=.2, fill="#FF6666") +
  geom_vline(aes(xintercept=mean(ash, na.rm=T)),
               color="red", linetype="dashed", size=1)

 

 

EQA "Estimated Minimum" and "Estimated Maximum" in Win ISI

One of the statistics we can see after an ".eqa" file is generated are these : "Estimated Minimum "  and the "Estimated Maximum".

 

These values are calculated from the ".cal" file (file with the laboratory values distribution) used to develop the calibration for a certain constituent.

 

Take into account that the  ".cal" used to develop the equation can be different than the original ".cal" due that some samples can be taken out as outliers (the ones specified in the ".oll" file).

 

So from the ".cal" file used for the equation file we calculate the Mean Value and the Standard Deviation Value and the Estimated Minimum is calculates as Mean - 3.SD and the "Estimated Maximum is calculated as Mean + 3.SD.

 

If we get a negative value for the Estimated Minimum we replaced it by cero.