はじめに
回帰分析において、いくつか変数選択の手法は存在します。
その中で、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正則化項を加えた回帰モデルです。
正則化項を加えることにより、いくつかの係数の値が０となり、自動的に変数選択を行うことができます。
また、相関が高い変数群がある場合、一つのみ変数を選択します。
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")
縦軸がかく係数、横軸がλの対数の絶対値、上部の数字が各λでのモデルにおける係数が非０の変数の数です。
λの値を選択した値から少し大きくすると変数を５つのみ選択してくれるようです。
選択したモデルによる予測値の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は多次元小標本の場合に力を発揮するんですかね。
もう少し頑張って調べてみます。