PLS 回帰, Partial Least Squares Regression（部分的最小二乗回帰）とは、目的変数 Yを予測するために係数を最適化する手法のひとつです。
業務でこの手法を応用したツールを使っているのですが、導き出した結果を検証する環境を探していたところ、R の pls パッケージであれば十分な検証ができそうだということが判りました。そこで、この pls パッケージの使い方を備忘録的にまとめました。
動作環境は次の通りです。
- OS: Fedora 26 x86_64 (beta)
- R-core-3.4.0-2.fc26.x86_64
- rstudio-1.0.143-1.x86_64
ここでは参考サイト [1] に従って、自分で動作を確認して備忘録にすることが目的なので、統計的解釈に深くは立ち入っていないことをご了承下さい。また、あとで得た知見で書き直したり書き足したりすることもあります。
pls パッケージのインストール
R を起動して、以下のコマンドで pls パッケージをインストールします。
> install.packages("pls")
パッケージを ‘/home/bitwalk/R/x86_64-redhat-linux-gnu-library/3.4’ 中にインストールします
(‘lib’ が指定されていないため)
--- このセッションで使うために、CRAN のミラーサイトを選んでください ---
URL 'https://cran.ism.ac.jp/src/contrib/pls_2.6-0.tar.gz' を試しています
Content type 'application/x-gzip' length 809111 bytes (790 KB)
==================================================
downloaded 790 KB
* installing *source* package ‘pls’ ...
** パッケージ ‘pls’ の解凍および MD5 サムの検証に成功しました
** R
** data
*** moving datasets to lazyload DB
** inst
** preparing package for lazy loading
** help
*** installing help indices
converting help for package ‘pls’
finding HTML links ... 完了
biplot.mvr html
...
...
** building package indices
** installing vignettes
** testing if installed package can be loaded
* DONE (pls)
ダウンロードされたパッケージは、以下にあります
‘/tmp/RtmpBE2Z3b/downloaded_packages’
>
最初に、pls パッケージを library でロードしておきます。
> library(pls)
次のパッケージを付け加えます: ‘pls’
以下のオブジェクトは ‘package:stats’ からマスクされています:
loadings
>
使用するサンプル
以下のサンプルデータを用いますので、data でロードします。
> data(gasoline)
データの精度を4桁に設定します。
> options(digits = 4)
サンプルデータ gasoline の解析
> ?gasoline
gasoline package:pls R Documentation
Octane numbers and NIR spectra of gasoline
Description:
A data set with NIR spectra and octane numbers of 60 gasoline
samples. The NIR spectra were measured using diffuse reflectance
as log(1/R) from 900 nm to 1700 nm in 2 nm intervals, giving 401
wavelengths. Many thanks to John H. Kalivas.
Usage:
gasoline
Format:
A data frame with 60 observations on the following 2 variables.
‘octane’ a numeric vector. The octane number.
‘NIR’ a matrix with 401 columns. The NIR spectrum.
Source:
Kalivas, John H. (1997) Two Data Sets of Near Infrared Spectra
_Chemometrics and Intelligent Laboratory Systems_, *37*, 255-259.
サンプルの gasoline データ（オクタン価とガソリンの拡散反射の近赤外線分光データ、NIR スペクトル）を解析します。データは以下のような構造になっています。
> names(gasoline) [1] "octane" "NIR" > gasoline$octane [1] 85.30 85.25 88.45 83.40 87.90 85.50 88.90 88.30 88.70 88.45 88.75 88.25 [13] 87.30 88.00 88.70 85.50 88.65 88.75 85.40 88.60 87.00 87.15 87.05 87.25 [25] 86.85 88.65 86.60 86.00 86.10 86.50 86.30 84.40 84.70 84.60 84.50 88.10 [37] 85.25 88.40 88.20 88.40 88.55 88.35 88.20 85.30 88.50 88.25 88.00 88.85 [49] 88.45 88.70 88.10 87.60 88.35 85.10 85.10 84.70 87.20 86.60 89.60 87.10 > gasoline$NIR 900 nm 902 nm 904 nm 906 nm 908 nm 910 nm 912 nm 1 -0.050193 -0.045903 -0.042187 -0.037177 -0.033348 -0.031207 -0.030036 2 -0.044227 -0.039602 -0.035673 -0.030911 -0.026675 -0.023871 -0.022571 ... ... [ reached getOption("max.print") -- 58 行を無視しました ]
gasoline$octane の要素は、波長ごとに NIR スペクトルのリスト（ベクトル）になっています。もう少し詳しく調べてみます。
> dim(gasoline$NIR) [1] 60 401 > dimnames(gasoline$NIR) [[1]] [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15" "16" [17] "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28" "29" "30" "31" "32" [33] "33" "34" "35" "36" "37" "38" "39" "40" "41" "42" "43" "44" "45" "46" "47" "48" [49] "49" "50" "51" "52" "53" "54" "55" "56" "57" "58" "59" "60" [[2]] [1] "900 nm" "902 nm" "904 nm" "906 nm" "908 nm" "910 nm" "912 nm" "914 nm" [9] "916 nm" "918 nm" "920 nm" "922 nm" "924 nm" "926 nm" "928 nm" "930 nm" ... ... [393] "1684 nm" "1686 nm" "1688 nm" "1690 nm" "1692 nm" "1694 nm" "1696 nm" "1698 nm" [401] "1700 nm"
イメージしやすいように、gasoline$octane のデータについて、波長を横軸にとってプロットしてみました。もっとスマートなプロットの仕方があると思いますが、とりあえず、これでお許しください。
> y.values <- gasoline$NIR > x.label <- colnames(y.values) > x.value <- as.numeric(substring(x.label, 1, nchar(x.label) - 3)) > plot(x.value, y.values[1,], type = "n", main = "NIR spectra of gasoline samples", xlab = "wavelength", ylab = "log(1/R)") > for (i in 1:nrow(y.values)) lines(x.value, y.values[i, ], type = "l")
このデータ解析の目的は、
octane = f(NIR) = f(NIRwavelength1, NIRwavelength2, ...)
という関数関係を求めて、スペクトルからオクタン価を予測することなのですが、401 個の変数（NIR スペクトル）、60 組のデータでは、従来の重回帰分析を使おうとしても自由度が全然足りなくて解析できません。しかし、401 個のデータは互いに独立した関係にないので、これら NIR スペクトルの代わりに、主成分分析で主成分に分解し、目的変数と相関のある成分 compn から関数関係を求めます。
octane = g(compwavelength1, compwavelength2, ...)
まず、データをトレーニング用とテスト用の二つに分けます。
> gasTrain <- gasoline[1:50,] > gasTest <- gasoline[51:60,]
フィッティングはトレーニング用のデータ gasTrain を用いて次のようにします。
> gas1 <- plsr(octane ~ NIR, ncomp = 10, data = gasTrain, validation = "LOO")
このフィッティングでは 10 個の成分で行われ、一個抜き交差検証, leave-one-out (LOO) cross-validated predictions が実施されています。
フィッティングと検証結果は summary で確認できます。
> summary(gas1)
Data: X dimension: 50 401
Y dimension: 50 1
Fit method: kernelpls
Number of components considered: 10
VALIDATION: RMSEP
Cross-validated using 50 leave-one-out segments.
(Intercept) 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps 9 comps 10 comps
CV 1.545 1.357 0.2966 0.2524 0.2476 0.2398 0.2319 0.2386 0.2316 0.2449 0.2673
adjCV 1.545 1.356 0.2947 0.2521 0.2478 0.2388 0.2313 0.2377 0.2308 0.2438 0.2657
TRAINING: % variance explained
1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps 9 comps 10 comps
X 78.17 85.58 93.41 96.06 96.94 97.89 98.38 98.85 99.02 99.19
octane 29.39 96.85 97.89 98.26 98.86 98.96 99.09 99.16 99.28 99.39>
ちなみに、gas1 は以下のような構造になっています。
> names(gas1)
[1] "coefficients" "scores" "loadings" "loading.weights"
[5] "Yscores" "Yloadings" "projection" "Xmeans"
[9] "Ymeans" "fitted.values" "residuals" "Xvar"
[13] "Xtotvar" "fit.time" "ncomp" "method"
[17] "validation" "call" "terms" "model"
実測値と予測値の最小二乗誤差 RMSEP, Root Mean Square Error of Prediction が小さいほどフィッティングが良いと考えることができます。gas1 では成分 (comps) が 10 までの関係を調べてありますので、成分の数とフィッティング (RMSEP) との関係をプロットしてみます。
> plot(RMSEP(gas1), legendpos = "topright")
このプロットによると、octane を表現するのに二成分あればほぼ十分であるということが判ります。二成分のときの RMSEP は 0.2966 です。
成分の数が決まれば、その成分で予測した値と実測値との相関を確認することができます。
> plot(gas1, ncomp = 2, asp = 1, line = TRUE)
次に最初の3つの成分の得点 score の一覧表を作ります。
> plot(gas1, plottype = "scores", comps = 1:3)
この例では、特に目立った傾向が見られません。成分の得点一覧は次のように explvar で表示できます。
> explvar(gas1)
Comp 1 Comp 2 Comp 3 Comp 4 Comp 5 Comp 6 Comp 7 Comp 8 Comp 9 Comp 10
78.1707683 7.4122245 7.8241556 2.6577773 0.8768214 0.9466384 0.4921537 0.4723207 0.1688272 0.1693770
>
成分の負荷量 loading をプロットすることは、成分を解釈するためによく使われます。
> plot(gas1, "loadings", comps = 1:2, legendpos = "topleft", labels = "numbers", xlab = "nm")
解析を始める前にデータをトレーニング用とテスト用の二つに分けましたが、テスト用のデータ gasTest を用いて、二つの成分で予測値を計算するには次のようにします。
> predict(gas1, ncomp = 2, newdata = gasTest)
, , 2 comps
octane
51 87.94125
52 87.25242
53 88.15832
54 84.96913
55 85.15396
56 84.51415
57 87.56190
58 86.84622
59 89.18925
60 87.09116
>
テスト用のデータ gasTest は octane の実測値が含まれていますので、各成分数ごとの予測値との RMSEP を算出することができます。
> RMSEP(gas1, newdata = gasTest)
(Intercept) 1 comps 2 comps 3 comps 4 comps 5 comps 6 comps 7 comps 8 comps 9 comps 10 comps
1.5369 1.1696 0.2445 0.2341 0.3287 0.2780 0.2703 0.3301 0.3571 0.4090 0.6116
>
まとめ
参考サイト [1] の一部を、自分で確認してきましたが、これでおおよその使い方を掴めました。あとは、自分のデータを使って、業務で使っているツールが導き出す結果を検証できるように、もう少しこのパッケージを使い込んでいきたいと考えています。
参考サイト
にほんブログ村
