One of the advantages of R, is that you can develop a function to apply to your calculations to get your purpose. In this case I wanted to apply the RMS calculation to every row of the difference matrix we have seen in the previous post.
We can see with an histogram the distribution of the RMS values
Once we have sorted and grouped all the samples and subsamples, we can create two matrices one that contains the odd samples (first subsample) and other with the even subsamples (second subsample), this way we can calculate the spectra difference between subsamples and check the spectra.
This difference
spectra matrix is important to calculate the RMS between subsamples, that gave
us an idea how similar are the subsamples between them. Just sum the squared values
for every spectrum difference, divide by the number of wavelengths and
calculate the square root.
Now we have four calibrations for clay in soil developed with a selection of samples from the LUCAS database (Spanish crop soils). Now we want to see if those calibrations predict with certain accuracy a new set of soil samples (155) from a Spanish region, acquired in a different instrument and at a different laboratory.
In these cases, it is normal to expect a bias or a slope
in the predictions, so we can use the model with those adjustments applied, until
the database is updated and a new expanded method with new variability
(instrument, laboratory, region,) developed.
Well, these are the results of the XY plots "Lab vs
Predictions" for this independent data set:
PLS predictions:
As we see in all the cases some adjustment or
calibration update is needed. We can try to reprocess everything trying to find
a better configuration which improves these values, but in that case this new
set will never be independent again like it is now.
Let´s select a
seed (to fix the training and test set) and develop the regression with four
algorithms (PLS, Random Forest, Cubist and Memory Based Learning.
These are the
Test Validation XY plots and statistics:
PLS Regression:
Cubist Regression:
To have a good
traceability in our data bases is important to develop accurate calibrations.
We have seen with the Soil LUCAS database that we can filter the data by sample
origin (Spain for my case) and after that filter it by land type (I choose for
this example “Croplands). The samples are split into training and test set
randomly.
Can we improve
the results with another type of regression for these cropland samples? This is
what we will see in the next coming posts.
Another
question can be: May I use this database to predict samples from another database?
These can be samples from a different area in the same country, taken with a
different instrument and analysed in a different laboratory. Check all this is
important to see the robustness of the calibrations.
As we saw in the post "Try to find high content gypsum samples (part 3)" , when we develop the principal component analysis as much variance as possible (we determine the explained variance limit) is explained, and unless we saw all the sores maps, we may have not a good idea about what is happening to our data.
The average spectrum represents all the groups, so we cannot expect that with the Mahalanobis distance all the samples with gypsum will be marked as outliers, because the average spectrum represents as well those samples. Only the ones with high and very high gypsum content can probably be marked as outliers by Mahalanobis distance. Now that we have make two samples sets ("gypsum content" and "non or low gypsum content"), we can check the sample sets separately to understand better our data.
Let´s check again the Mahalanobis distance plot once we have, by spectra visualization and discrimination by correlation, the two sample sets ("No" and "Yes" gypsum content):
If we go back to the post: "PCA with the first derivative" we saw how there was some groupings (in the PC2 vs. the rest PC maps) in the scores, after, when looking to second loading we though that could be the samples with gypsum content.
In the last two posts : "Trying to find high content gypsum samples (part 1)" and "Trying to find high content gypsum samples (part 2)", we found those samples by correlation with a Gypsum reference spectra, and now we can create a new variable called "gypsum" which takes the factor value "Yes" or "No", so now we can see easier in the score maps if that grouping was due to the content of gypsum:
pairs(scores_1df[ ,1:6], col = scores_1df$gypsum)
How can I find the high content gypsum samples in the LUCAS spanish database?. We have the spectrum of pure gypsum, so we can try to correlate (or meassure a distance) every spectrum of the database with the pure gypsum spectra and see if we find some threshold or gap that splits these high content gypsum samples.
I try in this case the correlation, because it is the easiest algorithm, and I don´t have a population of gypsum soil spectra (just one reference pure spectrum).
Let´s run the correlation:
corSG1 <- as.numeric(cor(t(lucas_spain$spcnir_SG[,21:546]),
In a previous post we had calculated the PCA with the
math treatment SNV+Detrend and we calculated a first sample set of outliers
with the Mahalanobis distance.
When calculating PCA we have to treat as best the
spectra as possible in order to detect populations or boundaries and if we
treat the spectra with math treatments which help to do this task is great. So indeed,
to apply the SNV + Detrend, I apply this time the first derivative to the
spectra (soil from Spanish soil from LUCAS database) and calculate the PCA.
We can have a look to the score’s
maps (six PC recommended) to find if there are boundaries on them:
Look at the maps which include
PC2 as one of the axes. There are a certain number of samples which takes a
different direction than the rest. This second PC term can be useful to find
something interesting on the spectra.
To see what is happening we can see the loading spectra for this second PC term:
This loading spectra has the
first derivative math treatment applied, so we can compare it with a library of
reference known spectra (minerals in this case) to see which is the best match,
and in this case the best match is with the gypsum mineral, so this second
term is explaining part of the variance included in our spectra database due to
the addition of gypsum to the soil.
In order to find if the soil has traces of a certain mineral, it is useful to overplot our soil samples (in this case the soil Spanish samples included in the LUCAS database, with the pure mineral spectrum (in this case Gypsum). We must overplot them with the same math treatment and in the same scale. I do it, in this case, with the spectra treated with the SG second derivative.
Second derivative is quite helpfull to find these matches. If we do it, for example with SNV+Detrend, some of the bands are hiden by other samples and the assumption that there are samples with gypsum is less clear.
In this first plot we see the spectra of the LUCAS spanish database treated with the SNV and Detrend math treatment of the "Prospectr" package, where we remove the quadratic trend (as we saw in the last post):
In the last post, I talk about the use of the function "detrend" in the package "prospectr" to remove the quadratic trend lines in the soil NIR spectra. We can see these trendlines overplotted to the spectra (I do it just for five of the soil spectra from the Spanish soil LUCAS database).
In red is the spectra treated with SNV and in blue the quadratic trend lines to apply to the SNV treated spectra, to remove them and convert them in a SNV + Detrend spectra.
The spectra from LUCAS database (ESDAC*) seem to come from a XDS (VIS + NIR) instrument, and this instrument give two options when exporting the spectra (every 0.5 nm, and every 2 nm), the data comes in 0.5 nm, and that means that we have 4200 data poits per spectrum, and that means a huge spectra matrix.
We can trim the spectra keeping just the spectral data every two nanometers, so we will have the reflectance values from 400 to 2498 nm every 2nm, so we have a less heavy matrix of 1050 data points per spectrum. This fuction will give us the spectra like a NIR6500.
Just create a sequence to select one column of every four, and call that function (like in my case) "trim05to2":
dim(lucas_spain$spc)This way we have
the samples IDs, but we don´t need the column name (wavelengths), because Win
ISI will create them when importing the spectra.
Once we have the txt file, we convert it to a ".nir" file which is
the spectra format for Win ISI with the Win ISI tool "Convert",
selecting the configuration from TXT to Win ISI.
After
filling the options Win ISI ask, we wait for the conversion succeeded
message and we have the ".nir" file ready to view:
