PLS部分最小二乗法

Partial Least Squares(PLS)

部分最小二乗法,偏最小二乗法

説明変数の相関がきついんですけど…

説明変数に回帰させる一般的な方法がOLSなのだが,に相関があると回帰がうまくいかない.

そこで,を主成分分析して,相関する変数をまとめた主成分を抽出する主成分回帰分析(PCR)というのがある.しかし,このの主成分分析にはの情報を含まないので,必ずしも主成分分析の結果が意味のあるものでないかもしれない.

そこでそこで,の情報も含めて主成分分析的な回帰をやろうというのがPLS.

もともと心理学の分野で開発されたらしいが,最近は化学の分野で多用されているらしい.特定の反応を説明するのに,サンプル数が数十個しかないのに対して,説明変数の候補が100個以上あったりする場合に,説明変数の候補を探し出すのに使われているようだ.

我々が行う購買行動調査のようなアンケートも,性別や年齢などの属性で説明できるような多純なものではなく,様々な消費性向のようなことを重ねて聞くから,どれとどれがどれと関連しているかを突き止めるのが大変.

PLSはあくまで連続変数間の関係なので,アンケート分析にどれだけ使えるかわからないが,なんかの参考にはなるだろう.PLSは被説明変数に複数の変数をとれるというのも魅力的だ.

気になったので勉強しようとしてみた.ところが,どうも教科書が無いのか?ネット上の解説をあさってみた.上記のような話が出てきたが,具体的に何をしているかがよくわからない.たくさん抽出できる主成分の中からどうやって主成分を抽出するかという手順が説明してある.しかし,その前にそもそもどんな回帰をしているのか?その説明が飛んでいる.PLSの場合,この手順が醍醐味なのはよくわかったが,言ってしまえば,そこは計算機の仕事でしょ,ということになる.それよりも,結果として出力される値の意味が知りたい.少なくとも分かった気になりたい.

それで,自分でデータを作成してやってみた.ネットの説明は,Rのパッケージに添付されているデータを使っているので,出力された値の理解がむずかしい.推定モデルというのは,まずは真の関係を想定する.そこにいろんなノイズが入り込む.なかなか真の関係にはたどり着かないのだけれど,できるだけそれに近い,あるいはもっともらしい解を見つけようとする.それはそれでよいのだが,しかし,例題で使われているのは,現実のデータであったりするから,真の関係が果たしてどうなのかというのは本当のところは分からない.だから推定された値が妥当なのか,的外れたものなのかがついぞわからん.

そこで,ひとまずノイズを含まない真の関係にあるデータを人為的に作り出して,それを分析にかけてみて,分析手法で何を明らかにできるのか,推定される値が何を意味しているのかをを確かめてみよう.

それからちょっとずつノイズを加えていき,その推定方法の信頼性や強みを見ていったらよいのではないか.

以下,最初に,ノイズを含まない被説明変数が2つ,説明変数が2つで,それぞれきつい相関があるデータを作成し,それを分析して,PLSの意味を理解していく.

次に,2×2変数はそのまま,攪乱項を加えた場合に,PLSがどう作動するかを確認する.

そして,さらに余計な説明変数をたくさん加えて,そうした雑音の中から,真の関係をどの程度抽出できるものなのかについて検証してみる.

library(tidyverse)
library(MASS)
## Warning: package 'MASS' was built under R version 4.0.4
library(pls)
library(scatterplot3d)

2×2変数でノイズなし

まず,が2変数,が2変数の対応を確認する

データの作成

の作成

PLSでは全ての変数を標準化(平均0,分散1に)する.
平均0,分散1,共分散0.8の強めの相関がある51個の2変数x1を作成する.

x1_0<-c( 0, 0,
         1, 1,
        -1,-1,
         1,-1,
         -1,1)%>%
  matrix(.,ncol=2,byrow = T)
set.seed(123)
x1_1<-mvrnorm(46,
              mu=c(0,0),
              Sigma=matrix(c(46/45,0.8888889,
                             0.8888889,46/45),
                           nrow=2),
              empirical=T)
x1<-rbind(x1_0,x1_1)

平均は2変数とも0

apply(x1_1,2,mean)%>%
  round(4)
## [1] 0 0

分散・共分散は

var(x1)
##      [,1] [,2]
## [1,]  1.0  0.8
## [2,]  0.8  1.0

・・・で,のデータは,うまいことできた!

プロット図に書いてみよう

データの変換をプロット図で追いかけられるように,2~5行目のデータは,特殊な点を無理やり入れてある.その点をつなぐ線を描きたい.

その関数を定義

f_grid1<-function(x1){
  segments(x1[2,1],x1[2,2],x1[4,1],x1[4,2],col=3)
  segments(x1[2,1],x1[2,2],x1[5,1],x1[5,2],col=3)
  segments(x1[3,1],x1[3,2],x1[4,1],x1[4,2],col=2)
  segments(x1[3,1],x1[3,2],x1[5,1],x1[5,2],col=2)
}

のプロット図.1が原点,2が(1,1),3が(-1,-1),4が(1,-1),5が(-1,1).

plot(x1,asp=1)
grid()
text(x=x1[1:5,1],y=x1[1:5,2],1:5,pos = 4)
f_grid1(x1[1:5,])

の主成分のスコアの計算

の主成分分析をすると・・・

prc_x1<-prcomp(x1)

固有ベクトルは,第1主成分が右上がりの45度線,第2主成分がその直交(右下がり45度線)

prc_x1$rotation
##            PC1        PC2
## [1,] 0.7071068  0.7071068
## [2,] 0.7071068 -0.7071068

固有値は

prc_x1$sdev^2%>%
  round(4)
## [1] 1.8 0.2

スコアを求める.

Px1<-prc_x1$x

スコアのプロット図:横軸が第1主成分,縦軸が第2主成分.
のプロットを45度度マイナスにした形となる.

plot(Px1,asp=1)
grid()
text(x=Px1[1:5,1],y=Px1[1:5,2],1:5,pos=4)
f_grid1(Px1)

のプロット図の45度線でスコアを測っていることがわかる.
の2変数とも標準化してしまうと,どうあがいても,固有ベクトルの方向は,45度線になってしまう.

第1主成分がの方のプロット図の右上がり45度線.

第2主成分が同図の右下がり45度線

つまり,固有ベクトルは,45度()線で,長さ(ノルム)1なので・・・

ということ.

各主成分で符号が逆になる場合があるが,どっち向きにスコアを測るかが違うだけで,本質的に同じ.

なので,と置いておくと・・・

a<-1/sqrt(2)

スコアは,(第1主成分),(第2主成分)

cbind(a*x1[,1]+a*x1[,2],a*x1[,1]-a*x1[,2])%>%
  round(4)%>%
  head()
##         [,1]    [,2]
## [1,]  0.0000  0.0000
## [2,]  1.4142  0.0000
## [3,] -1.4142  0.0000
## [4,]  0.0000  1.4142
## [5,]  0.0000 -1.4142
## [6,] -1.0460  0.1174

ちゃんと主成分分析のスコアと一致しますね.

Px1%>%
  round(4)%>%
  head()
##          PC1     PC2
## [1,]  0.0000  0.0000
## [2,]  1.4142  0.0000
## [3,] -1.4142  0.0000
## [4,]  0.0000  1.4142
## [5,]  0.0000 -1.4142
## [6,] -1.0460  0.1174

固有値は,それぞれスコアの分散なんですね,これが!

さらに,なので,
第1主成分の固有値は,この場合,
第2主成分の固有値は,この場合,

var(prc_x1$x)%>%
  round(4)
##     PC1 PC2
## PC1 1.8 0.0
## PC2 0.0 0.2

を作成する

方針

のダイレクトな規定関係ではなく,複数ののスコアと複数ののスコアとの関係があるものと考えよう.

観測できる変数(顕在変数)は,観測できない変数(潜在変数,のスコアとのスコア)がノイズをもって観測されたものと考える.

最も簡単な例として,が上記2変数で,も2変数.
が1変数という場合もあるが,PLSが2変数以上を被説明変数にできるというのだから,そこは2変数としよう.)

は,それぞれ2変数とも平均0,標準偏差1に標準化.

ひとまず,にノイズは一切入らないとする.

のスコアをのスコアをと表し・・・

・・・でのスコアが決まるとする.

は,スコアから45度逆回転させればよいのだから.

で, を求めることができる.

手順

まずの2つの主成分のスコアを求める.

第1主成分X$の主成分のスコアと次のような関係にあるとすると・・・(0.7,
1.3という値は適当に決めた)

Py1_1<-0.7*Px1[,1]+1.3*Px1[,2]

この時,の分散は

var(Py1_1)
## [1] 1.22

の第2主成分

は,これと直交するから

となるようにを決めたい.

さらに,2つのスコアの分散の合計を2にしたいから

b2pb1<--0.7*90/(1.3*10)

b1<-sqrt(0.78/(1.8+0.2*b2pb1^2))

b2<-b2pb1*b1

つまり,

c(b1,b2)
## [1]  0.346489 -1.679139

の係数と合わせてb0に代入しておこう.

b0<-matrix(c(0.7,1.3,b1,b2),nrow=2)

が完成!

Py1_2<-b1*Px1[,1]+b2*Px1[,2]

と合わせて,

Py1<-cbind(Py1_1,Py1_2)

直交しているか確認.

var(Py1)%>%
  round(4)
##       Py1_1 Py1_2
## Py1_1  1.22  0.00
## Py1_2  0.00  0.78

ちゃんと直交しているし,分散も,合わせて2になっている.
ができた.

次に,を求める.

これは簡単.を45度逆回転させればよいだけ.

a1<-matrix(c(a,a,a,-a),ncol=2)
y1<-Py1%*%solve(a1)

の平均は,

apply(y1,2,mean)%>%
  round(4)
## [1] 0 0

の分散・共分散

var(y1)%>%
  round(4)
##      [,1] [,2]
## [1,] 1.00 0.22
## [2,] 0.22 1.00

ちゃんと標準化されてる!

のデータ完成.

のプロット図

の分布を確認する

plot(y1,asp=1)
grid()
text(x=y1[1:5,1],y=y1[1:5,2],1:5,pos = 4)
f_grid1(y1)

のスコアのプロット.のプロットを45度マイナスにした形となる.

plot(Py1,asp=1)
grid()
text(x=Py1[1:5,1],y=Py1[1:5,2],1:5,pos=4)
f_grid1(Py1)

・・・という感じ.

の主成分となっていることを確認しておく.

prc_y1<-prcomp(y1)
par(mfrow=c(1,2))
plot(prc_y1$x[,1]~Py1_1,asp=1)
grid()
plot(prc_y1$x[,2]~Py1_2,asp=1)
grid()

・・・と,まあ,データの作成に手間取ったが,通常は,「こんなデータがありました!」ということで,ここからが分析のスタートとなる.

OLS回帰とPCR

関係式

の関係式は以下

→主成分

とすると,

は51行×2列,も51行2列,は2行×2列)

→主成分

これも・・・

は51行×2列,も51行2列,は2行×2列)

主成分→主成分

とすると,

は51行×2列,も51行2列,は2行×2列)

OLS回帰

作成したでOLS回帰してみる.

への回帰は・・・

lm(y1[,1]~x1-1)%>%
  summary()
## Warning in summary.lm(.): essentially perfect fit: summary may be unreliable

## 
## Call:
## lm(formula = y1[, 1] ~ x1 - 1)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -4.951e-16 -6.125e-17  2.570e-18  9.978e-17  1.918e-15 
## 
## Coefficients:
##      Estimate Std. Error   t value Pr(>|t|)    
## x11 3.337e-01  7.162e-17 4.659e+15   <2e-16 ***
## x12 7.128e-01  7.162e-17 9.953e+15   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.039e-16 on 49 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 2.708e+32 on 2 and 49 DF,  p-value: < 2.2e-16

データにかく乱項を全く入れなかったので推定誤差は0.以下同様なので,係数だけ見る.

lm(y1~x1-1)%>%
  coef()%>%
  round(4)
##       [,1]    [,2]
## x11 0.3337  1.6663
## x12 0.7128 -1.3128

(被説明変数を2変数にすると,それぞれの回帰をいっぺんにやってくれる.1列目がy1[,1]に対する回帰の結果,2列目がy1[,2]`に対する回帰の結果)

一見,に対してが1.38,がその半分ぐらいの強さで影響しているように見えるが,プロットしてみると・・・

d1_1_OLS<-data.frame(x_2=x1[,2],
                     x_1=x1[,1],
                     y_1=y1[,1])

d1_2_OLS<-data.frame(x_2=x1[,2],
                     x_1=x1[,1],
                     y_2=y1[,2])

par(mfrow=c(1,2))

g3d1_1_OLS<-scatterplot3d(d1_1_OLS,
              type="h",
              angle=60,
              highlight.3d = T)
fit1_1_OLS<-lm(y_1~x_2+x_1-1,d1_1_OLS)
g3d1_1_OLS$plane3d(fit1_1_OLS)

g3d1_2_OLS<-scatterplot3d(d1_2_OLS,
              type="h",
              angle=60,
              highlight.3d = T)
fit1_2_OLS<-lm(y_2~x_2+x_1-1,d1_2_OLS)
g3d1_2_OLS$plane3d(fit1_2_OLS)

の相関が強く,平面を回帰するのは無理な感じがするが,誤差はないので,これはこれ.

という関係にあるので,を直接回帰したら,

であり,の係数はということだ.

これを上の例で計算してみると・・・

a1%*%b0%*%solve(a1)%>%
  round(4)
##        [,1]    [,2]
## [1,] 0.3337  1.6663
## [2,] 0.7128 -1.3128

OLSの結果と一致する.

主成分回帰(PCR)

のスコアととの回帰で,主成分回帰にあたる.

lm(y1~Px1-1)%>%
  coef()
##              [,1]     [,2]
## Px1PC1  0.7399795 0.249970
## Px1PC2 -0.2680918 2.106569

であり,の係数は

上の例だと・・・

b0%*%solve(a1)
##            [,1]     [,2]
## [1,]  0.7399795 0.249970
## [2,] -0.2680918 2.106569

一致する.

d1_1_PCR<-data.frame(Px2=Px1[,2],
                     Px1=Px1[,1],
                     y1=y1[,1])
d1_2_PCR<-data.frame(Px2=Px1[,2],
                     Px1=Px1[,1],
                     y2=y1[,2])

par(mfrow=c(1,2))

g3d1_1_PCR<-scatterplot3d(d1_1_PCR,
              type="h",
              angle=80,
              highlight.3d = T)
fit1_1_PCR<-lm(y1~Px2+Px1-1,d1_1_PCR)
g3d1_1_PCR$plane3d(fit1_1_PCR)

g3d1_2_PCR<-scatterplot3d(d1_2_PCR,
              type="h",
              angle=80,
              highlight.3d = T)
fit1_2_PCR<-lm(y2~Px2+Px1-1,d1_2_PCR)
g3d1_2_PCR$plane3d(fit1_2_PCR)

この場合,y_1y_2$はPx2方向の傾きだけで説明できるような形となっている.

PLS的な回帰

主成分分析のスコアどうしでOLS回帰してみる.

lm(Py1~Px1-1)%>%
  coef()%>%
  round(4)
##        Py1_1   Py1_2
## Px1PC1   0.7  0.3465
## Px1PC2   1.3 -1.6791

であり,の係数は定義通りで,. 上の例だと・・・

b0
##      [,1]      [,2]
## [1,]  0.7  0.346489
## [2,]  1.3 -1.679139

の第1主成分のスコアをの第1主成分のスコアのみで説明していて,の第2主成分のスコアをの第2主成分のスコアのみで説明している.

d1_1_scr<-data.frame(Px_2=Px1[,2],
                     Px_1=Px1[,1],
                     Py_1=Py1[,1])
d1_2_scr<-data.frame(Px_2=Px1[,2],
                     Px_1=Px1[,1],
                     Py_2=Py1[,2])

par(mfrow=c(1,2))

g3d1_1_scr<-scatterplot3d(d1_1_scr,
              type="h",
              angle=60,
              highlight.3d = T)
fit1_1_scr<-lm(Py_1~Px_2+Px_1-1,d1_1_scr)
g3d1_1_scr$plane3d(fit1_1_scr)


g3d1_2_scr<-scatterplot3d(d1_2_scr,
              type="h",
              angle=60,
              highlight.3d = T)
fit1_2_scr<-lm(Py_2~Px_2+Px_1-1,d1_2_scr)
g3d1_2_scr$plane3d(fit1_2_scr)

しかし,PLSはこうした主成分分析のスコアどうしのOLS回帰ではない

PLS回帰

PLSの考え方

のスコアとのスコアの内積が最大化するように各スコアを決める.

のスコアとのスコアの空間で主成分分析しているのと同じ考え方.

しかし,このスコアはをそれぞれ別々に主成分分析した場合のスコアとは違う.

が最大化するようにを決めた場合のスコア.

ただし,に制約しておく.

式で表すと以下ということになる.

s.t.

なので,ラグランジュ乗数法で

最大化の一階の条件は,で偏微分して・・・

で偏微分して・・・

であるなら,おお!これは特異値分解の必要条件そのものではないか!!!

・・・ということで,の特異値分解から得られる?

は・・・

xy1<-t(x1)%*%y1
round(xy1,4)
##         [,1]    [,2]
## [1,] 45.1963 30.8037
## [2,] 48.9877  1.0123

の特異値分解をすると・・・

svd_xy1<-svd(xy1)
print(svd_xy1)
## $d
## [1] 70.43287 20.77510
## 
## $u
##            [,1]       [,2]
## [1,] -0.7518087 -0.6593812
## [2,] -0.6593812  0.7518087
## 
## $v
##            [,1]       [,2]
## [1,] -0.9410458  0.3382792
## [2,] -0.3382792 -0.9410458

が得られる・・・はず

u1<-svd_xy1$u
v1<-svd_xy1$v
d1<-svd_xy1$d

(符号が逆方向で気持ち悪いので,逆にしておく)

u1<--u1
v1<--v1

ちなみに$dの内積

t(x1%*%u1)%*%y1%*%v1%>%
  round(4)
##         [,1]    [,2]
## [1,] 70.4329  0.0000
## [2,]  0.0000 20.7751

で割ると共分散).

の値は?

$d(つまり,の内積)に一致する.

が・・・

t(x1)%*%y1%*%v1%>%
  round(4)
##         [,1]     [,2]
## [1,] 52.9520  13.6987
## [2,] 46.4421 -15.6189

が・・・

u1%*%diag(d1)
##          [,1]      [,2]
## [1,] 52.95205  13.69871
## [2,] 46.44211 -15.61890

ぴったり.

が・・・

t(t(x1)%*%y1)%*%u1
##          [,1]      [,2]
## [1,] 66.28055 -7.027784
## [2,] 23.82597 19.550319

が・・・

v1%*%diag(d1)
##          [,1]      [,2]
## [1,] 66.28055 -7.027784
## [2,] 23.82597 19.550319

ぴったり.つまり,

ということでもある.

PCRとPLS

ノイズ(かく乱項)を入れていないので,を,それぞれ主成分分析したスコアと,PLSの考え方で得られたスコアは一致するのか?

par(mfrow=c(1,2))
plot(Px1[,1]~x1%*%u1[,1],asp=1)
grid()
plot(Px1[,2]~x1%*%u1[,2],asp=1)
grid()

おおよそは一致しているが,ぴったりとは一致していない

ではさらにぶれる.

par(mfrow=c(1,2))
plot(Py1[,1]~y1%*%v1[,1],asp=1)
grid()
plot(Py1[,2]~x1%*%v1[,2],asp=1)
grid()

PLSの考え方で最大化されたが・・・

xyPx1<-x1%*%u1
xyPy1<-y1%*%v1
t(xyPx1)%*%xyPy1%>%
  round(4)
##         [,1]    [,2]
## [1,] 70.4329  0.0000
## [2,]  0.0000 20.7751

主成分分析で,で個別に抽出した主成分の内積は・・・

t(Px1)%*%Py1%>%
  round(4)
##     Py1_1    Py1_2
## PC1    63  31.1840
## PC2    13 -16.7914

PLSの内積よりも小さい!

主成分分析では,それぞれに,最も分散の大きな主成分を抽出しているはずだが,との関係についての情報を含んでいないので,全体としては(PLSの意味での)最適ではないということ.

次のプロットは,左がPCRのスコア.右がPLSのスコア.

のスコアはあまり変わらないが・・・

par(mfrow=c(1,2))
plot(Px1,asp=1)
grid()
text(x=Px1[1:5,1],y=Px1[1:5,2],1:5,pos=4)
f_grid1(Px1)

plot(xyPx1,asp=1)
grid()
text(x=xyPx1[1:5,1],y=xyPx1[1:5,2],1:5,pos=4)
f_grid1(xyPx1)

のスコアは,よりに直交する形でとられているように見える.

par(mfrow=c(1,2))
plot(Py1,asp=1)
grid()
text(x=Py1[1:5,1],y=Py1[1:5,2],1:5,pos=4)
f_grid1(Py1)

plot(xyPy1,asp=1)
grid()
text(x=xyPy1[1:5,1],y=xyPy1[1:5,2],1:5,pos=4)
f_grid1(xyPy1)

ののスコアは直交していない

var(xyPx1)%>%
  round(4)
##         [,1]    [,2]
## [1,]  1.7932 -0.1043
## [2,] -0.1043  0.2068
var(xyPy1)%>%
  round(4)
##        [,1]   [,2]
## [1,] 1.1401 0.1696
## [2,] 0.1696 0.8599

スコアと,との関係も,45度回転ではない.

xya1_x<-lm(xyPx1~x1-1)%>%
  coef()
xya1_x
##          [,1]       [,2]
## x11 0.7518087  0.6593812
## x12 0.6593812 -0.7518087
xya1_y<-lm(xyPy1~y1-1)%>%
  coef()
xya1_y
##          [,1]       [,2]
## y11 0.9410458 -0.3382792
## y12 0.3382792  0.9410458

スコアどうしの関係も定義したとは異なる

xyb1<-lm(xyPy1~xyPx1-1)%>%
  coef()
xyb1
##             [,1]      [,2]
## xyPx11 0.8093295 0.1204335
## xyPx12 0.4083002 2.0696220
b0
##      [,1]      [,2]
## [1,]  0.7  0.346489
## [2,]  1.3 -1.679139

の主成分どうしの関係はこちら側で定義して作成したのに・・・
のスコア, のスコア, のスコア
もOLSで評価すると間違いないのに,PLSで一緒に分析すると,(PLSの意味で)もっとよい解が見つかった!・・・ということ.

ただし,PLSの係数を使ってに直接を回帰したの係数にあたるを計算すると,・・・

xya1_x%*%xyb1%*%solve(xya1_y)
##           y11       y12
## x11 0.3336750  1.666325
## x12 0.7128141 -1.312814

(途中が違っても)PCRの場合と(OLSとも)一致する.

lm(y1~x1-1)%>%
  coef()
##          [,1]      [,2]
## x11 0.3336750  1.666325
## x12 0.7128141 -1.312814

(ノイズが一切ないので当然だが)

Rのplsパッケージを使ってみる.

データ

データは,data.frameにするが,をそれぞれmatrixとして,2変数data.frameにするという荒業を使う.

matrixのx1y1をそのままdata.frameに突っ込むには,I()関数(そのまんまにしてね!という関数)を使う.

それぞれxyと名前を付けよう.

d01<-data.frame(x=I(x1),y=I(y1))

データ構造を見てみると・・・

str(d01)
## 'data.frame':    51 obs. of  2 variables:
##  $ x: 'AsIs' num [1:51, 1:2] 0 1 -1 1 -1 ...
##  $ y: 'AsIs' num [1:51, 1:2] -1.06e-18 1.05 -1.05 -3.79e-01 3.79e-01 ...

51行の2変数で,xyがあり,それぞれ,[1:51,1:2]の変数であると説明にある.

pls回帰

回帰には,plsr()関数を使う. 引数は,formula=の後に・・・
ncomp=は抽出する主成分の数(この場合,2変数しかないので,それ以上は無理)
data=はさっきつくったdata.frame
validation=は,クロスバリデーション(CV)の方法.“LOO”がLeave-one-out法.この場合,ややこしいので="none"(CV)にした.

fit01<-plsr(y~x,
            ncomp=2,
            data=d01,
            validation="none")

要約を表示しても,主成分の寄与率しか表示してくれない.

summary(fit01)
## Data:    X dimension: 51 2 
##  Y dimension: 51 2
## Fit method: kernelpls
## Number of components considered: 2
## TRAINING: % variance explained
##     1 comps  2 comps
## X     89.96      100
## Y1    98.00      100
## Y2    12.66      100

PLS回帰を使う人が,回帰係数よりも,たくさんの説明変数の中でどの変数との関係が強いか
に関心がある場合に使われているということか?

適合度とモデル選択

適合度は,Root Mean Squared Error of Prediction (RMSEP) で表す.
なんか長いが,要は,実測値と予測値との差の2乗平均の平方根をとったもの.当然小さい方がよいですね.

RMSEP関数が用意されている.

RMSEP(fit01)
## 
## Response: Y1 
## (Intercept)      1 comps      2 comps  
##   9.901e-01    1.401e-01    1.871e-16  
## 
## Response: Y2 
## (Intercept)      1 comps      2 comps  
##   9.901e-01    9.253e-01    8.450e-16

(Intercept)は,何も説明変数を使わなかったのRMSEP.ほぼ標準偏差ということになるが,ではなく,で割られるので,1となっていない.

sd(y1[,1])*sqrt(50/51)
## [1] 0.9901475

1 compsは,第1主成分だけを使った場合のRMSEP.

2 compsは,第2主成分まで使った場合のRMSEP.2変数しかないので,当然ここで0となる.

プロットもできる.

plot(RMSEP(fit01))

フィットの様子を確認することもできる.

横軸は抽出した主成分の数だが,1.2とか,余計な目盛がついているので,それは無視.主成分名を数字に置き換えてあるので,主成分の数が少ないとこうなるみたいだ.

そのまま,plot()関数を使うと,の実際値と予測値との相関を出力する.

第1主成分だけで予測した場合

plot(fit01,ncomp=1,line=TRUE,asp=1)
grid()

第2主成分も使って予測した場合

plot(fit01,ncomp=2,line=TRUE,asp=1)
grid()

ちなみに,予測値はpredict関数で出力できる.

predict(fit01,comps=1)%>%
  round(.,4)%>%
  head()
##        Y1      Y2
## 1  0.0000  0.0000
## 2  1.0432  0.3750
## 3 -1.0432 -0.3750
## 4  0.0683  0.0246
## 5 -0.0683 -0.0246
## 6 -0.7659 -0.2753

主成分のスコア

のスコアはscores()関数で

Scrx1<-scores(fit01)
Scrx1%>%
  round(4)%>%
  head()
##    Comp 1  Comp 2
## 1  0.0000  0.0000
## 2 -1.4112  0.0103
## 3  1.4112 -0.0103
## 4 -0.0924 -1.4166
## 5  0.0924  1.4166
## 6  1.0361 -0.1252

のスコアはから出したxyScrx1に(符号を除いて)一致するはずだが,

par(mfrow=c(1,2))
plot(Scrx1[,1]~xyPx1[,1],asp=1)
grid()
plot(Scrx1[,2]~xyPx1[,2],asp=1)
grid()

第2主成分で多少の違いが・・・

のスコアはYscores()関数で

Scry1<-Yscores(fit01)
Scry1%>%
  round(4)%>%
  head()
##    Comp 1  Comp 2
## 1  0.0000  0.0000
## 2 -1.4058  0.0103
## 3  1.4058 -0.0103
## 4 -0.8287 -1.4166
## 5  0.8287  1.4166
## 6  0.9710 -0.1252

については,比例はするが,スケールが全然ちがう. なぜ????

par(mfrow=c(1,2))
plot(Scry1[,1]~xyPy1[,1],asp=1)
grid()
plot(Scry1[,2]~xyPy1[,2],asp=1)
grid()

plsr()のスコアのばらつきをのスコアのばらつきに対応させている?

lm(Yscores(fit01)~scores(fit01)-1)%>%
  coef()%>%
  round(4)
##                     Comp 1 Comp 2
## scores(fit01)Comp 1 1.0000      0
## scores(fit01)Comp 2 0.5198      1

のスコアの分散は

var(Scry1)%>%
  round(4)
##        Comp 1 Comp 2
## Comp 1 1.8474 0.1043
## Comp 2 0.1043 0.2008

から求めたのスコアの分散は

var(xyPy1)%>%
  round(4)
##        [,1]   [,2]
## [1,] 1.1401 0.1696
## [2,] 0.1696 0.8599

plsrで求めたの第1主成分のスコアは,約1.27倍されている.

sqrt(1.8474/1.1401)
## [1] 1.272943

第2主成分のスコアは約0.48倍されている

sqrt(0.2008/0.8599)
## [1] 0.4832344

plsrの2つのスコアを,それぞれ割ってやると,の最大化から求めたのスコアとほぼ一致する.

Scry1_c<-sweep(Scry1,2,c(sqrt(1.8474/1.1401),sqrt(0.2008/0.8599)),FUN="/")
par(mfrow=c(1,2))
plot(Scry1_c[,1]~xyPy1[,1],asp=1)
grid()
plot(Scry1_c[,2]~xyPy1[,2],asp=1)
grid()

なんでそんなことするかは謎?

のスコアのプロット

plot関数に引数plottype="scores"で. scoreplot()関数でもOK

plot(fit01,plottype="scores",comps=1:2,asp=1)
grid()
text(x=Scrx1[1:5,1],y=Scrx1[1:5,2],1:5,pos=4)
f_grid1(Scrx1)

多変数の場合は,pairs()の図が出力されるようだ・・・

Yのスコアは関数が用意されていないが,この場合のスコアと全く同じなので,あんまり意味ない.

スコアととの相関のプロット

plot()関数のplottype="correlation"引数でも同じ.

corrplot(fit01)

cor(x1,Scrx1)%>%
  round(4)
##       Comp 1  Comp 2
## [1,] -0.9554 -0.2954
## [2,] -0.9416  0.3369

・・・というのをプロットしているようだ.左下がで,左上が.

も,第1主成分の方に相関が高いようだ・・・ということはわかる.

PLSのOLSやPCRとの違い

以上は全くノイズの無い世界の話.ノイズがなければ,OLSでを直接に回帰させても,あるいは,主成分回帰(PCR)で,いったんの主成分を抽出して,主成分をに回帰させても,それはどう解釈するかが違うだけで,情報の損失はない.

しかし,は2変数よりもすごく多いと,そうはいかなくなる.

まず,OLSだとの変数間でマルチコが生じる可能性が高まる.

だからといって,PCRを使うとすると,から主成分を抽出して,主な主成分を2個とか3個とか選んで使うことになる.

ところが,この主成分は,あくまでもの中だけの情報で抽出されるし,どれがどれが主な主成分かもの中だけで決まる.これにはの情報は入らない.

しかし,PLSだとの情報も入れて主成分を抽出するので,PCRよりもうまく主成分を抽出できる・・・という話.

以下,ノイズを入れてやってみよう.

2×2変数にノイズあり

データ作成

攪乱項を加える

2変数×2変数というのは変えないで,上記で作成したx1y1にちょっとだけかく乱項を入れてみる.

x2<-x1
set.seed(111)
x2[6:51,1]<-x1[6:51,1]+rnorm(46,mean=0,sd=0.5)
set.seed(124)
x2[6:51,2]<-x1[6:51,2]+rnorm(46,mean=0,sd=0.5)
x2<-apply(x2,2,scale)

にどのくらいノイズが入ったかプロットしてみる

par(mfrow=c(1,2))
plot(x2[,1]~x1[,1],asp=1)
grid()
plot(x2[,2]~x1[,2],asp=1)
grid()

このぐらい・・・

にもノイズを加える

y2<-y1
set.seed(125)
y2[6:51,1]<-y2[6:51,1]+rnorm(46,0,0.5)
set.seed(126)
y2[6:51,2]<-y2[6:51,2]+rnorm(46,0,0.5)

にどれだけノイズが入ったか見てみよう

par(mfrow=c(1,2))
plot(y2[,1]~y1[,1],asp=1)
grid()
plot(y2[,2]~y1[,2],asp=1)
grid()

このぐらい.

OLSやPCRとの比較

OLS回帰

作成したでOLS回帰してみる.

lm(y2~x2-1)%>%
  coef()
##          [,1]       [,2]
## x21 0.3545037  0.9186370
## x22 0.5574486 -0.6727718

真の値はこれ・・・

a1%*%b0%*%solve(a1)
##           [,1]      [,2]
## [1,] 0.3336750  1.666325
## [2,] 0.7128141 -1.312814

の係数以外はあてはまりがいまいち.こういうのは,乱数を多数発生させて,推定値の分布を見て判断すべきで,1つの例でどうのこうの言えることではないのだが,をプロットしてみると・・・

d2_1_OLS<-data.frame(x_1=x2[,1],
                   x_2=x2[,2],
                   y_1=y2[,1])
d2_2_OLS<-data.frame(x_1=x2[,1],
                   x_2=x2[,2],
                   y_2=y2[,2])
par(mfrow=c(1,2))
g3d1_1_OLS<-scatterplot3d(d2_1_OLS,
              type="h",
              angle=60,
              highlight.3d = T)
fit2_1_OLS<-lm(y_1~x_1+x_2-1,d2_1_OLS)
g3d1_1_OLS$plane3d(fit2_1_OLS)

g3d1_2_OLS<-scatterplot3d(d2_2_OLS,
              type="h",
              angle=60,
              highlight.3d = T)
fit2_2_OLS<-lm(y_2~x_1+x_2-1,d2_2_OLS)
g3d1_2_OLS$plane3d(fit2_2_OLS)


の相関に対して回帰平面が不安定そうなのは確認できる.

係数のt値では,そんな不安定であるようには見えないのが問題.

lm(y2~x2-1)%>%
  summary()%>%
  coefficients()
## Response Y1 :
##      Estimate Std. Error  t value     Pr(>|t|)
## x21 0.3545037  0.1118888 3.168358 2.640143e-03
## x22 0.5574486  0.1118888 4.982167 8.223860e-06
## 
## Response Y2 :
##       Estimate Std. Error   t value     Pr(>|t|)
## x21  0.9186370  0.1255169  7.318831 2.130580e-09
## x22 -0.6727718  0.1255169 -5.360010 2.226719e-06

PCR 回帰

を主成分分析

prc_x2<-prcomp(x2)
prc_x2
## Standard deviations (1, .., p=2):
## [1] 1.2830922 0.5947054
## 
## Rotation (n x k) = (2 x 2):
##             PC1        PC2
## [1,] -0.7071068  0.7071068
## [2,] -0.7071068 -0.7071068

主成分分析のスコア(方向を逆に)

Px2<-prc_x2$x
Px2[,1]<--prc_x2$x[,1]

のスコアで回帰

lm(y2~Px2-1)%>%
  coef()%>%
  round(4)
##           [,1]   [,2]
## Px2PC1  0.6448 0.1739
## Px2PC2 -0.1435 1.1253

真の値は・・・

b0%*%solve(a1)
##            [,1]     [,2]
## [1,]  0.7399795 0.249970
## [2,] -0.2680918 2.106569

これも1つの乱数だけで云々できないが,に対する係数がだいぶ外れているようだ.

プロットしてみると・・・

d2_1_PCR<-data.frame(Px_2=Px2[,2],
                     Px_1=Px2[,1],
                     y_1=y2[,1])
d2_2_PCR<-data.frame(Px_1=Px2[,1],
                     Px_2=Px2[,2],
                     y_2=y2[,2])
par(mfrow=c(1,2))

g3d1_1_PCR<-scatterplot3d(d2_1_PCR,
              type="h",
              angle=60,
              highlight.3d = T)
fit2_1_PCR<-lm(y_1~Px_2+Px_1-1,d2_1_PCR)
g3d1_1_PCR$plane3d(fit2_1_PCR)

g3d1_2_PCR<-scatterplot3d(d2_2_PCR,
              type="h",
              angle=60,
              highlight.3d = T)
fit2_2_PCR<-lm(y_2~Px_1+Px_2-1,d2_2_PCR)
g3d1_2_PCR$plane3d(fit2_2_PCR)


説明変数の相関がましになったことは確認できる.

lm(y1~Px2-1)%>%
  summary()%>%
  coefficients()
## Response Y1 :
##          Estimate Std. Error   t value     Pr(>|t|)
## Px2PC1  0.7360053 0.03174716 23.183344 4.577494e-28
## Px2PC2 -0.2756515 0.06849531 -4.024385 1.975324e-04
## 
## Response Y2 :
##         Estimate Std. Error  t value     Pr(>|t|)
## Px2PC1 0.2334963 0.07711698 3.027820 3.921328e-03
## Px2PC2 1.1032709 0.16638187 6.630956 2.475341e-08

PLS回帰

主成分の抽出

の特異値分解

svd_xy2<-svd(t(x2)%*%y2)
svd_xy2
## $d
## [1] 55.05462 19.84587
## 
## $u
##            [,1]       [,2]
## [1,] -0.7462628 -0.6656515
## [2,] -0.6656515  0.7462628
## 
## $v
##            [,1]       [,2]
## [1,] -0.9599648  0.2801207
## [2,] -0.2801207 -0.9599648
d2<-svd_xy2$d
u2<--svd_xy2$u
v2<--svd_xy2$v

の主成分

scr_x2<-x2%*%u2
scr_y2<-y2%*%v2

の主成分のスコア

plot(scr_x2,asp=1)
grid()
text(x=scr_x2[1:5,1],y=scr_x2[1:5,2],1:5,pos=4)
f_grid1(scr_x2)

の主成分のスコア

plot(scr_y2,asp=1)
grid()
text(x=scr_y2[1:5,1],y=scr_y2[1:5,2],1:5,pos=4)
f_grid1(scr_y2)

ちなみに,,別々に主成分分析して主成分を抽出すると...

prc_x2<-prcomp(x2)
prc_y2<-prcomp(y2)

からの主成分と比較すると,

par(mfrow=c(1,2))
plot(prc_x2$x[,1]~scr_x2[,1],asp=1)
grid()
plot(prc_x2$x[,2]~scr_x2[,2],asp=1)
grid()

の主成分分析は若干ずれていることがわかる.

の主成分については,もっとずれてくる.

par(mfrow=c(1,2))
plot(prc_y2$x[,1]~scr_y2[,1],asp=1)
grid()
plot(prc_y2$x[,2]~scr_y2[,2],asp=1)
grid()

スコアどうしで回帰する

lm(scr_y2~scr_x2-1)%>%
  coef()
##              [,1]       [,2]
## scr_x21 0.6767605 0.05014599
## scr_x22 0.1391105 1.11940449

真値は

xyb1
##             [,1]      [,2]
## xyPx11 0.8093295 0.1204335
## xyPx12 0.4083002 2.0696220

第2主成分が今一つだが,第1主成分は安定している.

プロットしてみると・・・

d2_1_PLS<-data.frame(Scrx_2=scr_x2[,2],
                     Scrx_1=scr_x2[,1],
                     Scry_1=scr_y2[,1])
d2_2_PLS<-data.frame(Scrx_2=scr_x2[,2],
                     Scrx_1=scr_x2[,1],
                     Scry_2=scr_y2[,2])

par(mfrow=c(1,2))

g3d2_1_PLS<-scatterplot3d(d2_1_PLS,
              type="h",
              angle=60,
              highlight.3d = T)
fit2_1_PLS<-lm(Scry_1~Scrx_2+Scrx_1-1,d2_1_PLS)
g3d2_1_PLS$plane3d(fit2_1_PLS)

g3d2_2_PLS<-scatterplot3d(d2_2_PLS,
              type="h",
              angle=60,
              highlight.3d = T)
fit2_2_PLS<-lm(Scry_2~Scrx_2+Scrx_1-1,d2_2_PLS)
g3d2_2_PLS$plane3d(fit2_1_PLS)

被説明変数がある程度ばらついており,の傾向もPCRより確認できやすい気がする?

lm(scr_y2~scr_x2-1)%>%
  summary()%>%
  coefficients()
## Response Y1 :
##          Estimate Std. Error    t value     Pr(>|t|)
## scr_x21 0.6767605 0.06498591 10.4139573 5.147398e-14
## scr_x22 0.1391105 0.13920572  0.9993159 3.225513e-01
## 
## Response Y2 :
##           Estimate Std. Error   t value     Pr(>|t|)
## scr_x21 0.05014599 0.07678037 0.6531095 5.167386e-01
## scr_x22 1.11940449 0.16447051 6.8061106 1.324866e-08

の第1主成分について,の第1主成分がだけが効いており,の第2主成分について,の第2主成分だけが効いている.

の主成分が2つの主成分が2つの場合,内積を最大化しようとすると,主成分1対1で対応するような気がする(要確認).

Rのplsr()関数で

データ作成

d02<-data.frame(x=I(x2),y=I(y2))

PLS回帰

fit02<-plsr(y~x,
            ncomp=2,
            data=d02,
            validation = "none")

要約

summary(fit02)
## Data:    X dimension: 51 2 
##  Y dimension: 51 2
## Fit method: kernelpls
## Number of components considered: 2
## TRAINING: % variance explained
##     1 comps  2 comps
## X    82.271   100.00
## Y1   65.185    66.29
## Y2    6.181    53.09

適合度とモデル選択

RMSEP(fit02)
## 
## Response: Y1 
## (Intercept)      1 comps      2 comps  
##      1.0116       0.5969       0.5874  
## 
## Response: Y2 
## (Intercept)      1 comps      2 comps  
##      0.9586       0.9285       0.6565

(Intercept)は,何もを使わなかったときのRMSEP.
1 compsは,第1主成分だけを使った場合のRMSEP.
2 compsは,第2主成分まで使った場合のRMSEP.2変数しかないが,ノイズがあるので0にはならない.

プロット

plot(RMSEP(fit02))


フィットの様子

第1主成分だけで予測した場合

plot(fit02,ncomp=1,line=TRUE,asp=1)
grid()

第2主成分も使って予測した場合

plot(fit02,ncomp=2,line=TRUE,asp=1)
grid()

主成分のスコア

の主成分のスコアは(最初の6個だけ)

px2<-scores(fit02)
px2%>%
  round(4)%>%
  head()
##    Comp 1  Comp 2
## 1 -0.1106 -0.0836
## 2 -1.3586 -0.1048
## 3  1.1375 -0.0624
## 4 -0.2181 -1.3323
## 5 -0.0030  1.1652
## 6  1.1188 -0.6712

の主成分のスコアは(最初の6個だけ)

py2<-Yscores(fit02)
py2%>%
  round(4)%>%
  head()
##    Comp 1  Comp 2
## 1  0.0619 -0.0882
## 2 -1.5840 -0.1580
## 3  1.7078 -0.0184
## 4 -0.6398 -2.7412
## 5  0.7637  2.5647
## 6  0.4764 -0.3743

のスコアはから出したscr_x2に(符号を除いて)ほぼ一致する

par(mfrow=c(1,2))
plot(px2[,1]~scr_x2[,1],asp=1)
grid()
plot(px2[,2]~scr_x2[,2],asp=1)
grid()

しかし,のスコアはとスケールが一致しない!

par(mfrow=c(1,2))
plot(py2[,1]~scr_y2[,1],asp=1)
grid()
plot(py2[,2]~scr_y2[,2],asp=1)
grid()

やはりのスコアの分散の調整が行われている.

lm(Yscores(fit02)~scores(fit02)-1)
## 
## Call:
## lm(formula = Yscores(fit02) ~ scores(fit02) - 1)
## 
## Coefficients:
##                      Comp 1      Comp 2    
## scores(fit02)Comp 1   1.000e+00  -8.270e-17
## scores(fit02)Comp 2   2.075e-01   1.000e+00

のスコアの分散は

py2_v<-var(py2)
py2_v%>%
  round(4)
##        Comp 1 Comp 2
## Comp 1 2.4023 0.1101
## Comp 2 0.1101 0.7093

の最大化から求めたのスコアの分散は

xy2_v<-var(scr_y2)
xy2_v%>%
  round(4)
##        [,1]   [,2]
## [1,] 1.0801 0.0414
## [2,] 0.0414 0.9009

plsr()の2つのスコアを,それぞれ割ってやると,の最大化から求めたのスコアと一致する.

py2_c<-sweep(py2,2,
      c(sqrt(py2_v[1,1]/xy2_v[1,1]),sqrt(py2_v[2,2]/xy2_v[2,2])),
      FUN="/")%>%
  round(4)
par(mfrow=c(1,2))
plot(py2_c[,1]~scr_y2[,1],asp=1)
grid()
plot(py2_c[,2]~scr_y2[,2],asp=1)
grid()

のスコアのプロット

scoreplot()でもOK

plot(fit02,plottype="scores",comps=1:2,asp=1)
grid()
text(x=px2[1:5,1],y=px2[1:5,2],1:5,pos=4)
f_grid1(px2)

スコアととの相関のプロット

plot()関数のplottype="correlation"引数でも同じ.

corrplot(fit02)

cor(x2,px2)
##          Comp 1     Comp 2
## [1,] -0.9180898 -0.3963724
## [2,] -0.8958415  0.4443737

・・・というのをプロットしているようだ.左下がで,左上が.

も,第1主成分の方に相関が高そう・・・ということはわかる.

の出力はloadings()で.

loadings(fit02)
## 
## Loadings:
##    Comp 1 Comp 2
## x1 -0.716 -0.666
## x2 -0.699  0.746
## 
##                Comp 1 Comp 2
## SS loadings     1.002  1.000
## Proportion Var  0.501  0.500
## Cumulative Var  0.501  1.001

Comp 1は,-0.716*x1+0.699*x2 だよと言っている.
Comp 2は,-0.666*x1+0.746*x2 だよと言っている.

2変数だけなので,図示する意味はあまりないが,loadingplot()でいける

loadingplot(fit02)

の1,2を数字として扱っているため,途中に1.2とか余計な目盛が入っているが,これは無視する.

は,Yloadings()で得られるはずだが

Yloadings(fit02)
## 
## Loadings:
##    Comp 1 Comp 2
## Y1 -0.644  0.180
## Y2 -0.188 -1.114
## 
##                Comp 1 Comp 2
## SS loadings     0.450  1.272
## Proportion Var  0.225  0.636
## Cumulative Var  0.225  0.861

・・・と,loadingが1を超えていたりとよくわからない???

との積をとれば,ちゃんとスコアと比例はするが,スケールが全然違う.

par(mfrow=c(1,2))
plot(y2%*%Yloadings(fit02)[,1]~fit02$Yscores[,1],asp=1)
grid()
plot(y2%*%Yloadings(fit02)[,2]~fit02$Yscores[,2],asp=1)
grid()

によるのスコアの分散は

svd_xy2<-t(x2)%*%y2%>%
  svd()
v_scr_xy2<-y2%*%svd_xy2$v%>%
  var()
v_scr_xy2
##            [,1]       [,2]
## [1,] 1.08010171 0.04140877
## [2,] 0.04140877 0.90089439

plsr()で出したスコアの分散は,が調整されている分だけ変わっている.

v_scr_pls2<-fit02$Yscores%>%
  var()
v_scr_pls2
##           Comp 1    Comp 2
## Comp 1 2.4023141 0.1100564
## Comp 2 0.1100564 0.7093003

PLSのLoadingsに,第1主成分の分散の比の平方根を掛けてやると,の1列目に一致する.

sqrt(v_scr_pls2[1,1]/v_scr_xy2[1,1])*Yloadings(fit02)[1,1]
## [1] -0.9599648
sqrt(v_scr_pls2[1,1]/v_scr_xy2[1,1])*Yloadings(fit02)[2,1]
## [1] -0.2801207

しかし,第2主成分の分散の比の平方根を掛けても,の2列目と近くはなるが,全くは一致しない.

sqrt(v_scr_pls2[2,2]/v_scr_xy2[2,2])*Yloadings(fit02)[1,2]
## [1] 0.159741
sqrt(v_scr_pls2[2,2]/v_scr_xy2[2,2])*Yloadings(fit02)[2,2]
## [1] -0.9880762
v2
##           [,1]       [,2]
## [1,] 0.9599648 -0.2801207
## [2,] 0.2801207  0.9599648

無相関の説明変数を加える

以上,2×2変数の場合で説明したが,PLSの強みは,数多くの説明変数の中から相関する主成分を抽出してくること.

8つのランダムな変数をに加えてみる.

データ作成

set.seed(123)
x3_c<-rnorm(51*8)%>%
  matrix(ncol=8)%>%
  scale()
x3<-cbind(x2,x3_c)

相関を確認してみる

pairs(x3,asp=1)


だけが相関していることを確認.

はそのまま使う

y3<-y2

plsr()関数用のデータセット

d03<-data.frame(x=I(x3),y=I(y3))

plsr()関数で

fit03<-plsr(y~x,
            ncomp=10,
            data=d03,
            validation="CV")

要約

summary(fit03)
## Data:    X dimension: 51 10 
##  Y dimension: 51 2
## Fit method: kernelpls
## Number of components considered: 10
## 
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
## 
## Response: Y1 
##        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## CV           1.032   0.7149   0.6844   0.6857   0.7112   0.6963   0.7012
## adjCV        1.032   0.7059   0.6780   0.6798   0.6996   0.6886   0.6935
##        7 comps  8 comps  9 comps  10 comps
## CV      0.7109   0.7188   0.7215    0.7220
## adjCV   0.7029   0.7100   0.7125    0.7129
## 
## Response: Y2 
##        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## CV          0.9777   0.9649   0.8875   0.8269   0.8100   0.7647   0.7612
## adjCV       0.9777   0.9642   0.8801   0.8096   0.8024   0.7559   0.7520
##        7 comps  8 comps  9 comps  10 comps
## CV      0.7596   0.7603   0.7599    0.7594
## adjCV   0.7506   0.7513   0.7509    0.7506
## 
## TRAINING: % variance explained
##     1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps  8 comps
## X    19.526    34.28    40.26    50.45    60.00    70.42    76.90    85.01
## Y1   60.092    66.74    67.06    69.49    70.04    70.06    70.07    70.09
## Y2    9.682    38.86    58.05    59.00    61.19    61.63    61.69    61.69
##     9 comps  10 comps
## X     94.23    100.00
## Y1    70.09     70.09
## Y2    61.69     61.69

主成分数をいくつにするか決める

その前にRMSEPって?

の実際値と推定値との誤差の平方平均の平方根がRMSEP.
しかし,計算してみると一致しない.

第1主成分までのRMSEPは

RMSEP(fit03,ncomp=1)
## 
## Response: Y1 
##        (Intercept)  1 comps
## CV           1.032   0.7149
## adjCV        1.032   0.7059
## 
## Response: Y2 
##        (Intercept)  1 comps
## CV          0.9777   0.9649
## adjCV       0.9777   0.9642

ところが,残差の平方平均の平方根は

fit03$residuals[,,1]^2%>%
  apply(.,2,mean)%>%
  sqrt()
##        Y1        Y2 
## 0.6390451 0.9109793

なんかちょっと違う.
実際には,クロスバリデーションで得られたRMSEPを計算している.

すべてのサンプルを使って推定すると(引数esimate="all"で,train1 compのところを見ると・・・),単純に計算したRMSEPと一致する.

RMSEP(fit03,ncomp=1,estimate="all")
## 
## Response: Y1 
##        (Intercept)  1 comps
## train        1.012   0.6390
## CV           1.032   0.7149
## adjCV        1.032   0.7059
## 
## Response: Y2 
##        (Intercept)  1 comps
## train       0.9586   0.9110
## CV          0.9777   0.9649
## adjCV       0.9777   0.9642

RMSEPで

plot(RMSEP(fit03), legendpos = "topright")

成分数を選ぶヒューリスティックスとして,selectNcomp()関数が用意されている.
使てみたいが,これはが1変数のときしか使えないようだ.

どんなものか知るだけのために,だけにPLS回帰してみよう.

fit03_1<-plsr(y[,1]~x,ncomp=10,data=d03,validation="CV")

selectNcomp()関数を使ってみる.判断基準は2通りあり,引数method="onsigma"と引数method="randomization"で使い分ける.

"onesigma"は,まず,RMSEPが最低の点を求め,そこから1標準誤差以内で最も少ない主成分数を採用するという方法.

"randomization"は,やはりRMSEPが最低の点を求め,そこから有意差がない主成分数を採用するという方法.

par(mfrow=c(1,2))
selectNcomp(fit03_1, method = "onesigma",plot=TRUE)
## [1] 1
selectNcomp(fit03_1, method = "randomization",plot=TRUE)

## [1] 1

両方とも1主成分が最適と出してきた.

についてもやってみよう.

fit03_2<-plsr(y[,2]~x,ncomp=10,data=d03,validation="CV")
par(mfrow=c(1,2))
selectNcomp(fit03_2, method = "onesigma",plot=TRUE)
## [1] 1
selectNcomp(fit03_2, method = "randomization",plot=TRUE)

## [1] 0

1か0なのだそうだ.

しかし,これはあくまで,を直接対応させた場合の結果.

予測値で

もとのと予測値とをプロットして確認

1主成分だけだとはいい感じでfitしているが,はいまいち.

plot(fit03, ncomp = 1, asp = 1, line = TRUE)
grid()

第2主成分まで使うとのfitはあんまりかわらんが,はだいぶよくなった.

plot(fit03, ncomp = 2, asp = 1, line = TRUE)
grid()

第3主成分まで使うと,両者ともあんまり変わらない.

plot(fit03, ncomp = 3, asp = 1, line = TRUE)
grid()


2主成分ということか?

のスコアのプロット

の第4主成分までスコアのプロットしてみると...

scoreplot(fit03,asp=1,comps = 1:4)
grid()

う~ん,この情報はどうなんだろう.

説明力の高いの変数を選択する

第1主成分だけでいくとして,のうち,どの変数を説明変数として採用するか・・・

(loading value)を確認する

loadingplot(fit03,legendpos = "topright",comps=1:4)

Comp 1loading valueを見ると,変数1と2だけが大きくぶれている.

Comp 2loading valueは,特に大きな値はないなあ.

これを見る限り,だけで1つの主成分を抽出し,を説明するのが適切か?という仮説が出される・・・たぶん.

ダメ押しで,スコアととの相関も.

次は,第1主成分と第2主成分のスコアとの各変数との相関

loading
valueのノルムで,左端のは1に近いが,他は0.5あたりか,それ未満.あんまり相関していないということ.

補足:シミュレーション

係数の不偏性を確かめる.

51×2のについて,OLS()とPCR(のスコア),それにの特異値分解から求めたスコア間のOLS回帰の結果を比較する.

3通りの係数を計算する関数

とりあえず1回の推定をする関数を定義する.

攪乱項はの両方に与え,その標準偏差は両者で同じとする.
各推定4つずつ,3種類で,12個の推定係数が出力.

f_pls01<-function(x1,y1,sd1=0.1){
  x2<-x1+matrix(rnorm(nrow(x1)*2,mean=0,sd=sd1),ncol=2)
  y2<-y1+matrix(rnorm(nrow(y1)*2,mean=0,sd=sd1),ncol=2)
  fit_OLS<-lm(y2~x2-1)%>%
    coef()
  prc01<-prcomp(x1)
  Px<-prc01$x
  fit_PCR<-lm(y2~Px-1)%>%
    coef()
  svd_xy<-svd(t(x2)%*%y2)
  Scrx<-x2%*%svd_xy$u
  Scry<-y2%*%svd_xy$v
  fit_PLS<-lm(Scry~Scrx-1)%>%
    coef()
  coef01<-c(fit_OLS,fit_PCR,fit_PLS)
  coef01
}

シミュレーション

攪乱項の標準偏差を0.2として,1000回推定

fit_sim01<-matrix(nrow=1000,ncol=12)
for(i in 1:1000){
  fit_sim01[i,]<-f_pls01(x1,y1,sd1=0.2)
}

fit_sim_OLS<-fit_sim01[,1:4]
fit_sim_PCR<-fit_sim01[,5:8]
fit_sim_PLS<-fit_sim01[,9:12]

OLS

のOLS推定

fit_OLS0<-a1%*%b0%*%solve(a1)
plot(density(fit_sim_OLS[,1]),xlim=c(-1.5,2),ylim=c(0,20),col=1)
par(new=T)
plot(density(fit_sim_OLS[,2]),xlim=c(-1.5,2),ylim=c(0,20),col=2)
par(new=T)
plot(density(fit_sim_OLS[,3]),xlim=c(-1.5,2),ylim=c(0,20),col=3)
par(new=T)
plot(density(fit_sim_OLS[,4]),xlim=c(-1.5,2),ylim=c(0,20),col=4)
abline(v=fit_OLS0[1,1],col=1)
abline(v=fit_OLS0[2,1],col=2)
abline(v=fit_OLS0[1,2],col=3)
abline(v=fit_OLS0[2,2],col=4)


垂直線が真値. 不偏性がない! 推定値の分散も大きい.

PCR

の主成分分析のスコア

fit_PCR0<-b0%*%solve(a1)
plot(density(fit_sim_PCR[,1]),xlim=c(-0.5,2.5),ylim=c(0,20),col=1)
par(new=T)
plot(density(fit_sim_PCR[,2]),xlim=c(-0.5,2.5),ylim=c(0,20),col=2)
par(new=T)
plot(density(fit_sim_PCR[,3]),xlim=c(-0.5,2.5),ylim=c(0,20),col=3)
par(new=T)
plot(density(fit_sim_PCR[,4]),xlim=c(-0.5,2.5),ylim=c(0,20),col=4)
abline(v=fit_PCR0[1,1],col=1)
abline(v=fit_PCR0[2,1],col=2)
abline(v=fit_PCR0[1,2],col=3)
abline(v=fit_PCR0[2,2],col=4)

不偏性がありそうだ.

の第1主成分の係数(黒と緑)は標準誤差が小さいが,第2主成分の標準が長が大きくなっている.

PLS

の特異値分解から,スコア間の関係

fit_PLS0<-xyb1
plot(density(fit_sim_PLS[,1]),xlim=c(-0.5,2.5),ylim=c(0,20),col=1)
par(new=T)
plot(density(fit_sim_PLS[,2]),xlim=c(-0.5,2.5),ylim=c(0,20),col=2)
par(new=T)
plot(density(fit_sim_PLS[,3]),xlim=c(-0.5,2.5),ylim=c(0,20),col=3)
par(new=T)
plot(density(fit_sim_PLS[,4]),xlim=c(-0.5,2.5),ylim=c(0,20),col=4)
abline(v=fit_PLS0[1,1],col=1)
abline(v=fit_PLS0[2,1],col=2)
abline(v=fit_PLS0[1,2],col=3)
abline(v=fit_PLS0[2,2],col=4)

の第1主成分(黒と緑)については,だいぶいいところまで来ているが,不偏性は怪しい.

第2主成分(赤と青)はだいぶ怪しい???

実際のPLSでは,このように一気にを算出するのではなく,まず,第1主成分を算出し,次いで,それを前提にを算出・・・というような手順を踏むようだ.(そっから先は,一般的なPLSの解説がたくさんある)

そもそも,PLSの関心が,係数の値そのものにあるのではなく,説明力の高い主成分と,その主成分と相関の高いの特定にある.plsr()関数などは,に合わせてのスケールそのものを変えてしまう.

複数のを主成分として扱えるところに魅力を感じたが,そこはあんまり意味なさそう.

たくさんのの中から,PLSで主成分とそれに相関の高いを抽出し,PCRで回帰するのが妥当な使い方のような気がする.

結論

  1. PLSはを最大化するようにを決める.
  2. このは,の特異値分解から得られる.
  3. 推定の考え方がOLSとは違うので,にかく乱項が加わった場合,係数の不偏性は怪しい.
  4. PLSの実際の利用は,多数の変数を持つからの説明力の高いの主成分を抽出し,それと相関の高いを特定すること.
  5. 変数間の係数に関心がある場合は,PLSで特定したでPCRをするのがよさそう.
  6. 実際のPLSは,の特異値分解という単純なものでなく,主成分を1つずつ算出していく複雑なアルゴリズムが複数ある.
  7. Rのパッケージplsにあるplsr()関数で出力されるのスコア(Yscores)はのスコアに合わせてスケールを調整しているので要注意.