RでPLSを初めて行う時、正直使い方がわからなくて困りました。
なぜかRjpWikiでPLSを検索しても引っかからないし...。
ネット上でRでのPLSの仕方を紹介しているサイトもあるんですが、少し昔のバージョンが紹介されているので、新しいバージョンの説明をここでして、皆さんの助けになればいいかなと思ってます。
library(pls)
data(yarn)
Rでのヘルプと同じように今回はRにもともとついているデータセットyarnを用いてPLSを行いたいと思います。
まず、library(pls)でパッケージplsを使えるようにします。パッケージplsがRに入ってない人はダウンロードしておきましょう。
続いて、data(yarn)でyarnのデータセットを扱いやすくしておきます。
なおRでPLSを行う際は分散が０の説明変数をもちいるとエラーとなります。必ず分散０の説明変数は取り除いておきましょう！
実際、私は初めてPLSを行った時このエラーで困り果てました(苦笑)。分散０の説明変数の取り除き方がわからないかたはR言語１を参照してください。
yarn.plsr <- plsr(density ~ NIR, 6, data = yarn, validation = "CV")
plsr()でPLSを行います、パッケージはplsで関数はplsr、ややこしいです(笑)
yarn.plsr <- plsr(density ~ NIR, 6, data = yarn, validation = "CV")
densityを目的変数、NIRを説明変数とするという意味です。
NIR1,NIR2,NIR3,・・・,NIR100と説明変数が100個あって、その中のNIR25,NIR40,NIR60だけを使う、という場合
あるいは、データセットyarnのdensity以外の変数をすべて用いるという場合はそれぞれ
yarn.plsr <- plsr(density ~ NIR25+NIR40+NIR60, 6, data = yarn, validation = "CV")
yarn.plsr <- plsr(density ~ ., 6, data = yarn, validation = "CV")
と、書けばOKです。
２行目は少し見にくいですがピリオドを書いています。２行目の方が重宝することが多いです。
yarn.plsr <- plsr(density ~ NIR, 6, data = yarn, validation = "CV")
この６という数字はPLSの潜在変数の数をいくつに指定するか、ということです。
RでPLSをこのように６で行った場合、１〜６のすべての潜在変数のパターンでPLSを行ってくれるようです。
つまり６というのは最大の潜在変数の数ということですね。当然ですが、説明変数以上の数をここで設定するとエラーとなるので注意しましょう。
yarn.pls <- plsr(density ~ NIR, 6, data = yarn, validation = "CV")
summary(yarn.pls)
Data: X dimension: 28 268
Y dimension: 28 1
Fit method: kernelpls
Number of components considered: 6
VALIDATION: RMSEP
Cross-validated using 10 random segments.
(Intercept) 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps
CV 27.46 4.826 4.242 2.052 0.8743 0.4864 0.4506
adjCV 27.46 4.440 4.325 2.034 0.8309 0.4759 0.4427
TRAINING: % variance explained
1 comps 2 comps 3 comps 4 comps 5 comps 6 comps
X 46.83 98.38 99.46 99.67 99.85 99.97
density 98.12 98.25 99.64 99.97 99.99 99.99
plsr()で計算した結果をyarn.plsに代入します。その後summary(yarn.pls)を行うことによって、PLSの結果を確認できます。
VALIDATION: RMSEP
Cross-validated using 10 random segments.
(Intercept) 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps
CV 27.46 4.826 4.242 2.052 0.8743 0.4864 0.4506
adjCV 27.46 4.440 4.325 2.034 0.8309 0.4759 0.4427
求めた各潜在変数に対して10-fold cross-validationを行い、RMSEP値を求めることによって評価することができます。 10-fold cross-validationはランダムにデータを10分割するので、得られる結果はPLSを行うたびに若干変わります。 RMSEP値は小さい程よい回帰式が得られているといえます。表の値を図として確認したい時には
plot(yarn.pls,"validation")
と、することで図にプロットできます。
今回の場合は、潜在変数の数が５あるいは６がうまく回帰できてるのではないかと判断できます。
PLSを行った時の各説明変数の係数を表示したい場合は次のように入力します。
yarn.pls$coefficients
, , 1 comps
density
NIR1 0.0632451417
NIR2 0.2454484945
NIR3 0.4759864013
NIR4 0.6553115941
NIR5 0.7010550584
NIR6 0.6216402781
NIR7 0.4500745938
NIR8 0.2051650494
NIR9 -0.0802673562
NIR10 -0.3721046802
・
・
・
各潜在変数の係数が表示されますが、潜在変数の数が１の場合,２の場合,３の場合…と各潜在変数の係数が表示されてしまうので、 自分が見たい係数を探すのが大変になることがあります。今回は潜在変数の数が５の時の係数を見たいとします。その場合は、
yarn.pls$coefficients[ , ,5]
というように、[ , , 見たい潜在変数] と入力することによって潜在変数の数が５の時の係数を表示できます。
図に係数を表示させたい場合は
plot(yarn.pls,"coefficients",type="l">
と、することでプロットすることができます。ただし、この図は潜在変数の数が最大の時(今回の場合は６)の係数を表示することしかできないようです。 潜在変数の数が５のときのプロットをしたい場合は、
hoge <- yarn.pls$coefficients[,,5]
plot(1:length(hoge),hoge,type="l",xlab="variable",ylab="regression coefficients",main="density")
abline(h=0,col="grey")
とすることで、同じ図を表示させることができます。
最後に、PLSで予測された目的変数を得るためには以下のようにすればいいようです。
yarn.pls$fitted.value[,,5]
こうすることで、潜在変数５の場合の予測値を得ることができます。
実測値に対する予測値をプロットするために用意されている関数がありますが、やはり潜在変数の数が最大である６の場合しかプロットできないようです。
潜在変数の数が５の場合のプロットをするためには、
plot(yarn.pls,"prediction")
hoge <- yarn.pls$fitted.value[,,5]
plot(yarn$density,hoge,xlab="measured",ylab="predicted")
abline(0,1,col="grey")
最初の１行は潜在変数の数が最大の場合(今回の場合は６)のプロットをするときのものです。
５の場合は２行目以降を使えばOKです。最後の行でy=xの直線を引いています。予測がうまくいってるかの基準になっていいと思います。
必要がない場合は飛ばしてもらって構いません。
以上で、簡単ではありますがPLSの説明は終わりです。
公開日：４月２２日