About connecting the dots.

data science related trivial things

glmnetで正則化を試してみる

タイトルの通り,よく考えたら今までL1/L2正則化を知識としては知ってるけど,実際に試したことはなかったことに気がついたので試してみましたよという話.L1/L2正則化にの理屈については,TJOさんのエントリとか,unnounnoさんのエントリとかをみてもらえれば良いのではと思います.それより詳しいことが知りたければ,PRMLでも読めば良いのではないでしょうか(適当*1).

まずはデータを眺める

使用したデータは,caretパッケージのcarsパッケージです*2.中古車販売のデータっぽくて,価格と,走行距離とか気筒とかドア数とかの車に関するカラムが並んでます.データを読み込んで,可視化して,とりあえず lm() してみます.

> library(glmnet)
> library(caret)
> library(psych)
> 
> # load data
> data(cars)
> t(head(cars))
                   1        2        3        4        5        6
Price       22661.05 21725.01 29142.71 30731.94 33358.77 30315.17
Mileage     20105.00 13457.00 31655.00 22479.00 17590.00 23635.00
Cylinder        6.00     6.00     4.00     4.00     4.00     4.00
Doors           4.00     2.00     2.00     2.00     2.00     2.00
Cruise          1.00     1.00     1.00     1.00     1.00     1.00
Sound           0.00     1.00     1.00     0.00     1.00     0.00
Leather         0.00     0.00     1.00     0.00     1.00     0.00
Buick           1.00     0.00     0.00     0.00     0.00     0.00
Cadillac        0.00     0.00     0.00     0.00     0.00     0.00
Chevy           0.00     1.00     0.00     0.00     0.00     0.00
Pontiac         0.00     0.00     0.00     0.00     0.00     0.00
Saab            0.00     0.00     1.00     1.00     1.00     1.00
Saturn          0.00     0.00     0.00     0.00     0.00     0.00
convertible     0.00     0.00     1.00     1.00     1.00     1.00
coupe           0.00     1.00     0.00     0.00     0.00     0.00
hatchback       0.00     0.00     0.00     0.00     0.00     0.00
sedan           1.00     0.00     0.00     0.00     0.00     0.00
wagon           0.00     0.00     0.00     0.00     0.00     0.00
> pairs.panels(cars)
> 
> # lm
> fit.lm <- glm(Price ~ ., data = cars)
> summary(fit.lm)

Call:
glm(formula = Price ~ ., data = cars)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-9513.5  -1540.9    125.4   1470.3  13619.7  

Coefficients: (3 not defined because of singularities)
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.124e+03  9.926e+02  -1.133  0.25773    
Mileage     -1.842e-01  1.256e-02 -14.664  < 2e-16 ***
Cylinder     3.659e+03  1.133e+02  32.286  < 2e-16 ***
Doors        1.567e+03  2.589e+02   6.052  2.2e-09 ***
Cruise       3.409e+02  2.960e+02   1.152  0.24978    
Sound        4.409e+02  2.345e+02   1.880  0.06043 .  
Leather      7.908e+02  2.497e+02   3.167  0.00160 ** 
Buick        9.477e+02  5.525e+02   1.715  0.08670 .  
Cadillac     1.336e+04  6.248e+02  21.386  < 2e-16 ***
Chevy       -5.492e+02  4.397e+02  -1.249  0.21203    
Pontiac     -1.400e+03  4.868e+02  -2.875  0.00414 ** 
Saab         1.228e+04  5.546e+02  22.139  < 2e-16 ***
Saturn              NA         NA      NA       NA    
convertible  1.102e+04  5.419e+02  20.340  < 2e-16 ***
coupe               NA         NA      NA       NA    
hatchback   -6.362e+03  6.104e+02 -10.422  < 2e-16 ***
sedan       -4.449e+03  4.463e+02  -9.969  < 2e-16 ***
wagon               NA         NA      NA       NA    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 8445957)

    Null deviance: 7.8461e+10  on 803  degrees of freedom
Residual deviance: 6.6639e+09  on 789  degrees of freedom
AIC: 15122

Number of Fisher Scoring iterations: 2

http://f.st-hatena.com/images/fotolife/S/SAM/20160221/20160221162514_original.png

相関は各変数そこそこある感じで,極端に高いものはないみたいです.でもlm()の結果をみると,Saturnとcoupeとwagonは落ちちゃってますね.まぁその辺りは気にせず次へ.

glmnetでL1正則化

Rで正則化をするパッケージとして {glmnet} を使います.読み込みにfactor型は使えず,すべて数値型で渡す必要があります.かつYとXは別のmatrixとして渡してあげる必要があるとのこと.実行して,plotしてみます.式中の alpha が1ならL1正則化で,0ならL2正則化になります.そしてその間ならElastic Netです.

> pairs.panels(cars)
> # lasso
> fit.glmnet.lasso <- glmnet(as.matrix(cars[, -1]),
+                            as.matrix(cars[, 1]),
+                            alpha = 1)
> plot(fit.glmnet.lasso)

http://f.st-hatena.com/images/fotolife/S/SAM/20160221/20160221162509_original.png

この図は,左にいくほどL1正則化が強く効いている状態とのことです.さらに cv.glmnet() で5 foldのクロスバリデーションを行います.

