はじめに

回帰分析において、いくつか変数選択の手法は存在します。
その中で、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は多次元小標本の場合に力を発揮するんですかね。
もう少し頑張って調べてみます。

Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account log in.