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

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

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

```math
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正則化項を加えた回帰モデルです。
正則化項を加えることにより、いくつかの係数の値が０となり、自動的に変数選択を行うことができます。
また、相関が高い変数群がある場合、一つのみ変数を選択します。

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

詳しくは、tekenukoさんの[記事](http://tekenuko.hatenablog.com/entry/2017/11/18/214317)等をご覧ください。

パッケージ 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")
```

　<img src="https://qiita-image-store.s3.amazonaws.com/0/197919/2bf94603-b636-d348-9abf-e277f7a19f24.png" width=70%　>
縦軸がかく係数、横軸がλの対数の絶対値、上部の数字が各λでのモデルにおける係数が非０の変数の数です。
λの値を選択した値から少し大きくすると変数を５つのみ選択してくれるようです。

選択したモデルによる予測値の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は多次元小標本の場合に力を発揮するんですかね。
もう少し頑張って調べてみます。