Lasso 回帰は、データに多重共線性が存在するときに、回帰モデルを適合させるために使用できる手法である。
簡単に言うと、最小2乗回帰は，残差2乗和 (RSS) を最小化する係数推定を見つける。
RSS = Σ(yi - ŷi)2
ここで
Σ : 和を意味するギリシャ記号 yi: 番目の観測の実際の応答値 ŷi: 重回帰モデルに基づく予測応答値
逆にlasso 回帰は、次を最小化しようとする：
RSS + λΣ|βj|．
ここで j は1 から p の予測変数で、λ ≥ 0 である。
式中のこの第2項は、縮小制約として知られている。このチュートリアルは、Rでlasso 回帰を実行する方法のステップ・バイ・ステップの例を提供する。
ステップ1：データのロード
この例では、mtcarsというR組み込みデータセットを使用する。hpを応答変数として、以下の変数を予測変数として使用する： mpg
- hp : Gross horsepower：馬力
- mpg : Miles/(US) gallon：燃費
- wt : weight(1000 lbs)（重量）
- drat : Rear axle ratio（ドライブシャフトの回転数と車軸の回転数との比率）
- qsec : 1/4 mile time（停止状態から約400m進むまでの時間）
lasso 回帰を実行するために、glmnet パッケージの関数を使うことにする。
#応答変数の定義
y <- mtcars$hp
#予測変数の行列を定義
x <- data.matrix(mtcars[, c('mpg', 'wt', 'drat', 'qsec')])
ステップ 2: Lasso回帰モデルのフィット 次に、glmnet()関数を使用してLasso回帰モデルをフィットし、alpha=1を指定する。
alphaを0に設定することは、リッジ回帰を使用することと同じであり、alphaを0と1の間の値に設定することは、elastic netを使用することと同じであることに注意する。
λにどの値を使うかを決定するために、k-fold 交差検証（k-fold cross-validation）を実行し、テストの平均2乗誤差（MSE）が最小になるλ値を特定する。
関数cv.glmnet()は、k = 10 foldsを用いて自動的にk-fold交差検証を行うことに注意する。
library(glmnet) #perform k-fold cross-validation to find optimal lambda value cv_model <- cv.glmnet(x, y, alpha = 1) #find optimal lambda value that minimizes test MSE best_lambda <- cv_model$lambda.min best_lambda
[1] 2.431182
#produce plot of test MSE by lambda value plot(cv_model)
テストのMSEを最小化するラムダ値は5.616345であることが判明した。
ステップ3：最終モデルの分析
最後に、最適なラムダ値によって生成された最終モデルを分析する。
次のコードを使って、このモデルの係数推定値を得ることができる：
#最適モデルの係数を求める best_model <- glmnet(x, y, alpha = 1, lambda = best_lambda) coef(best_model)
5 x 1 sparse Matrix of class "dgCMatrix"
s0
(Intercept) 485.152675
mpg -2.936266
wt 21.698919
drat .
qsec -19.569135
lasso 回帰が係数をゼロまで縮小したので、予測変数dratの係数は表示されていません。これがリッジ回帰とlasso 回帰の重要な違いであることに注意。リッジ回帰は、すべての係数をゼロに向かって縮めますが、lasso 回帰は、係数を完全にゼロに縮めることによって、モデルから予測変数を除去する可能性があります。
我々はまた、新しいオブザベーションで予測を行うために、最終的なlasso 回帰モデルを使用できる。
たとえば、次の属性を持つ新しい車があるとしよう：
- mpg: 24
- wt: 2.5
- drat: 3.5
- qsec: 18.5
次のコードは、この新しい観測のhpの値を予測するために、フィットしたlasso回帰モデルを使用する方法を示している：
#define 新しい観測 new = matrix(c(24, 2.5, 3.5, 18.5), nrow=1, ncol=4) #応答値を予測するためにlasso回帰モデルを使用する predict(best_model, s = best_lambda, newx = new)
[1,] 106.9006
入力値に基づき、モデルはこの車のhp値を106.9006と予測する。
最後に、トレーニングデータに対するモデルのR2乗を計算する：R2乗を計算すると、以下のようになる：
#フィットしたベストモデルを使って予測を行う y_predicted <- predict(best_model, s = best_lambda, newx = x) #SSTとSSEを求める sst <- sum((y - mean(y))^2) sse <- sum((y_predicted - y)^2) #R二乗を求める rsq <- 1 - sse/sst rsq
[1] 0.8043193
R2乗は0.8043193となった。つまり、最良のモデルは、トレーニング・データの応答値の変動の80.47%を説明することができた。
この例で使用したRコード一式は、ここで見ることができる。