> fit.glmnet.lasso.cv <- cv.glmnet(as.matrix(cars[, -1]),
+                                  as.matrix(cars[, 1]),
+                                  nfold = 5,
+                                  alpha = 1)
> plot(fit.glmnet.lasso.cv)
> fit.glmnet.lasso.cv$lambda.min
[1] 18.54925
> fit.glmnet.lasso.cv$lambda.1se
[1] 275.4505
> coef(fit.glmnet.lasso.cv, s = fit.glmnet.lasso.cv$lambda.min)
18 x 1 sparse Matrix of class "dgCMatrix"
                        1
(Intercept)  3332.5402054
Mileage        -0.1816504
Cylinder     3634.6147087
Doors        -628.9022241
Cruise        340.6710317
Sound         382.0774957
Leather       753.7844280
Buick         924.8954890
Cadillac    13401.5393074
Chevy        -482.8713605
Pontiac     -1294.4163546
Saab        12273.6389734
Saturn          .        
convertible 11016.0739756
coupe           .        
hatchback   -1881.6575830
sedan           .        
wagon        4317.2107830

http://f.st-hatena.com/images/fotolife/S/SAM/20160221/20160221162510_original.png

MSE*3が最小になる lambda の値が18.55です.上の図において,log(lambda) = 3 あたりに点線が縦に引かれていますが,これがminの位置になります.で,真ん中ぐらいにある点線が,lambdaの1se地点の値です.lm() の結果と,lambda.minにおける各係数のあたいをみてみると,微妙に違いがありますね.lassoの場合はwagonではなくsedanが係数0に収束しているようです.

ちなみに,coef() を実行したときにオプションで指定している s というのは,coef を表示する際の lambda の値を指します. cranのドキュメントには,p18の predict.cv.glmnet に

Value(s) of the penalty parameter lambda at which predictions are required. Default
is the value s="lambda.1se" stored on the CV object. Alternatively
s="lambda.min" can be used. If s is numeric, it is taken as the value(s) of
lambda to be used.

と記述されています.

L2正則化

同じように,今度は alpha を0にして,ridge回帰をおこなってみます.

> # ridge
> fit.glmnet.ridge <- glmnet(as.matrix(cars[, -1]),
+                            as.matrix(cars[, 1]),
+                            alpha = 0)
> plot(fit.glmnet.ridge)
> ## cv
> fit.glmnet.ridge.cv <- cv.glmnet(as.matrix(cars[, -1]),
+                                  as.matrix(cars[, 1]),
+                                  nfold = 5,
+                                  alpha = 0)
> plot(fit.glmnet.ridge.cv)
> fit.glmnet.ridge.cv$lambda.min
[1] 714.8007
> fit.glmnet.ridge.cv$lambda.1se
[1] 1370.924
> coef(fit.glmnet.ridge.cv, s = fit.glmnet.ridge.cv$lambda.min)
18 x 1 sparse Matrix of class "dgCMatrix"
                        1
(Intercept) 10337.0169783
Mileage        -0.1712252
Cylinder     3170.1294044
Doors        -831.5817913
Cruise       1026.4850774
Sound         224.6847147
Leather       956.7665236
Buick       -1936.7129470
Cadillac    10320.3212670
Chevy       -3629.5572140
Pontiac     -4111.5436304
Saab         7869.4887792
Saturn      -3239.5070500
convertible  9384.0693961
coupe       -1729.7150676
hatchback   -2962.1017215
sedan       -1250.9639471
wagon        2760.2080200

http://f.st-hatena.com/images/fotolife/S/SAM/20160221/20160221162511_original.png

http://f.st-hatena.com/images/fotolife/S/SAM/20160221/20160221162512_original.png

lassoのときとは結果が大きく変わりました.

ElasticNetによる正則化

最後に,ElasticNetを試してみます.基本的には,0 < alpha < 1 であればよいのですが,この alpha の値を決める方法は,{glmnet} 内では提供されていません*4ので,今回はえいやで alpha=0.5 にしちゃいます.ちゃんとやりたいのであれば,{caret} あたりを使って,lambda とalpha を共に変動させてみるのがよいかと思います.

> fit.glmnet.elasticnet.cv <- cv.glmnet(as.matrix(cars[, -1]),
+                                       as.matrix(cars[, 1]),
+                                       nfold = 5,
+                                       alpha = 0.5)
> plot(fit.glmnet.elasticnet.cv)
> fit.glmnet.elasticnet.cv$lambda.min
[1] 33.80277
> fit.glmnet.elasticnet.cv$lambda.1se
[1] 416.7364
> coef(fit.glmnet.elasticnet.cv, s = fit.glmnet.elasticnet.cv$lambda.min)
18 x 1 sparse Matrix of class "dgCMatrix"
                        1
(Intercept)  3864.0228103
Mileage        -0.1815926
Cylinder     3625.4650658
Doors        -626.2371069
Cruise        371.2344641
Sound         386.8548536
Leather       767.8711841
Buick         402.2514994
Cadillac    12865.1630823
Chevy       -1013.4253903
Pontiac     -1821.5463454
Saab        11707.2380808
Saturn       -464.8502014
convertible 11028.5218995
coupe           .        
hatchback   -1877.0655854
sedan          -3.3327249
wagon        4323.9610872

ということで,{glmnet} でL1/L2正則化を試してみました.なお,コードはgistにあげてあるので,興味がある方はどうぞ.

*1:単に「正則化 amazon」で検索したら一番上に出てきたのがビショップ本だっただけです.

*2:MiluHatsuneさんのエントリがわかりやすかったので,それをベースになぞってる感じです.

*3:Mean Squared Error,つまり予測値と実際の値の差の二乗和を意味します.

*4:cross validation - Choosing optimal alpha in elastic net logistic regression - Cross Validated