PLSの計算
PLSを手計算し、その後関数での結果と比べてみる。
x1 <- c(0.966,0.207,-0.552,-0.931,0.587,1.55,-0.931,0.587,-1.69,0.207)
x2 <- c(0.619,-1.009,-0.358,-0.684,0.944,1.921,0.293,0.293,-1.335,-0.684)
x3 <- c(-0.474,-0.474,-1.684,0.735,-0.474,1.849,-0.474,0.735,-0.474,0.735)
y1 <- c(1.065,-0.821,0.061,-0.957,0.311,1.813,-0.452,0.529,-1.55,0.001)
y2 <- c(-0.767,-0.067,-0.764,0.195,-1.062,1.933,-0.055,0.443,-1.093,1.236)
標準化されているので平均0、分散1。
colMeans(cbind(x1,x2,x3,y1,y2))
x1 x2 x3 y1 y2 -1.665335e-17 0.000000e+00 1.110223e-17 -5.637851e-18 -1.000000e-04
cov(cbind(x1,x2,x3,y1,y2))
x1 x2 x3 y1 y2 x1 1.0005353 0.7521066 0.4558819 0.8955430 0.4747129 x2 0.7521066 1.0000464 0.3995542 0.8875398 0.3523491 x3 0.4558819 0.3995542 0.9998569 0.4023940 0.8534377 y1 0.8955430 0.8875398 0.4023940 0.9997969 0.4593081 y2 0.4747129 0.3523491 0.8534377 0.4593081 1.0000501
引っ付ける。
X <- cbind(x1,x2,x3)
x1 x2 x3 [1,] 0.966 0.619 -0.474 [2,] 0.207 -1.009 -0.474 [3,] -0.552 -0.358 -1.684 [4,] -0.931 -0.684 0.735 [5,] 0.587 0.944 -0.474 [6,] 1.550 1.921 1.849 [7,] -0.931 0.293 -0.474 [8,] 0.587 0.293 0.735 [9,] -1.690 -1.335 -0.474 [10,] 0.207 -0.684 0.735
Y <- cbind(y1,y2)
y1 y2 [1,] 1.065 -0.767 [2,] -0.821 -0.067 [3,] 0.061 -0.764 [4,] -0.957 0.195 [5,] 0.311 -1.062 [6,] 1.813 1.933 [7,] -0.452 -0.055 [8,] 0.529 0.443 [9,] -1.550 -1.093 [10,] 0.001 1.236
この10×3のXデータと、10×2のYデータに対し、NIPALSアルゴリズムでPLSモデルを作成してみる。
第二スコアまで取り出すとする。
まず、第一スコアを取り出す。
X'Yに特異値分解を適用。
svd_res <- svd(t(X)%*%Y)
$d [1] 14.364201 4.781158 $u [,1] [,2] [1,] -0.6282229 -0.2793742 [2,] -0.5787447 -0.4559675 [3,] -0.5199909 0.8450111 $v [,1] [,2] [1,] -0.8054403 -0.5926769 [2,] -0.5926769 0.8054403
第一Xスコア(t1)の重み(w1)。
w1 <- svd_res$u[,1]
[1] -0.6282229 -0.5787447 -0.5199909
第一Yスコア(u1)の重み(c1)。
c1 <- svd_res$v[,1]
[1] -0.8054403 -0.5926769
第一Xスコアの計算。
t1 <- X %*% w1
[,1] [1,] -0.7186306 [2,] 0.7003870 [3,] 1.4296343 [4,] 0.5985436 [5,] -0.6686262 [6,] -3.0469773 [7,] 0.6617790 [8,] -0.9205324 [9,] 2.0807966 [10,] -0.1163741
第一Yスコアの計算。
u1 <- Y %*% c1
[,1] [1,] -0.4032108 [2,] 0.7009759 [3,] 0.4036733 [4,] 0.6552344 [5,] 0.3789309 [6,] -2.6059077 [7,] 0.3966563 [8,] -0.6886338 [9,] 1.8962283 [10,] -0.7333540
第一Xスコアに対するX負荷量の計算。
p1_t <- solve( t(t1)%*%t1 ) %*% t(t1) %*% X
p1 <- t(p1_t)
[,1] x1 -0.6238024 x2 -0.6037275 x3 -0.4975260
第一Yスコアに対するY負荷量の計算。
q1_t <- solve( t(u1)%*%u1 ) %*% t(u1) %*% Y
q1 <- t(q1_t)
[,1] y1 -0.7490939 y2 -0.6692510
第一Xスコアまでで説明されない残差部分。
t1 p1' が、第一XスコアまででのXデータの近似。
Ex1 <- X - t1 %*% t(p1)
x1 x2 x3 [1,] 0.51771653 0.18514296 -0.8315374 [2,] 0.64390304 -0.58615715 -0.1255393 [3,] 0.33980926 0.50510951 -0.9727197 [4,] -0.55762710 -0.32264280 1.0327910 [5,] 0.16990943 0.54033202 -0.8066589 [6,] -0.35071160 0.08145614 0.3330496 [7,] -0.51818070 0.69253416 -0.1447477 [8,] 0.01276975 -0.26275067 0.2770112 [9,] -0.39199419 -0.07876595 0.5612504 [10,] 0.13440558 -0.75425822 0.6771009
第一Yスコアまでで説明されない残差部分を計算する。
第一YスコアまででのYデータの近似だが、u1 q1'ではない。ここで、Xスコアを説明変数、Yスコアを目的変数とした単回帰を行う。
D行列の第一対角成分(d1)を計算。(切片なしの回帰係数)
d1 <- solve(t(t1)%*%t1) %*% t(t1) %*% u1
[,1] [1,] 0.7653233
Ey1 <- Y - as.numeric(d1) * t1 %*% t(q1)
y1 y2 [1,] 0.653009762 -1.13507785 [2,] -0.419468828 0.29173357 [3,] 0.880607987 -0.03175074 [4,] -0.613855536 0.50157007 [5,] -0.072322735 -1.40446590
以上、第一スコア関連の計算。
第二スコアの取り出しに進む。
第一スコアで行った計算の繰り返しになる。ただし、最後に計算した残差部分に対して計算を行う。
(X残差行列)'(Y残差行列) への特異値分解。
svd_res2 <- svd(t(Ex1)%*%Ey1)
$d [1] 4.7994939 0.2475326 $u [,1] [,2] [1,] -0.2760840 0.7274020 [2,] -0.4590238 -0.6740562 [3,] 0.8444376 -0.1285871 $v [,1] [,2] [1,] -0.5199760 0.8541809 [2,] 0.8541809 0.5199760
第二Xスコア(t2)の重み(w2)。
w2 <- svd_res2$u[,1]
[1] -0.2760840 -0.4590238 0.8444376
第二Yスコア(u2)の重み(c2)。
c2 <- svd_res2$v[,1]
[1] -0.5199760 0.8541809
第二Xスコアの計算。
t2 <- Ex1 %*% w2
[,1] [1,] -0.93009969 [2,] -0.01472129 [3,] -1.14707426 [4,] 1.17418014 [5,] -0.97610760 [6,] 0.34067511 [7,] -0.29705871 [8,] 0.35100196 [9,] 0.61831969 [10,] 0.88088466
第二Yスコアの計算。
u2 <- Ey1 %*% c2
[,1] [1,] -1.30911119 [2,] 0.46730695 [3,] -0.48501586 [4,] 0.74762169 [5,] -1.16206186 [6,] 0.28365415 [7,] 0.28030389 [8,] -0.02499112 [9,] 0.16241458 [10,] 1.03902459
第二Xスコアに対するX負荷量の計算。
p2_t <- solve( t(t2)%*%t2 ) %*% t(t2) %*% Ex1
p2 <- t(p2_t)
[,1] x1 -0.2978230 x2 -0.4388790 x3 0.8482805
第二Yスコアに対するY負荷量の計算。
q2_t <- solve( t(u2)%*%u2 ) %*% t(u2) %*% Ey1
q2 <- t(q2_t)
[,1] y1 -0.3707696 y2 0.9450091
第二Xスコアまでで説明されない残差部分。
Ex2 <- Ex1 - t2 %*% t(p2)
x1 x2 x3 [1,] 0.240711453 -0.223058287 -0.0425519806 [2,] 0.639518696 -0.592618023 -0.1130514849 [3,] -0.001815838 0.001682669 0.0003209963 [4,] -0.207929247 0.192680246 0.0367568772 [5,] -0.120797870 0.111938862 0.0213541507 [6,] -0.249250712 0.230971299 0.0440615157 [7,] -0.606651611 0.562161326 0.1072413768 [8,] 0.117306205 -0.108703266 -0.0207369084 [9,] -0.207844368 0.192601591 0.0367418726 [10,] 0.396753293 -0.367656417 -0.0701364154
D行列の第二対角成分(d2)を計算。
d2 <- solve(t(t2)%*%t2) %*% t(t2) %*% u2
[,1] [1,] 0.8001358
第二Yスコアまでで説明されない残差部分。
Ey2 <- Ey1 - as.numeric(d2) * t2 %*% t(q2)
y1 y2 [1,] 0.3770808 -0.4317964 [2,] -0.4238361 0.3028649 [3,] 0.5403100 0.8355929 [4,] -0.2655162 -0.3862693 [5,] -0.3619007 -0.6663962 [6,] 0.1672376 0.1147618 [7,] -0.1607300 0.5085755 [8,] 0.1053900 -0.2938957 [9,] -0.1736469 -0.4947633 [10,] 0.1956114 0.5103257
以上、第二スコア関連の計算。第一、第二スコアの取出し完了。
第二スコアまでの結果をひっつけ行列にする。
Xスコア(T)。
T <- cbind(t1, t2)
[,1] [,2] [1,] -0.7186306 -0.93009969 [2,] 0.7003870 -0.01472129 [3,] 1.4296343 -1.14707426 [4,] 0.5985436 1.17418014 [5,] -0.6686262 -0.97610760 [6,] -3.0469773 0.34067511 [7,] 0.6617790 -0.29705871 [8,] -0.9205324 0.35100196 [9,] 2.0807966 0.61831969 [10,] -0.1163741 0.88088466
Yスコア(U)。
U <- cbind(u1, u2)
[,1] [,2] [1,] -0.4032108 -1.30911119 [2,] 0.7009759 0.46730695 [3,] 0.4036733 -0.48501586 [4,] 0.6552344 0.74762169 [5,] 0.3789309 -1.16206186 [6,] -2.6059077 0.28365415 [7,] 0.3966563 0.28030389 [8,] -0.6886338 -0.02499112 [9,] 1.8962283 0.16241458 [10,] -0.7333540 1.03902459
X負荷量(P)。
P <- cbind(p1, p2)
[,1] [,2] x1 -0.6238024 -0.2978230 x2 -0.6037275 -0.4388790 x3 -0.4975260 0.8482805
Y負荷量(Q)。
Q <- cbind(q1, q2)
[,1] [,2] y1 -0.7490939 -0.3707696 y2 -0.6692510 0.9450091
D行列(D)。
D <- diag(c(d1,d2))
[,1] [,2] [1,] 0.7653233 0.0000000 [2,] 0.0000000 0.8001358
X重み(W)。
W <- cbind(w1, w2)
w1 w2 [1,] -0.6282229 -0.2760840 [2,] -0.5787447 -0.4590238 [3,] -0.5199909 0.8444376
Y重み(C)。
C <- cbind(c1, c2)
c1 c2 [1,] -0.8054403 -0.5199760 [2,] -0.5926769 0.8541809
NIPALSアルゴリズムでは、最終的な回帰係数(B)は、B=W (P'W)^-1 C' と教科書では表記されていると思う。ただし、このCはY重みのことではないようである。
Y重み(C)ではないので、C*とし、B=W (P'W)^-1 C*' と書き直す。
C*'は、Xスコア(T)でYを回帰したときの係数である。
Y=T C*'+Error
よって、C*'= (T' T)^-1 T' Y
B <- W %*% solve(t(P) %*% W) %*% solve( t(T) %*% T ) %*% t(T) %*% Y
y1 y2 [1,] 0.49447906 0.1088077 [2,] 0.54069342 -0.0396553 [3,] -0.03711641 0.8233860
よって、Yの予測値はXBで計算される。
X %*% B
y1 y2 [1,] 0.8299492 -0.3097233 [2,] -0.4256093 -0.3277496 [3,] -0.4040167 -1.4324472 [4,] -0.8574749 0.5310129 [5,] 0.8182670 -0.3638494 [6,] 1.7364864 1.6149147 [7,] -0.2843437 -0.5032039 [8,] 0.4214018 0.6574398 [9,] -1.5399022 -0.5212301 [10,] -0.2947577 0.6548361
Rの関数を使用する。
chemometricsライブラリのpls2_nipals()関数を利用。この関数はNIPALSアルゴリズム。
plsライブラリのmvr()関数を使用すると、SIMPLSやKernel法を使えるとのこと。ちなみにNIPALSとKernel法は一致するようである。
library(chemometrics)
第二スコアまで取り出す。引数aはスコア数。もともと標準化しているデータなのでscale(標準化)は行わない。
pls2_res <- pls2_nipals(X, Y, a=2, scale = FALSE)
$P [,1] [,2] x1 0.6238024 0.2978230 x2 0.6037275 0.4388790 x3 0.4975260 -0.8482805 $T [,1] [,2] [1,] 0.7186306 0.93009969 [2,] -0.7003870 0.01472129 [3,] -1.4296343 1.14707426 [4,] -0.5985436 -1.17418014 [5,] 0.6686262 0.97610760 [6,] 3.0469773 -0.34067511 [7,] -0.6617790 0.29705871 [8,] 0.9205324 -0.35100196 [9,] -2.0807966 -0.61831969 [10,] 0.1163741 -0.88088466 $Q [,1] [,2] y1 0.7490939 0.3665540 y2 0.6692510 -0.9475753 $U [,1] [,2] [1,] 0.4032701 1.25693839 [2,] -0.7009166 -0.41662732 [3,] -0.4036140 0.58855239 [4,] -0.6551752 -0.70432381 [5,] -0.3788716 1.11351346 [6,] 2.6059670 -0.50458887 [7,] -0.3965970 -0.23242261 [8,] 0.6886931 -0.04181580 [9,] -1.8961691 -0.01168087 [10,] 0.7334133 -1.04754497 $D [1] 0.7653233 0.8001358 $W [,1] [,2] x1 0.6282229 0.2760840 x2 0.5787447 0.4590238 x3 0.5199909 -0.8444376 $C [,1] [,2] y1 0.8054403 0.5199760 y2 0.5926769 -0.8541809 $B y1 y2 x1 0.64000890 0.15218608 x2 0.69603317 -0.03463704 x3 -0.02816547 1.04246642
X負荷量。
P_r <- pls2_res$P
[,1] [,2] x1 0.6238024 0.2978230 x2 0.6037275 0.4388790 x3 0.4975260 -0.8482805
Xスコア。
T_r <- pls2_res$T
[,1] [,2] [1,] 0.7186306 0.93009969 [2,] -0.7003870 0.01472129 [3,] -1.4296343 1.14707426 [4,] -0.5985436 -1.17418014 [5,] 0.6686262 0.97610760 [6,] 3.0469773 -0.34067511 [7,] -0.6617790 0.29705871 [8,] 0.9205324 -0.35100196 [9,] -2.0807966 -0.61831969 [10,] 0.1163741 -0.88088466
Y負荷量。
Q_r <- pls2_res$Q
[,1] [,2] y1 0.7490939 0.3665540 y2 0.6692510 -0.9475753
Yスコア。
U_r <- pls2_res$U
[,1] [,2] [1,] 0.4032701 1.25693839 [2,] -0.7009166 -0.41662732 [3,] -0.4036140 0.58855239 [4,] -0.6551752 -0.70432381 [5,] -0.3788716 1.11351346 [6,] 2.6059670 -0.50458887 [7,] -0.3965970 -0.23242261 [8,] 0.6886931 -0.04181580 [9,] -1.8961691 -0.01168087 [10,] 0.7334133 -1.04754497
これの第二Yスコアが手計算からずれるが、理由はわからない。
D行列。
D_r <- diag(pls2_res$D)
[,1] [,2] [1,] 0.7653233 0.0000000 [2,] 0.0000000 0.8001358
X重み。
W_r <- pls2_res$W
[,1] [,2] x1 0.6282229 0.2760840 x2 0.5787447 0.4590238 x3 0.5199909 -0.8444376
Y重み。
C_r <- pls2_res$C
[,1] [,2] y1 0.8054403 0.5199760 y2 0.5926769 -0.8541809
最終的な回帰係数。
B_r <- pls2_res$B
y1 y2 x1 0.64000890 0.15218608 x2 0.69603317 -0.03463704 x3 -0.02816547 1.04246642
この最終的な回帰係数であるが、上で手計算したものと異なる。この結果は、B=W (P'W)^-1 C' で、CをC*とせず計算したもののようである。なぜこうなっているかは不明。
当然、Yの予測値を X %*% B_r とすることはできない。
念のため、plsパッケージのmvr()関数でPLSでの予測値を出しておく。
library(pls)
NIPALSと結果が一致するとされるkernel法を使用
a1 <- mvr(Y~X, ncomp=2, method="kernelpls")
予測値。
predict(a1)
y1 y2 1 0.8299492 -0.3098233 2 -0.4256093 -0.3278496 3 -0.4040167 -1.4325472 4 -0.8574749 0.5309129 5 0.8182670 -0.3639494 6 1.7364864 1.6148147 7 -0.2843437 -0.5033039 8 0.4214018 0.6573398 9 -1.5399022 -0.5213301 10 -0.2947577 0.6547361
これは手計算に一致する。
PLSはNIPALS、SIMPLSなどのアルゴリズムによっても数値がずれるし、ソフトウェアによっても途中での基準化や制約のかけ方が異なったりで数値がずれると思う。