PLSR Demonstration

This is an example of PLSR on simulated spectroscopic data to determine the concentration of red and blue dye in various solutions. The M-file may be downloaded here. The M-file implementing the NIPALS algorithm may be downloaded here (html version).

Contents

Constructing the data

Let's contruct spectra for 4 dyes (2 red, 1 blue, and 1 yellow) at 40 points (x, that presumabely correspond to wavelengths):

randn('state', 0);
numInputs = 40;
x = linspace(-8,12,numInputs);
red1Spectrum = 10*evpdf(-x,-2,2);
red2Spectrum = 9.8*evpdf(-x,-1.5,2);
blueSpectrum = 4*evpdf(-x,4,1.);
yellowSpectrum = 5*evpdf(-x,1,1.2);

Spectra plots

Combine the spectra into one array and plot them

dyeSpectra = [
    red1Spectrum
    red2Spectrum
    blueSpectrum
    yellowSpectrum];
set(gcf,'DefaultAxesColorOrder',[1 0 0; 0.5 0 0; 0 0 1; 1 1 0]);
plot(1:numInputs,dyeSpectra);
xlabel('Wavelength')
ylabel('Optical Density')
legend('red1','red2','blue','yellow')
title('Spectra for dyes')

Create stock solutions

Let's have 3 stock solutions:

stocks= [
    0.67 0.27 0.00 0.06
    0.36 0.56 0.00 0.08
    0.00 0.00 0.95 0.05];

Create sample solutions

Let's make samples of each stock at different concentrations, and a few mixtures:

samples1 = [0.01 0.02 0.04 0.08]'*stocks(1,:);
samples2 = [0.03 0.06 0.09 0.12]'*stocks(2,:);
samples3 = [0.01 0.02 0.04 0.08 0.16]'*stocks(3,:);
samples13 = ([0.03 0.05 0.08]'*stocks(1,:) + [0.07 0.05 0.02]'*stocks(3,:));
samples23 = ([0.02 0.03 0.04]'*stocks(2,:) + [0.03 0.02 0.01]'*stocks(3,:));
samples = [samples1; samples2; samples3; samples13; samples23];

Create sample spectra

sampleSpectra = samples*dyeSpectra;
sampleSpectra = max(1e-6, sampleSpectra + 0.002*randn(size(sampleSpectra)) );

plot(sampleSpectra(end-2,:),'k-')
xlabel('Wavelength')
ylabel('Optical Density')
title('Spectrum for a mixture of red + blue')

Construct the raw inputs (X) and outputs (Y).

The inputs (Xraw) are just the sample spectra.

The two outputs (Yraw) are

Xraw = sampleSpectra;
Yraw = [samples(:,1)+samples(:,2), samples(:,3)];

Write the outputs to a text file for analysis in Simca-P

csvwrite('spec.txt',[Yraw Xraw])

Apply PLSR to the data

First mean-center the data. However, we should not variance scale the data because all the are all in the same units, and larger data points really should be weighted more heavily (because they are more informative).

Then, let's extract three principal components

numSamples = size(Xraw,1);
numOutputs = size(Yraw,2);
Xmean = mean(Xraw,1);
X = Xraw - repmat(Xmean, [numSamples 1]);
% Xmean = mean(Xraw,2);
% X = Xraw - repmat(Xmean, [1 numInputs]);
Ymean = mean(Yraw,1);
Y = Yraw - repmat(Ymean, [numSamples 1]);
% Ymean = mean(Yraw,2);
% Y = Yraw - repmat(Ymean, [1 numOutputs]);
[P,Q,W,B] = nipals(X,Y,3);
T = X*W;
U = Y*Q;
Ypred = X*W*B*Q';
YrawPred = Ypred + repmat(Ymean, [numSamples 1]);
% YrawPred = Ypred + repmat(Ymean, [1 numOutputs]);

Analysis of predictive ability

plot(Yraw(:,1),Yraw(:,2),'kx',YrawPred(:,1),YrawPred(:,2),'bo')
xlabel('Red Conc.')
ylabel('Blue Conc.')
legend('Expt.', 'Pred.')
title('Comparison of Predictive Ability')

Analysis of the Y loadings

PC1 is has a positive "red" weight and a negative "blue" weight. PC2 is mostly "blue", but some "red."

clf
colormap([1 0 0; 0 0 1]);
bar(Q')
set(gca,'XTickLabel',{'PC1','PC2','PC3'})
legend('Red','Blue');
title('Y Loadings')

Analysis of the X loadings

The X loadings are linear combinations of "red" and "blue" spectra. As indicated by the Y loadings, PC1 is "red" minus "blue", and PC2 is "blue" plus a little "red". PC3 looks like a lot of noise. As we'll later see using Simca-P, PC3 is not statistically significant.

clf
plot(1:numInputs,P);
legend('PC1','PC2','PC3')
xlabel('Arbitrary Wavelength')
title('X Loadings');

Analysis of the scores

s = {samples1, samples2, samples3, samples13, samples23};
legendText = {'1','2','3','1+2','2+3'};
symbols = 'xo+s^';
colors = [1 0 0; 0.5 0 0; 0 0 1; 0.7 0 0.7; 0.3 0 0.3];
clf
hold on
idx1 = 1;
for i=1:numel(s)
    idx2 = idx1 + size(s{i},1) - 1;
    Ti = T(idx1:idx2,:);
    plot(Ti(:,1),Ti(:,2),symbols(i),'MarkerFaceColor',colors(i,:),'MarkerEdgeColor',colors(i,:));
    idx1 = idx2+1;
end
xlabel('PC1');
ylabel('PC2');
title('X Scores')
legend(legendText{:});
hold off

Training the PLSR model mostly with blue solutions

At first, I was surprised to find that each principal component doesn't correspond to a single color. Then I realized it has a lot to do with the data you use to train the model. The first principal component is the one that explains the most variance in your data. So if most of your samples are just "blue", then PC1 will be a blue spectrum. If most of your samples are mixtures, then PC1 will be a mixture of the spectra.

Let's confirm this by adding a bunch more blue solutions to our samples.

n = 100;
blueConc = 0.5*randn(n,1);
samplesB = zeros(n,4);
samplesB(:,3) = blueConc;
samplesB = [samples; samplesB];
YrawB = [samplesB(:,1)+samplesB(:,2), samplesB(:,3)];
sampleSpectraB = blueConc*blueSpectrum;
sampleSpectraB = sampleSpectraB + max(1e-6,0.002*randn(size(sampleSpectraB)));
sampleSpectraB = [sampleSpectra; sampleSpectraB];
XrawB = sampleSpectraB;

csvwrite('specB.txt',[YrawB XrawB])

numSamplesB = size(XrawB,1);
numOutputsB = size(YrawB,2);
XmeanB = mean(XrawB,1);
XB = XrawB - repmat(XmeanB, [numSamplesB 1]);
YmeanB = mean(YrawB,1);
YB = YrawB - repmat(YmeanB, [numSamplesB 1]);
[PB,QB,WB,BB] = nipals(XB,YB,3);
clf
colormap([1 0 0; 0 0 1]);
bar(QB')
set(gca,'XTickLabel',{'PC1','PC2','PC3'})
legend('Red','Blue');
title('Y Loadings')