統計ER

統計ソフトRの使い方を中心に、統計解析方法の解説をするブログ。ありそうでなかなか見つからないサンプルサイズ計算などニッチな方法について紹介しています。

主成分回帰 PCR と部分最小二乗回帰 PLS はどうやる?

ブログランキングに参加しています。
まずはぽちぽちっとお願いします。
↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
にほんブログ村 科学ブログ 数学へ

主成分回帰(Principal Component Regression, PCR)は、主成分分析と回帰分析の融合。主成分分析で情報の集約をして、変数を減らしてから回帰分析を行う方法。多重共線性が心配な変数同士が含まれていても、主成分得点に集約されるため問題がなくなる。

部分最小二乗回帰(Partial Least Squares Regression, PLS Regression)は、主成分回帰の発展版。独立変数は主成分回帰と同じ。PLSは独立変数だけでなく、従属変数にも新たな得点を想定する。独立変数も従属変数も違う世界に投射して潜在的な関係を見る。

主成分分析はこちら。

toukeier.hatenablog.com

データとパッケージの準備

データは、ISLRパッケージのCollegeデータを用いる。米国の777大学の出願者や学生数、書籍購入にかかるコストなどという18個のプロファイルデータ。

解析パッケージはplsパッケージ。主成分回帰(PCR)も部分最小二乗回帰(PLS)もできる。

最初の一回だけインストールする。

install.packages("ISLR")
install.packages("pls")

使うときはlibrary()で呼び出す。

library(ISLR)
library(pls)

例題は、登録する学生数(Enroll)を予測するモデルを作ること。

説明変数同士のVIFを見てみる。

toukeier.hatenablog.com

round(1/(1-cor(College[,c(-1,-4)])^2))

Apps(出願者数)とAccept(出願者受付人数)が10に近く多重共線性が疑われる。

> round(1/(1-cor(College[,c(-1,-4)])^2))
            Apps Accept Top10perc Top25perc F.Undergrad P.Undergrad Outstate Room.Board Books
Apps         Inf      9         1         1           3           1        1          1     1
Accept         9    Inf         1         1           4           1        1          1     1
Top10perc      1      1       Inf         5           1           1        1          1     1
Top25perc      1      1         5       Inf           1           1        1          1     1
F.Undergrad    3      4         1         1         Inf           1        1          1     1
P.Undergrad    1      1         1         1           1         Inf        1          1     1
Outstate       1      1         1         1           1           1      Inf          2     1
Room.Board     1      1         1         1           1           1        2        Inf     1
Books          1      1         1         1           1           1        1          1   Inf
Personal       1      1         1         1           1           1        1          1     1
PhD            1      1         1         1           1           1        1          1     1
Terminal       1      1         1         1           1           1        1          1     1
S.F.Ratio      1      1         1         1           1           1        1          1     1
perc.alumni    1      1         1         1           1           1        1          1     1
Expend         1      1         2         1           1           1        2          1     1
Grad.Rate      1      1         1         1           1           1        1          1     1
            Personal PhD Terminal S.F.Ratio perc.alumni Expend Grad.Rate
Apps               1   1        1         1           1      1         1
Accept             1   1        1         1           1      1         1
Top10perc          1   1        1         1           1      2         1
Top25perc          1   1        1         1           1      1         1
F.Undergrad        1   1        1         1           1      1         1
P.Undergrad        1   1        1         1           1      1         1
Outstate           1   1        1         1           1      2         1
Room.Board         1   1        1         1           1      1         1
Books              1   1        1         1           1      1         1
Personal         Inf   1        1         1           1      1         1
PhD                1 Inf        4         1           1      1         1
Terminal           1   4      Inf         1           1      1         1
S.F.Ratio          1   1        1       Inf           1      2         1
perc.alumni        1   1        1         1         Inf      1         1
Expend             1   1        1         2           1    Inf         1
Grad.Rate          1   1        1         1           1      1       Inf

学習データとテストデータに分ける。5分の4を学習データ。5分の1をテストデータにする。

set.seed(20180924)
sub <- sample(1:777, round(777*4/5))

主成分回帰 PCR はどうやるか?

主成分回帰 PCRは統計ソフトRでどうやるか?

plsパッケージのpcr()を使う。クロスバリデーションを行いながらモデルを作るために、validation="CV"を指定する。

学習データで、モデルを作る。

pcr.res.cv <- pcr(Enroll ~ ., data=College[sub,], validation="CV")

結果のサマリーを確認。Enrollを除き全部で17変数なので、17の合成変数が作成される。

> summary(pcr.res.cv)
Data:   X dimension: 622 17 
        Y dimension: 622 1
Fit method: svdpc
Number of components considered: 17

