統計ER

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

リッジ、ラッソ、エラスティックネットで高性能な予測モデルを作る!

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

線形回帰モデルは、変数がたくさんあればあるほど、当てはまりのいい回帰式ができる。

当てはまりのいい回帰式の欠点は、新しいデータでの予測に使えないことだ。

予測性能に優れるモデルを作るのがRidge(リッジ)回帰とLASSO(ラッソ)回帰、Elastic net(エラスティックネット)だ。

リッジ、ラッソ、エラスティックネットは何をやっているか?

線形回帰モデルは、係数\beta(パラメータ)を推定するときに最小二乗法を用いる。

通常の最小二乗法は、従属変数の実測値とモデルから計算された値との差=残差を二乗して足し合わせる。この合計が最小になるパラメータのセットを見つける。

Ridge, LASSO, Elastic netは残差の二乗和に以下の式を足したものを最小にする方法を取る。「罰則」と呼ぶ。

 \displaystyle \lambda \left( \frac{1-\alpha}{2}||\beta||_2^{2}+\alpha||\beta||_1 \right)

 \displaystyle = \lambda \sum_{j=1}^{p} \left(\frac{1-\alpha}{2}\beta_j^{2} + \alpha|\beta_j| \right)

Ridgeは最初の項だけ、LASSOは後の項だけ、Elastic netはどちらの項も使う方法だ。

\alphaが0の時はRidge、\alphaが1の時はLASSO、0と1の間のときはElastic netになる。

最初の項は\beta^{2}の和になっている。後の項は\betaの絶対値の和になっている。

いずれも\betaが多くなる、大きくなると、全体が大きくなる。

つまり、どんなに残差二乗和が小さくなるようにしても、この項たちが大きくなったのでは、目的が達成できない。「罰則」と呼ばれるゆえんだ。

\betaを小さくする、いくつかを0にするなどして、「罰則」を小さくする。

結果として、当てはまりすぎ(過学習、overfitting)を防ぐことになる。   ちなみに||\beta||_2^{2}とか、||\beta||_1はノルムと呼ばれる。距離を表していると考えればよい。

 \displaystyle ||\beta||_2^{2} = \sqrt{|\beta|_1^{2}+|\beta|_2^{2}+ \cdots +|\beta|_n^{2}}

 \displaystyle ||\beta||_1 = |\beta|_1+|\beta|_2+ \cdots +|\beta|_n

Ridge、LASSO、Elastic netの名前の由来、意味は?

気になる人のために名前の由来、意味を記しておこう。

Ridgeは、一般的な意味は盛り上がった場所、例えば山の背、尾根など。二乗和の項を入れることで、上に凸型の二次反応関数に、数学的類似性があるところから名づけられた。

LASSOは、Least Absolute Shrinkage and Selection Operatorの頭文字語。係数の絶対値(Absolute)を用いる方法で、係数を小さくして(Shrinkage)罰則を小さくし、要らない変数の係数をゼロにして変数を選択する(Selection)方法だ。

Elastic netは、まるで伸縮性の投網のように「大きな魚(重要な係数)全部」を残すイメージからつけられた名前だ。

リッジ、ラッソ、エラスティックネットを統計ソフトRでやってみる

必要なパッケージの準備

最初に一回インストール。

glmnetパッケージとISLRパッケージ。

install.packages("glmnet")
install.packages("ISLR")

ともに呼び出しておく。

library(glmnet)
library(ISLR)

 

例題のデータを準備

ISLRパッケージ中のCollegeデータを使う。

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

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

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

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

リッジ、ラッソ、エラスティックネットの前に従来の最小二乗法

Enrollをその他の変数全部で予測する回帰式を作る。

lm.res <- lm(Enroll ~ ., data=College[sub,])
summary(lm.res)

6つの変数がp<0.05で統計学的有意。

> summary(lm.res)

Call:
lm(formula = Enroll ~ ., data = College[sub, ])

Residuals:
     Min       1Q   Median       3Q      Max 
