リッジ/Ridge回帰、Lasso回帰、Elastic Net (R - glmnet)
リッジ/Ridge回帰、Lasso回帰、Elastic Net に関して。
まず、モデルの複雑性とオーバーフィッティングに関して復習メモ。
複雑なモデル： バイアス(Bias)が小さく、バリアンス(Variance)が大きい
シンプルなモデル： バイアスが大きく、バリアンスが小さい
バイアスと言うのは、モデルによる予測値の平均が真の値からどれくらい乖離しているかの指標。小さいほうが良い。
バリアンスと言うのは、モデルの予測値自体のばらつき。予測の安定性みたいなもの。小さいほうが良い。
バイアスとバリアンスはトレードオフ関係にあり、モデル作成の際には複雑なモデルとシンプルなモデルの間の、ほどよく複雑なモデル（結果的に予測精度が高くなるモデル）を選ぶ必要がある。
各統計モデルによって、どのようにこの複雑性を調整するかが異なる。
例えば、重回帰であれば、説明変数を全部入れるだとか、たくさん入れると複雑なモデルになり、ステップワイズ法などにより少数の選択された説明変数のみを入れれば、シンプルなモデルになる。
|
モデル |
複雑性の調整方法 |
|
重回帰 |
投入する変数の数。多いと複雑。少ないとシンプル。 |
|
主成分回帰 |
作成する主成分の数。多いと複雑。少ないとシンプル。（主成分数=変数の数、が最大。この場合、全変数投入の重回帰に一致。） |
|
PLS |
作成する因子の数（因子数=変数の数、が最大。この場合、全変数投入の重回帰に一致。） |
|
ツリー（線形モデルでない） |
分岐する枝の数。多いと複雑。少ないとシンプル。 |
|
Ridge回帰 |
λ(Complexity Parameter)の大きさ。小さいと複雑。大きいとシンプル（Shrinkageが大きい）。 |
|
Lasso回帰 |
λ(Complexity Parameter)の大きさ。小さいと複雑。大きいとシンプル（Shrinkageが大きい）。 |
では、Ridge回帰、Lasso回帰、Elastic Netの説明。
Ridge回帰、Lasso回帰に関しては、参考文献(a)に詳しい説明がある。
例えば、p.64- には、主成分分析（特異値分解）と関連させてRidge回帰の係数の推定の特徴が説明されていたりと、興味深い。
これら3つのモデルはどれも線形回帰で、パラメータβは、以下の式によって推定される。
（参考文献(c), p.3 より）
式(2)（と(3)）が、ペナルティである。
α= 0 とした場合、Ridge回帰となる。
α= 1 とした場合、Lasso回帰となる。
0 < α< 1 の場合、Elastic Netである。
λはComplexity Parameterと呼ばれ、λが0だと、通常の最小二乗法に一致する。
パラメータβは式(1)から分かるとおり、大きいとペナルティの影響が強くなるので、パラメータβが小さく推定される。
このλにより、モデルの複雑性が調整される。小さいと複雑、大きいとシンプル（Shrinkageが大きい）なモデルとなる。
では、Rのglmnetパッケージを使ってみる。
library(glmnet)
外部データをインポート。
JMPのDiabetesデータを使用した。
LARSパッケージに、diabetesというデータがあるらしいが、おそらく同じものでないかと思う。
LASSO and Ridge regression - iAnalysis 〜おとうさんの解析日記〜
では、このデータを用いて同手法を説明しているが、このサンプルデータのデータ構造が良く分からないので、持ちデータをインポートして使用する。
data1に読み込んだデータを格納。
head(data1, 15)
Y Age Gender BMI BP Total.Cholesterol LDL HDL TCH LTG Glucose Validation
1 151 59 2 32.1 101 157 93.2 38 4.00 4.8598 87 Training
2 75 48 1 21.6 87 183 103.2 70 3.00 3.8918 69 Validation
3 141 72 2 30.5 93 156 93.6 41 4.00 4.6728 85 Training
4 206 24 1 25.3 84 198 131.4 40 5.00 4.8903 89 Training
5 135 50 1 23.0 101 192 125.4 52 4.00 4.2905 80 Training
6 97 23 1 22.6 89 139 64.8 61 2.00 4.1897 68 Training
7 138 36 2 22.0 90 160 99.6 50 3.00 3.9512 82 Training
8 63 66 2 26.2 114 255 185.0 56 4.55 4.2485 92 Validation
9 110 60 2 32.1 83 179 119.4 42 4.00 4.4773 94 Training
10 310 29 1 30.0 85 180 93.4 43 4.00 5.3845 88 Training
11 101 22 1 18.6 97 114 57.6 46 2.00 3.9512 83 Validation
12 69 56 2 28.0 85 184 144.8 32 6.00 3.5835 77 Training
13 179 53 1 23.7 92 186 109.2 62 3.00 4.3041 81 Training
14 185 50 2 26.2 97 186 105.4 49 4.00 5.0626 88 Validation
15 118 61 1 24.0 91 202 115.4 72 3.00 4.2905 73 Training
str(data1)
'data.frame': 442 obs. of 12 variables: $ Y : int 151 75 141 206 135 97 138 63 110 310 ... $ Age : int 59 48 72 24 50 23 36 66 60 29 ... $ Gender : int 2 1 2 1 1 1 2 2 2 1 ... $ BMI : num 32.1 21.6 30.5 25.3 23 22.6 22 26.2 32.1 30 ... $ BP : num 101 87 93 84 101 89 90 114 83 85 ... $ Total.Cholesterol: int 157 183 156 198 192 139 160 255 179 180 ... $ LDL : num 93.2 103.2 93.6 131.4 125.4 ... $ HDL : num 38 70 41 40 52 61 50 56 42 43 ... $ TCH : num 4 3 4 5 4 2 3 4.55 4 4 ... $ LTG : num 4.86 3.89 4.67 4.89 4.29 ... $ Glucose : int 87 69 85 89 80 68 82 92 94 88 ... $ Validation : Factor w/ 2 levels "Training","Validation": 1 2 1 1 1 1 1 2 1 1 ...
glmnet() 関数を使用するが、目的、説明変数はnumeric型のmatrixでないとエラーが出るので、変換しておく。カテゴリー型の変数は直接扱えない。
intをnumに変換したdata3を作る。
data3 <- data1
lstCol <- length(data1)
c_num <- 1
for(clm in data3[1:(lstCol-1)]){
data3[c_num] <- as.numeric(clm)
c_num <- c_num + 1
}
str(data3)
'data.frame': 442 obs. of 12 variables: $ Y : num 151 75 141 206 135 97 138 63 110 310 ... $ Age : num 59 48 72 24 50 23 36 66 60 29 ... $ Gender : num 2 1 2 1 1 1 2 2 2 1 ... $ BMI : num 32.1 21.6 30.5 25.3 23 22.6 22 26.2 32.1 30 ... $ BP : num 101 87 93 84 101 89 90 114 83 85 ... $ Total.Cholesterol: num 157 183 156 198 192 139 160 255 179 180 ... $ LDL : num 93.2 103.2 93.6 131.4 125.4 ... $ HDL : num 38 70 41 40 52 61 50 56 42 43 ... $ TCH : num 4 3 4 5 4 2 3 4.55 4 4 ... $ LTG : num 4.86 3.89 4.67 4.89 4.29 ... $ Glucose : num 87 69 85 89 80 68 82 92 94 88 ... $ Validation : Factor w/ 2 levels "Training","Validation": 1 2 1 1 1 1 1 2 1 1 ...
data.frameをmatrixにする。
Gender 、BMI、BP、Total.Cholesterol、LDL、HDL、TCH、LTG、Glucoseの9変数を説明変数とする。
exp_vars <- as.matrix(data3[3:11])
Yを目的変数。
target_var <- as.matrix(data3[1])
Ridge回帰を実行してみる。
fitRidge1 <- glmnet( x=exp_vars, y=target_var, family="gaussian", alpha=0 )
上でも説明したが、α=0の時がRidgeである。
fitRidge1$betaを実行すると、各Complexity Parameter（λ）におけるパラメータ推定値が確認できる。
デフォルトで、λは100分割されるように設定されている。nlambda引数に値を指定することにより、変更も可能なようである。
lambda.min.ratio引数で、どの程度の幅で分割するかを決められるようである。
fitRidge1$lambda[i-1] / fitRidge1$lambda[i] = 1.097499
1.097499^99 = 10000.21
(1.097499^99) * 0.0001 = 1.000021
0.0001は、glmnet()の引数であるlambda.min.ratioで指定される。
マニュアル（参考文献(b)）によれば、family="gaussian"は最小２乗法、それ以外の分布の場合は最尤法で行われるとのこと。
では、ペナルティと推定値のグラフ。
plot(fitRidge1, xvar="norm", label=TRUE)
label=TRUE 番号ラベルを付ける。
ラベルの凡例は以下。
1 Gender
2 BMI
3 BP
4 Total.Cholesterol
5 LDL
6 HDL
7 TCH
8 LTG
9 Glucose
軸上の目盛はnon zeroの変数の数。Ridgeはいずれも0にならない。
軸下の目盛、”L1 Norm”の意味はよく分からない。範囲は0-80となっている。
このグラフの方が有名だと思うが、横軸をλ（log(λ)）とする場合も描いてみる。
plot(fitRidge1, xvar="lambda", label=TRUE)
この方がλとしているので、違和感なければ、こちらの方が見やすい気がする。
クロスバリデーションによって、最適なλを決定する。要するに、上で説明した、最適なモデルの複雑性をクロスバリデーションによって決定する、と言うことである。
cv.glmnet()で実施。デフォルトの分割数は10分割とのこと。任意の数を指定したい場合は引数nfoldsで指定する。
fitRidgeCV1 <- cv.glmnet( x=exp_vars, y=target_var, family="gaussian", alpha=0 )
MSEのプロット。
plot(fitRidgeCV1)
左側の縦点線が、MSEが最小となるときのλの対数となる。
min(fitRidgeCV1$cvm)
3034.586
fitRidgeCV1$lambda.min
4.516003
log(fitRidgeCV1$lambda.min)
1.507627
右側の点線が、上のMSEが最小となるときのMSEの上側1seとなるときのλの対数。
選ばれたλを、ペナルティと推定値のグラフに描き入れてみる。
plot(fitRidge1, xvar="lambda", label=TRUE)
abline(v=log(fitRidgeCV1$lambda.min), lty=2)
MSEが最小となる時のλに対応するパラメータを出力するときは以下。
coef(fitRidgeCV1, s="lambda.min")
10 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) -235.51985766
Gender -20.86494597
BMI 5.43730943
BP 1.06492509
Total.Cholesterol -0.16518308
LDL -0.07801274
HDL -0.66474289
TCH 4.19250608
LTG 43.08980263
Glucose 0.33275805
MSEが最小となるときのMSEの上側1seとなるときのλに対応するパラメータを出力するときは以下。
coef(fitRidgeCV1, s="lambda.1se")
10 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) -1.482047e+02
Gender -9.306175e+00
BMI 3.600844e+00
BP 7.589827e-01
Total.Cholesterol 1.217223e-03
LDL -6.034893e-02
HDL -5.886791e-01
TCH 4.431805e+00
LTG 2.618412e+01
Glucose 4.785016e-01
推定したパラメータによる予測値を求める場合は以下。
pred_fitRidgeCV1 <- predict(fitRidgeCV1, s="lambda.min", newx=exp_vars)
glmnetでは、パラメータ推定の際に一度スケーリングが行われ、推定後に元のスケールに戻される処理が行われていることに注意。
Lasso回帰のあてはめ。
fitLasso1 <- glmnet( x=exp_vars, y=target_var, family="gaussian", alpha=1 )
glmnetで引数alpha=1と指定する。
Lidgeと同様にCVを実行。
fitLassoCV1 <- cv.glmnet( x=exp_vars, y=target_var, family="gaussian", alpha=1 )
MSEのプロット。
plot(fitLassoCV1)
上軸のラベルを見ると、λが大きくなるにつれてnon zero変数が減っていっているのが分かる。
各点線の位置。
MSEが最小。
log(fitLassoCV1$lambda.min)
-2.236981
1 SEの位置。
log(fitLassoCV1$lambda.1se)
1.949538
ペナルティと推定値のグラフにラインを加えたもの。
plot(fitLasso1, xvar="lambda", label=TRUE)
abline(v=log(fitLassoCV1$lambda.min), lty=2)
λが最小となる時の、パラメータ推定値。
coef(fitLassoCV1, s="lambda.min")
10 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) -297.31652686
Gender -22.36620204
BMI 5.63820180
BP 1.09805284
Total.Cholesterol -0.71269983
LDL 0.40363399
HDL -0.06213511
TCH 5.31599377
LTG 59.15181902
Glucose 0.27136481
予測の場合もRidgeと同様に。
pred_fitLassoCV1 <- predict(fitLassoCV1, s="lambda.min", newx=exp_vars)
Elastic Netも同様に実行してみる。
fitEN1 <- glmnet( x=exp_vars, y=target_var, family="gaussian", alpha=0.5 )
alpha=0.5としてとりあえず実行。
CVの実行。αが自動的に決定されるわけではない。
fitENCV1 <- cv.glmnet( x=exp_vars, y=target_var, family="gaussian", alpha=0.5 )
CVによるMSEのプロット。
plot(fitENCV1)
MSEが最小。
log(fitENCV1$lambda.min)
-3.776643
ペナルティと推定値のグラフにラインを加えたもの。
plot(fitEN1, xvar="lambda", label=TRUE)
abline(v=log(fitENCV1$lambda.min), lty=2)
予測値。
pred_fitENCV1 <- predict(fitENCV1, s="lambda.min", newx=exp_vars)
これらの手法の違いだが、JMPマニュアルに簡潔な記述があったので、それを引用しておく。
----- ここから -----
JMP Pro の「モデルのあてはめ」プラットフォームで提供されている「一般化回帰」手法は、パラメータ推定値を収縮（shrink）させた回帰分析を行います。この手法は、相関の高い説明変数が多数あるようなデータで役立ちます。用意されている方法のうち、Lasso と弾性ネットは、パラメータ推定値を収縮させるだけでなく、変数選択を行う側面もあります。
相関の高い説明変数が多数あるデータでは、一般的に、多重共線性の問題が生じます。最近のデータには、標本サイズよりも変数の数の方が多いこともよくあります。そのような場合は、変数を選択する必要があります。しかし、多重共線性や説明変数が多数あることにより、従来の方法ではうまく分析できないことがあります。
実験計画などの、小規模で相関があまりないようなデータでも、Lasso や弾性ネットは役立ちます。予測モデルの予測精度を上げる場合や、モデルに含める変数を選択する場合に役立ちます。
「一般化回帰」手法は、正則化した（罰則を課した）回帰分析を行います。正則化した回帰分析では、パラメータ推定値をゼロの方向に収縮（shrink）することで、モデルの予測精度を高めようとします。正則化によりパラメータ推定値を収縮させると、予測値にはバイアスが生じます。バイアスは増加するのですが、それと引き換えに、ばらつき（分散）が減れば、全体として予測誤差は小さくなります。
Lasso には2 つの欠点があります。第1 に、同じ程度に重要な変数のあいだに高い相関がある場合、Lasso はそこから変数を1 つだけしか選択しない傾向があります。第2 に、変数の個数p が標本サイズnよりも多い場合、Lasso は最大でn個しか説明変数を選択しません。
一方、弾性ネットは、重要な変数間に高い相関があっても、そこから複数の変数を選択しようとします。また、弾性ネットは、n < p のとき、n 個以上の説明変数を選択することもできます。ただし、弾性ネットは、通常、Lasso よりも計算時間が長くなります。
リッジ回帰は、罰則を課す回帰手法の中では初期のものの1 つです（Hoerl, 1962, Hoerl and Kennard,1970）。リッジ手法は係数を0 に収縮しないので、変数の選択は行いません。
「一般化回帰」手法では、Lasso と弾性ネットにおいて、適応型推定 (Adaptive Lasso) も行えます。適応型推定では、効果のある変数に対して、効果があまりない変数より少ない罰則が課されます。適応型推定は、モデルが預言的性質（oracle property）をもつようにするために提案されました。預言的性質とは、漸近的には効果のある説明変数がきちんと選択されるという性質です。より具体的には、次のような性質をもつことを指します。
? パラメータがゼロである説明変数を正確に識別すること。
? 得られるパラメータ推定値が、効果のある説明変数だけを使って構成したモデルのパラメータ推定値に、漸近的に一致すること。
----- ここまで -----
では、この記述のいくつかを試すため、70行のサンプルデータを作り、シミュレーションしてみる。
x1_1、x1_2、x1_3で高い相関を持たせる。
x1_1 <- rnorm(70)
x1_2 <- 0.9*x1_1 + 0.2*rnorm(70)
x1_3 <- 0.7*x1_1 + 0.2*rnorm(70)
x2_1、x2_2、x2_3で高い相関を持たせる。
x2_1 <- rnorm(70)
x2_2 <- 0.9*x2_1 + 0.2*rnorm(70)
x2_3 <- 0.7*x2_1 + 0.2*rnorm(70)
x3、x4は独立。
x3 <- rnorm(70)
x4 <- rnorm(70)
yを以下のように線形結合する。
y <- 0.7*x1_1 + 0.1*x2_1 + 0.7*x3 + 0.1*x4
cor(cbind(x1_1, x1_2, x1_3, x2_1, x2_2, x2_3, x3, x4, y))
x1_1 x1_2 x1_3 x2_1 x2_2 x2_3 x3 x4 y x1_1 1.0000000 0.9798883 0.9633291 0.29446205 0.30760303 0.269618825 0.28359103 0.273269025 0.8259803 x1_2 0.9798883 1.0000000 0.9424121 0.27852622 0.28981647 0.246803991 0.26735420 0.266495332 0.8023936 x1_3 0.9633291 0.9424121 1.0000000 0.24517903 0.25446480 0.208568264 0.29483674 0.237208146 0.8025598 x2_1 0.2944620 0.2785262 0.2451790 1.00000000 0.98061893 0.958538880 0.10406548 0.019824971 0.3218394 x2_2 0.3076030 0.2898165 0.2544648 0.98061893 1.00000000 0.941537391 0.06121275 0.015564267 0.3033108 x2_3 0.2696188 0.2468040 0.2085683 0.95853888 0.94153739 1.000000000 0.09870941 0.004176373 0.2987194 x3 0.2835910 0.2673542 0.2948367 0.10406548 0.06121275 0.098709407 1.00000000 0.064414685 0.7632235 x4 0.2732690 0.2664953 0.2372081 0.01982497 0.01556427 0.004176373 0.06441468 1.000000000 0.3026419 y 0.8259803 0.8023936 0.8025598 0.32183941 0.30331081 0.298719414 0.76322353 0.302641884 1.0000000
yはx1グループとx3と強い相関を持つ。
このデータに3手法を適用してみる。
all_x <- cbind(x1_1, x1_2, x1_3, x2_1, x2_2, x2_3, x3, x4)
Ridge回帰。
cvRidge <- cv.glmnet( x=all_x, y=y, family="gaussian", alpha=0 )
Lasso回帰。
cvLasso <- cv.glmnet( x=all_x, y=y, family="gaussian", alpha=1 )
Elastic Net（α=0.5）。
cvEN <- cv.glmnet( x=all_x, y=y, family="gaussian", alpha=0.5 )
Elastic Net（α=0.1）。
cvEN2 <- cv.glmnet( x=all_x, y=y, family="gaussian", alpha=0.1 )
これら4つのグラフを作成。
par(mfrow=c(2,2))
plot(cvRidge$glmnet.fit, xvar="lambda", label=TRUE, main="Ridge")
abline(v=log(cvRidge$lambda.min), lty=2)
plot(cvLasso$glmnet.fit, xvar="lambda", label=TRUE, main="Lasso")
abline(v=log(cvLasso$lambda.min), lty=2)
plot(cvEN$glmnet.fit, xvar="lambda", label=TRUE, main="Elastic Net, α=0.5")
abline(v=log(cvEN$lambda.min), lty=2)
plot(cvEN2$glmnet.fit, xvar="lambda", label=TRUE, main="Elastic Net, α=0.1")
abline(v=log(cvEN2$lambda.min), lty=2)
ラベルの凡例は以下。
1, 2, 3 = x1_1, x1_2, x1_3
4, 5, 6 = x2_1, x2_2, x2_3
7 = x3
8 = x4
Ridge回帰。
x1、x2はグループ内で固まっている。
また、x1はグループでyへの影響を分け合っているような形で推定値が出ている。
# coef(cvRidge, s="lambda.min")
Lasso回帰。
x1、x2はグループ内からyとの相関が強いものが極端に選ばれる形でパラメータ推定される。
# coef(cvLasso, s="lambda.min")
Elastic Net（α=0.5）。
グループから一つだけ変数が選ばれる傾向は、Lassoに近い推定が行われているようである。
# coef(cvEN, s="lambda.min")
Elastic Net（α=0.1）。
α=0.1と小さくしたが、これでもα=0.5に近い推定のようである。
# coef(cvEN2, s="lambda.min")
--------------------- 参考文献 ---------------------
(a) The Elements of Statistical Learning: Data Mining, Inference, and Prediction.
Trevor Hastie
Robert Tibshirani
Jerome Friedman
http://statweb.stanford.edu/~tibs/ElemStatLearn/
(b) ‘glmnet’ マニュアル
http://cran.r-project.org/web/packages/glmnet/glmnet.pdf
(c) Regularization Paths for Generalized Linear Models via Coordinate Descent
Jerome Friedman
Trevor Hastie
Rob Tibshirani
(d) JMPマニュアル
http://www.jmp.com/japan/training/jmp_manual.shtml
