東京に棲む日々

データ分析、統計、ITを勉強中。未だ世に出ず。

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などのアルゴリズムによっても数値がずれるし、ソフトウェアによっても途中での基準化や制約のかけ方が異なったりで数値がずれると思う。