VALIDATION: RMSEP
Cross-validated using 10 random segments.
       (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps  8 comps
CV           868.6    791.8    217.1    207.8    207.5    191.5    192.6    192.9    188.5
adjCV        868.6    806.7    216.9    207.6    207.3    191.2    192.2    192.6    188.0
       9 comps  10 comps  11 comps  12 comps  13 comps  14 comps  15 comps  16 comps
CV       188.7     188.9     187.9     188.4     186.4     186.5     187.1     187.6
adjCV    188.1     188.3     187.3     187.8     185.8     185.9     186.4     186.9
       17 comps
CV        188.0
adjCV     187.3

TRAINING: % variance explained
        1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps  8 comps  9 comps
X         47.34    87.44    94.47    96.84    98.42    99.15    99.62    99.97   100.00
Enroll    20.93    93.90    94.39    94.50    95.40    95.43    95.43    95.76    95.76
        10 comps  11 comps  12 comps  13 comps  14 comps  15 comps  16 comps  17 comps
X         100.00    100.00    100.00    100.00    100.00    100.00    100.00    100.00
Enroll     95.76     95.81     95.82     95.89     95.89     95.91     95.91     95.91

主成分数ごとのMean Squared Error of Prediction (MSEP)を確認。

plot(MSEP(pcr.res.cv), legendpos="topright")

第二主成分でMSEPは激減してそれ以降は大きくは変わらない。

f:id:toukeier:20180924181636p:plain

計算上の最適な主成分の数を確認。答えは5。

selectNcomp(pcr.res.cv, method="onesigma", plot=T)

f:id:toukeier:20180924181828p:plain

主成分が5の時のMSEPを計算する。主成分が2の時のMSEPも計算してみる。

MSEP(pcr.res.cv, ncomp=5)
MSEP(pcr.res.cv, ncomp=2)

5のほうがMSEPが小さいことは小さい。

> MSEP(pcr.res.cv, ncomp=5)
       (Intercept)  5 comps
CV          754513    36689
adjCV       754513    36571
> MSEP(pcr.res.cv, ncomp=2)
       (Intercept)  2 comps
CV          754513    47134
adjCV       754513    47046

テストデータに適用して、MSEPを計算する。MSEPは学習データに比べてだいぶ大きくなった。予測性能は高くない。

> MSEP(pcr.res.cv, ncomp=5, newdata=College[-sub,])
(Intercept)      5 comps  
    1317758        95394  

部分最小二乗回帰 PLS はどうやるか?

Partial Least Squaresは、Projection to Latent Structureとも呼ばれる。Projection to Latent Structure(潜在的構造への投射)のほうが、計算方法を反映した名前だ。

使う関数はplsr()で、書き方はpcr()と同じ。

plsr.res.cv <- plsr(Enroll ~ ., data=College[sub,], validation="CV")
summary(plsr.res.cv)

PCRと比較して、第一主成分からEnrollの%Varianceが高い。当てはまりがいいことが期待される。

> summary(plsr.res.cv)
Data:   X dimension: 622 17 
        Y dimension: 622 1
Fit method: kernelpls
Number of components considered: 17

VALIDATION: RMSEP
Cross-validated using 10 random segments.
       (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps  8 comps
CV           868.6    227.6    216.8    203.8    196.6    198.2    195.9    195.0    195.4
adjCV        868.6    227.1    216.3    203.3    195.9    197.2    194.9    194.1    194.5
       9 comps  10 comps  11 comps  12 comps  13 comps  14 comps  15 comps  16 comps
CV       195.3     195.1     193.8     194.4     195.1     195.3     195.5     195.6
adjCV    194.4     194.0     192.8     193.4     194.0     194.2     194.4     194.5
       17 comps
CV        196.3
adjCV     195.2

TRAINING: % variance explained
        1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps  8 comps  9 comps
X         42.16    86.69    93.55    96.13    97.72    98.80    99.49    99.85   100.00
Enroll    93.65    94.32    95.05    95.57    95.68    95.76    95.76    95.76    95.76
        10 comps  11 comps  12 comps  13 comps  14 comps  15 comps  16 comps  17 comps
X         100.00    100.00     100.0    100.00    100.00    100.00    100.00    100.00
Enroll     95.87     95.89      95.9     95.91     95.91     95.91     95.91     95.91

第一主成分だけでMSEPはガクンと減少する。そのあとはなだらかに減少。

plot(MSEP(plsr.res.cv), legendpos="topright")

f:id:toukeier:20180925215922p:plain

最適な主成分数は4と示された。

selectNcomp(plsr.res.cv, method="onesigma", plot=T)

f:id:toukeier:20180925220035p:plain

主成分数4のときのMSEPを計算する。テストデータのMSEPと比較する。

MSEP(plsr.res.cv, ncomp=4)
MSEP(plsr.res.cv, ncomp=4, newdata=College[-sub,])

テストデータのMSEPはぐんと増えてしまって当てはまりはよくなかった。

> MSEP(plsr.res.cv, ncomp=4)
       (Intercept)  4 comps
CV          754513    38670
adjCV       754513    38384
> MSEP(plsr.res.cv, ncomp=4, newdata=College[-sub,])
(Intercept)      4 comps  
    1317758        89507  

まとめ

独立変数に多重共線性が疑われるときに有利な主成分回帰(PCR)と部分最小二乗回帰(PLS)。

今回はたまたまテストデータへの当てはまりが悪かっただけと思う。

普通の最小二乗法ならどちらかをやめなければならないとても相関が強くて、多重共線性が疑われる変数同士を同時に扱いたい場合は、どちらの変数もやめなくていいPCRやPLSをお試しあれ。