-1577.18   -63.47   -12.88    45.64  1581.09 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 219.946039  87.176801   2.523 0.011892 *  
PrivateYes  -20.716365  29.254284  -0.708 0.479126    
Apps         -0.040325   0.007908  -5.099 4.57e-07 ***
Accept        0.162698   0.014738  11.040  < 2e-16 ***
Top10perc     4.381100   1.270779   3.448 0.000605 ***
Top25perc    -2.206551   0.943228  -2.339 0.019642 *  
F.Undergrad   0.142490   0.004467  31.901  < 2e-16 ***
P.Undergrad  -0.029804   0.008131  -3.666 0.000269 ***
Outstate     -0.004381   0.004120  -1.063 0.288023    
Room.Board   -0.018011   0.010351  -1.740 0.082361 .  
Books         0.018743   0.049424   0.379 0.704654    
Personal     -0.004377   0.014736  -0.297 0.766533    
PhD           0.353353   1.024941   0.345 0.730399    
Terminal     -1.275895   1.139472  -1.120 0.263276    
S.F.Ratio    -0.842953   2.686361  -0.314 0.753789    
perc.alumni   1.920409   0.879343   2.184 0.029353 *  
Expend        0.001368   0.002644   0.517 0.605075    
Grad.Rate     0.277165   0.639241   0.434 0.664745    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 199.2 on 604 degrees of freedom
Multiple R-squared:  0.9529,    Adjusted R-squared:  0.9516 
F-statistic: 718.9 on 17 and 604 DF,  p-value: < 2.2e-16

モデルから得られる予測値と実測値のかい離をmean squared error of prediction(MSEP)で確認する。

mean((predict(lm.res)-College[sub,]$Enroll)^2)

> mean((predict(lm.res)-College[sub,]$Enroll)^2)
[1] 38530.83

学習データから得られたモデルをテストデータに適用して、テストデータで得られる予測値と実測値のかい離をMSEPで確認する。

lm.test <- predict(lm.res, newdata=College[-sub,])
mean((lm.test - College[-sub,]$Enroll)^2)

> mean((lm.test - College[-sub,]$Enroll)^2)
[1] 47213.26

MSEPはテストデータで大きくなった。これは学習データでの予測モデルの性能があまり高くないことを示している。

リッジで解析するとどうなるか?

glmnetパッケージのglmnet()を使う。データは行列で与える。data.matrix()を使う。説明変数, 従属変数の順番に指定し、alpha=0を指定する。

ridge.res <- glmnet(data.matrix(College[sub,-4]), College[sub,]$Enroll, alpha=0)
mean((predict(ridge.res, newx=data.matrix(College[sub, -4]))-College[sub,]$Enroll)^2)
ridge.test <- predict(ridge.res, newx=data.matrix(College[-sub,-4]))
mean((ridge.test - College[-sub,]$Enroll)^2)

学習データで作られたモデルを、学習データとテストデータのMSEPで見てみると、テストデータでははるかにMSEPが大きくなっている。あまり予測性能はよくない。

> mean((predict(ridge.res, newx=data.matrix(College[sub, -4]))-College[sub,]$Enroll)^2)
[1] 442744.8

> mean((ridge.test - College[-sub,]$Enroll)^2)
[1] 558366.1

予測性能を上げる方法として、クロスバリデーション(cross-validation)という方法がある。学習データをさらに分割して、テストセット1つと学習セットとして繰り返し、最適な\lambdaを見つける方法。

学習データにcross-validation付きのglmnetを適用して、結果をプロットする。

ridge.res.cv <- cv.glmnet(data.matrix(College[sub,-4]), College[sub,]$Enroll, alpha=0)
plot(ridge.res.cv)

f:id:toukeier:20180923152237p:plain

MSEP(Y軸)が最も小さい値を確認する。左側の縦破線が通っている赤いドットのY軸値。

> min(ridge.res.cv$cvm)
[1] 61345.74

MSEPが最も小さいときの\lambda\log(\lambda)を確認する(左側の縦破線のX軸値)。

> ridge.res.cv$lambda.min
[1] 95.77041
> log(ridge.res.cv$lambda.min)
[1] 4.561954

MSEPが最も小さいときの\lambdaに対応した回帰係数を確認する。

> coef(ridge.res.cv, s="lambda.min")
18 x 1 sparse Matrix of class "dgCMatrix"
                        1
(Intercept)  2.234662e+02
Private     -7.114514e+01
Apps         1.738649e-02
Accept       1.106883e-01
Top10perc    2.040518e+00
Top25perc   -3.542996e-01
F.Undergrad  1.035560e-01
P.Undergrad  8.384248e-03
Outstate    -4.248129e-03
Room.Board  -3.121038e-02
Books        2.823821e-02
Personal     1.564955e-02
PhD          5.052552e-01
Terminal     2.447818e-02
S.F.Ratio    2.254079e+00
perc.alumni  1.854235e+00
Expend      -2.158011e-04
Grad.Rate    1.225929e-01

学習データでMSEPが最も小さいときの\lambdaのモデルをテストデータに適用する。

ridge.test.cv <- predict(ridge.res.cv, s="lambda.min", newx=data.matrix(College[-sub,-4]))
mean((ridge.test.cv-College[-sub,]$Enroll)^2)

テストデータの予測値と実測値のかい離をMSEPで求める。

> mean((ridge.test.cv-College[-sub,]$Enroll)^2)
[1] 44719.1

