読者です 読者をやめる 読者になる 読者になる

東京に棲む日々

データ分析、統計、ITを勉強中。未だ世に出ず。

リッジ/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つのモデルはどれも線形回帰で、パラメータβは、以下の式によって推定される。

f:id:High_School_Student:20150208115607j:plain

(参考文献(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"は最小2乗法、それ以外の分布の場合は最尤法で行われるとのこと。

 


では、ペナルティと推定値のグラフ。

 plot(fitRidge1, xvar="norm", label=TRUE)

label=TRUE 番号ラベルを付ける。

f:id:High_School_Student:20150208121400j:plain

ラベルの凡例は以下。

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)

f:id:High_School_Student:20150208121625j:plain

この方がλとしているので、違和感なければ、こちらの方が見やすい気がする。

 

クロスバリデーションによって、最適なλを決定する。要するに、上で説明した、最適なモデルの複雑性をクロスバリデーションによって決定する、と言うことである。

cv.glmnet()で実施。デフォルトの分割数は10分割とのこと。任意の数を指定したい場合は引数nfoldsで指定する。

fitRidgeCV1 <- cv.glmnet( x=exp_vars, y=target_var, family="gaussian", alpha=0 )

 

 

MSEのプロット。
plot(fitRidgeCV1)

f:id:High_School_Student:20150208121857j:plain

左側の縦点線が、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)

f:id:High_School_Student:20150208122128j:plain


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)

f:id:High_School_Student:20150208122725j:plain

上軸のラベルを見ると、λが大きくなるにつれて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)

f:id:High_School_Student:20150208123042j:plain


λが最小となる時の、パラメータ推定値。

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)

f:id:High_School_Student:20150208123430j:plain

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)

f:id:High_School_Student:20150208141539j:plain

ラベルの凡例は以下。

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

http://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&frm=1&source=web&cd=1&ved=0CCQQFjAA&url=http%3A%2F%2Fwww.jstatsoft.org%2Fv33%2Fi01%2Fpaper&ei=IH6TVOHQMoHZmAWp1YCIBg&usg=AFQjCNEHiTrirITAlESPVPCAh83ht4xLzg&bvm=bv.82001339,d.dGY

(d) JMPマニュアル

http://www.jmp.com/japan/training/jmp_manual.shtml