Help us understand the problem. What is going on with this article?

stepwiseとlasso回帰における変数選択を比較

More than 1 year has passed since last update.

はじめに

回帰分析において、いくつか変数選択の手法は存在します。
その中で、stepwise法とlasso回帰のどちらが適切な変数選択ができるのか、なんとなく気になりました。
検索したところ意見が別れている模様...(英語が苦手な私が頑張って英語のページを見た結果)
そのため、サンプルデータを用いて実験してみました。

設定

今回用いた真のモデルは以下のようなモデルです。

Y = 0.68 X_1-0.41 X_2+0.27 X_3-0.89 X_4+0.52 X_5+0.66 X_6+2.84+\varepsilon\\

ノイズεはN(0,1)に従うものとしました。
説明変数は、真のモデルの変数以外にも14種の変数が観測できることにしました。
各20種の変数はN(0,1)に従います。

サンプル数は1000個、そのうち800を学習データ、200個をテストデータとしました。

stepwise

stepwiseは、モデルに用いる説明変数を減少・増化をさせながら、指標が最良となる変数の組み合わせを探索する手法です。
今回は変数減少法を用い、指標はAICとしました。

> model.rg <- glm(Y.train~., data = X.train)
> step.result<-step(model.rg)
      "略"
> step.result

Call:  glm(formula = Y.train ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X13 + 
    X18 + X20, data = X.train)

Coefficients:
(Intercept)           X1           X2           X3           X4           X5           X6  
    2.87439      0.67332     -0.42552      0.29593     -0.90412      0.50839      0.64421  
         X7          X13          X18          X20  
   -0.03680      0.05682      0.02724      0.02971  

Degrees of Freedom: 799 Total (i.e. Null);  789 Residual
Null Deviance:      2188 
Residual Deviance: 186  AIC: 1127

ちなみに、真のモデルの係数と推定された係数の平均二乗誤差(MSE)は3.8916e^-4でした。

真のモデルの説明変数も選択し、推定された係数は真値に近い値を算出しています。
係数の有意確率の検定を行なってみると、真のモデルの説明変数は有意な差を示します。
しかし、それ以外の説明変数でもP<0.05というものが存在しました。

> summary(step.result)

Call:
glm(formula = Y.train ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X13 + 
    X18 + X20, data = X.train)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6503  -0.3270  -0.0165   0.3118   1.4407  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.87439    0.01724 166.757   <2e-16 ***
X1           0.67332    0.01672  40.272   <2e-16 ***
X2          -0.42552    0.01680 -25.324   <2e-16 ***
X3           0.29593    0.01674  17.680   <2e-16 ***
X4          -0.90412    0.01720 -52.576   <2e-16 ***
X5           0.50839    0.01743  29.174   <2e-16 ***
X6           0.64421    0.01825  35.292   <2e-16 ***
X7          -0.03680    0.01716  -2.145   0.0323 *  
X13          0.05682    0.02873   1.978   0.0483 *  
X18          0.02724    0.01725   1.580   0.1146    
X20          0.02971    0.01776   1.673   0.0947 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.2357639)

    Null deviance: 2188.27  on 799  degrees of freedom
Residual deviance:  186.02  on 789  degrees of freedom
AIC: 1127.3

Number of Fisher Scoring iterations: 2

テストデータで出力した予測値の平均二乗誤差は次にようになりました。

> result<-predict(model.rg.step, X.test)
> mse.stepwise<-sum((result-Y.test)^2)/200
> mse.stepwise
[1] 0.2517455

lasso回帰

lasso回帰は、目的関数にL1正則化項を加えた回帰モデルです。
正則化項を加えることにより、いくつかの係数の値が0となり、自動的に変数選択を行うことができます。
また、相関が高い変数群がある場合、一つのみ変数を選択します。

L= ||Y-\sum _{i = 0}^p  X_i \beta_i ||^2  + \lambda \sum _{i = 0}^p \left| \beta _i \right|   

詳しくは、tekenukoさんの記事等をご覧ください。

パッケージ glmnet を利用しました。
はじめに、交差検証を行い最適なパラメータλを得ました。
そして、最適なパラメータを利用したモデルを生成。
選択された変数はstepwiseよりも多いです。
また、推定された係数のいくつかは、真値よりも少し抑えられているようです。

> library("glmnet")
> lasso.model.cv <- cv.glmnet(x = X.train, y = Y.train, family = "gaussian", alpha = 1)
> lasso.model.cv$lambda.min #最適なパラメーター
[1] 0.01249896

> bestLamda<-lasso.model.cv$lambda.min
> lasso.model <- glmnet(x = X.mt, y = Y, family = "gaussian", lambda = bestLamda, alpha = 1)
> lasso.model$beta #係数を表示
20 x 1 sparse Matrix of class "dgCMatrix"
              s0
X1   0.666474030
X2  -0.415939795
X3   0.288122593
X4  -0.896548564
X5   0.499017901
X6   0.634270222
X7  -0.025958793
X8   .          
X9  -0.007309271
X10  0.003286149
X11  .          
X12  .          
X13  0.039858130
X14  0.010436849
X15  .          
X16 -0.004412281
X17  0.001093284
X18  0.018225320
X19  .          
X20  0.021279948  

> lasso.model$a0 #切片
      s0 
2.872088 

真のモデルの係数との平均二乗誤差(MSE)は3.8274e^-4でした。

ggfortifyパッケージを利用してpath図を作成してみます。
glmnet関数でlambdaを指定せずにモデルを生成し、プロット。

> library(ggfortify)
> lasso.allmodel <- glmnet(x = X.mt, y = Y, family = "gaussian", alpha = 1)
> autoplot(lasso.allmodel, xvar = "lambda")

 
縦軸がかく係数、横軸がλの対数の絶対値、上部の数字が各λでのモデルにおける係数が非0の変数の数です。
λの値を選択した値から少し大きくすると変数を5つのみ選択してくれるようです。

選択したモデルによる予測値のMSEは、stepwiseとほぼ一緒ですね。

> result<-predict(lasso.model.cv,newx = X.test,s=bestLamda)
> mse.lassso<-sum((result-Y.test)^2)/200
> mse.lasso
[1] 0.2543211

少しまとめ

stepwiseとlassoの変数選択を比較してみました。

stepwise lasso回帰
選択された変数 10 15
係数のMSE 3.8916e^-4 3.8274e^-4
予測のMSE 0.2517455 0.2543211

stepwiseの方が、より良い変数選択ができているようです。
真のモデルが、より複雑なモデルに対しては結果が違うのでしょうか。
また、lassoは多次元小標本の場合に力を発揮するんですかね。
もう少し頑張って調べてみます。

saltcooky
<今の主はHatenaです> 備忘録ぽいものを自動生成するGANです. 仕事でデータ分析していますが, 知識と経験はありません.
https://saltcooky.hatenablog.com/
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away