学習データのMSEPは、以下の通り。

> min(ridge.res.cv$cvm)
[1] 61345.74

テストデータのMSEPのほうが小さくなっているのがわかる。予測性能がよいことを示している。

ラッソで解析するとどうなるか?

glmnet()でalpha=1と指定する。最初からクロスバリデーションを使う。

lasso.res.cv <- cv.glmnet(data.matrix(College[sub,-4]), College[sub,]$Enroll, alpha=1)
plot(lasso.res.cv)

f:id:toukeier:20180923154452p:plain

MSEPが最小値を表示させる(Y軸)。MSEPが最小値の\log(\lambda)を表示させる(X軸)。MSEP最小値の\lambdaのときの回帰係数を表示させる。

> min(lasso.res.cv$cvm)
[1] 46512.83
> log(lasso.res.cv$lambda.min)
[1] -0.5781604
> coef(lasso.res.cv, s="lambda.min")
18 x 1 sparse Matrix of class "dgCMatrix"
                        1
(Intercept) 223.249334574
Private     -17.323492035
Apps         -0.034361012
Accept        0.152731795
Top10perc     3.758399461
Top25perc    -1.812982422
F.Undergrad   0.142965370
P.Undergrad  -0.029337386
Outstate     -0.003177355
Room.Board   -0.018717120
Books         0.012830680
Personal     -0.003259145
PhD           0.191660710
Terminal     -1.050251664
S.F.Ratio    -0.547854820
perc.alumni   1.834686143
Expend        0.001059255
Grad.Rate     0.177922346

学習データでもっともすぐれたモデルで、テストデータ予測をする。MSEPを計算する。

lasso.test.cv <- predict(lasso.res.cv, s="lambda.min", newx=data.matrix(College[-sub,-4]))
mean((lasso.test.cv-College[-sub,]$Enroll)^2)

学習データでのMSEPとほぼ同じ値。大きく増加していないので、予測に全くつかえないモデルではない。

> mean((lasso.test.cv-College[-sub,]$Enroll)^2)
[1] 46711.87
> min(lasso.res.cv$cvm)
[1] 46512.83

エラスティックネットはどうやる?

 0 \lt \alpha \lt 1がエラスティックネット elastic netだ。RidgeとLASSOの中間。いいとこどり。

 \alphaを振って試してくれる関数がglmnetUtilsパッケージのcva.glmnet()だ。

最初一回だけインストールする。使うときはlibrary()で呼び出す。

install.packages("glmnetUtils")
library(glmnetUtils)

cva.glmnet()はcv.glmnet()と書き方は一緒。\alphaはデフォルトで11個に振ってくれる。

elnet.res.cv <- cva.glmnet(data.matrix(College[sub,-4]), College[sub,]$Enroll)
plot(elnet.res.cv)

 \log(\lambda)とMSEPのプロットは以下の通り。\alphaが11個に振られているので、11本のグラフになる。 f:id:toukeier:20180923223058p:plain

\alphaを確認すると以下の通り。

> elnet.res.cv$alpha
 [1] 0.000 0.001 0.008 0.027 0.064 0.125 0.216 0.343 0.512 0.729 1.000

最適な\alphaを選ぶためMSEPがもっとも小さい\alphaを探す。各\alphaのうち最も小さいMSEPを取りだし、その中でもっとも小さいMSEPを取り出したところ、\alpha=1がもっともMSEPが小さかった。つまり、LASSOがもっとも当てはまりがよいモデルであった。

> sapply(sapply(elnet.res.cv$modlist, "[[", "cvm"),min)
 [1] 61274.57 61208.48 47948.64 47276.72 47268.60 47292.46 47324.64 47347.99 47177.38
[10] 47006.23 46913.66
> min(sapply(sapply(elnet.res.cv$modlist, "[[", "cvm"),min))
[1] 46913.66

\alpha=1のときの、MSEPが最小になる\lambdaのモデルを使って、テストデータで予測する。

elnet.test.cv <- predict(elnet.res.cv, alpha=1, s="lambda.min", newx=data.matrix(College[-sub,-4]))
mean((elnet.test.cv-College[-sub,]$Enroll)^2)

MSEPは学習データよりも小さくなり、予測性能が高いと言える。

> mean((elnet.test.cv-College[-sub,]$Enroll)^2)
[1] 45275.5

まとめ

学習データに当てはまりすぎない、予測性能に優れたモデルを作る、Ridge, LASSO, Elastic netを紹介した。

Ridgeがいいか、LASSOがいいか、は、Elastic netで\alphaを振ってみて、選ぶのがよい。

glmnetUtilsパッケージのcva.glmnet()は痒いところに手が届く関数。見つけられてよかった。