Introduction

In this vignette, we consider a supervised dimensional reduction method partial least squares (PLS).

Test data is available from toyModel.

library("guidedPLS")
data <- guidedPLS::toyModel("Easy")
str(data, 2)

## List of 8
##  $ X1      : int [1:100, 1:300] 86 101 95 106 113 85 88 103 106 84 ...
##  $ X2      : int [1:200, 1:150] 106 81 91 101 91 105 111 81 113 105 ...
##  $ Y1      : int [1:100, 1:50] 101 77 77 87 101 89 111 113 101 112 ...
##  $ Y1_dummy: num [1:100, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
##  $ Y2      : int [1:200, 1:50] 107 81 102 90 84 106 97 90 88 115 ...
##  $ Y2_dummy: num [1:200, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
##  $ col1    : chr [1:100] "#66C2A5" "#66C2A5" "#66C2A5" "#66C2A5" ...
##  $ col2    : chr [1:200] "#66C2A5" "#66C2A5" "#66C2A5" "#66C2A5" ...

You will see that there are three blocks in the data matrix as follows.

suppressMessages(library("fields"))
layout(c(1,2,3))
image.plot(data$Y1_dummy, main="Y1 (Dummy)", legend.mar=8)
image.plot(data$Y1, main="Y1", legend.mar=8)
image.plot(data$X1, main="X1", legend.mar=8)

Partial Least Squares by Singular Value Decomposition (PLS-SVD)

Here, suppose that we have two data matrices \(X\) (\(N \times M\)) and \(Y\) (\(N \times L\)), and their column vectors are assumed to be centered. Considering the projection of them to lower dimension, we have the scores \(XV\) and \(YW\), where \(V\) (\(M \times K\)) and \(W\) (\(L \times K\)) are called loading matrices for \(X\) and \(Y\), respectively. Depending on the loading matrice, these scores can take various values, but consider \(V\) and \(W\) such that the covariance is maximized as follows:

\[ \max_{V,W} \mathrm{tr} \left( V^{T}X^{T}YW \right)\ \mathrm{s.t.}\ V^{T}V = W^{T}W = I_{K} \]

This is known as Partial Least Squares (PLS). The widely used PLS algorithms are NIPALS and SIMPLS but here we introduce more simple one, Partial Least Squares by Singular Value Decomposition (PLS-SVD (Lê Cao 2008)). PLS-SVD is solved by svd against cross-product matrix \(X'Y\). This formulation is also known as Guided Principal Component Analysis (guided-PCA (Reese S E 2013)).

out1 <- PLSSVD(X=data$X1, Y=data$Y1, k=2)

plot(out1$scoreX, col=data$col1, pch=16)

out2 <- PLSSVD(X=data$X1, Y=data$Y1, k=2, deflation=TRUE)

plot(out2$scoreX, col=data$col1, pch=16)

Partial Least Squares Discriminant Analysis (PLS-DA)

In many cases, \(X\) and \(Y\) in PLS are assumed to be continuous. On the other hand, some matrices containing discrete values may be set for \(Y\). For example, a dummy variable matrix containing only the values 0 and 1, which means the correspondence of data and group in \(X\), is used as \(Y\). Such a formulation is called Partial Least Squares Discriminant Analysis (PLS-DA (Lê Cao 2008)). In the context of distinguishing it from PLS-DA, the former is called Partial Least Squares Regression (PLS-R).

out3 <- PLSSVD(X=data$X1, Y=data$Y1_dummy, k=2, fullrank=TRUE)

plot(out3$scoreX, col=data$col1, pch=16)

Sparse Partial Least Squares Discriminant Analysis (sPLS-DA)

Sparse Partial Least Squares Discriminant Analysis (sPLS-DA (Lê Cao 2008)) is an extension of PLS-DA to enhance the interpretability of the scores and loadings by making the values sparse. Soft-thresholding function forces the values within a certain range (\(-\lambda\) to \(\lambda\)) to be adjusted to zero.

out4 <- sPLSDA(X=data$X1, Y=data$Y1_dummy, k=2, lambda=30)

plot(out4$scoreX, col=data$col1, pch=16)

layout(cbind(1:2, 3:4))
barplot(out3$loadingX[,1], main="PLSDA (1st loading)")
barplot(out3$loadingX[,2], main="PLSDA (2nd loading)")
barplot(out4$loadingX[,1], main="sPLSDA (1st loading)")
barplot(out4$loadingX[,2], main="sPLSDA (2nd loading)")

Comparison with unsupervised learning

Next, we use another data data2.

data2 <- guidedPLS::toyModel("Hard")
str(data2, 2)

## List of 8
##  $ X1      : int [1:100, 1:300] 81 110 109 117 98 108 104 89 87 103 ...
##  $ X2      : int [1:200, 1:150] 103 118 97 88 116 99 87 97 79 108 ...
##  $ Y1      : int [1:100, 1:50] 89 119 104 96 92 100 98 99 107 87 ...
##  $ Y1_dummy: num [1:100, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
##  $ Y2      : int [1:200, 1:50] 111 109 96 97 108 104 121 112 102 82 ...
##  $ Y2_dummy: num [1:200, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
##  $ col1    : chr [1:100] "#66C2A5" "#66C2A5" "#66C2A5" "#66C2A5" ...
##  $ col2    : chr [1:200] "#66C2A5" "#66C2A5" "#66C2A5" "#66C2A5" ...

data2 has some blocks at the same coordinates as data in the matrix, but the background noise is larger and the block values are set smaller than data, resulting in a low S/N ratio like below.

layout(t(1:2))
image.plot(data$X1, main="X1", legend.mar=8)
image.plot(data2$X1, main="X1 (Noisy)", legend.mar=8)

For such data, supervised learning has difficulty extracting the signal from the noise. For example, let’s perform an unsupervised learning method principal component analysis (PCA) using only \(X\).

out.pca <- prcomp(data2$X1, center=TRUE, scale=FALSE)

Next, perform PLS-DA, which is a supervised learning using \(X\), and label information \(Y\) as well.

out5 <- PLSSVD(X=data2$X1, Y=data2$Y1_dummy, k=2, fullrank=TRUE)

Comparing the two, it is clear that PLS-DA, which uses labels, captures the structure of the data better like below.

layout(t(1:2))
plot(out.pca$x[, 1:2], col=data2$col1, main="PCA", pch=16)
plot(out5$scoreX, col=data2$col1, main="PLS-DA", pch=16)

Let’s compare the two with the PLS-DA loading performed on data. In PCA loading, PC1 correlates well with PLS-DA loading, but the correlation is a bit smaller for PC2. In contrast, in PLS-DA (Noisy), the first and second loadings correlate well with the PLS-DA loadings. In other words, the use of labels mitigates the effect of noise.

layout(t(1:2))
image.plot(cor(out5$scoreX, out.pca$x[,1:2]),
xlab="PLS-DA (Noisy)", ylab="PCA (Noisy)", legend.mar=8)
image.plot(cor(out5$scoreX, out3$scoreX),
xlab="PLS-DA (Noisy)", ylab="PLS-DA", legend.mar=8)