以前にStataでLasso、リッジ回帰の使い方は解説している。
今回はRでの実行例について。
Lasso Regression Model with R code
Tibshirani (1996) は、パラメータの選択と縮小のために、LASSO (Least Absolute Shrinkage and Selection Operator) モデルを導入している。リッジ・モデルは、重要でない変数の係数をゼロに近づけるが、正確にはゼロにしないので、縮小の点ではこれに似ているが、選択機能はない。
これらの回帰モデルは正則化回帰モデルまたはペナルティ付き回帰モデルと呼ばれる。特にLassoは非常に強力で、変数の数が100以上、1000以上、10000以上、...というような大きなデータセットにも対応できる。従来の線形回帰モデルでは、このようなビッグデータには対応できない。
ここでは、説明のためにYを標準化しているが、これは必ずしなければならないものではない。この制約項は、ハイパーパラメータλによって調整される。 ハイパーパラメータは、手動検索やクロスバリデーションのプロセスを通じて、ユーザが外生的に与える。
ある変数がLassoに含まれるが、RSS(Residual Sum of Squares) の減少が無視できるほど小さい場合（すなわち、RSSを0.000000001減少させる場合）、縮小ペナルティの影響は大きくなる。これは、この変数の係数がゼロ（Lasso）またはゼロに近い（Ridge）ことを意味する。 Ridge（convex and differentiable）と異なり、Lasso（nononvex and nondifferentiable）はほとんどの問題で閉形式解を持たないので、循環座標降下アルゴリズムを使用する。
データ作成
例えば、ある時系列データを作成し、Xを10個のランダムな時系列（変数）とし、Yは、あらかじめ決められた係数とランダムに描かれた誤差項を持つ変数である。標準的な線形回帰、Lasso回帰、Ridge回帰の違いを明確に理解するために、いくつかの係数はゼロに設定されている。
library(glmnet)
graphics.off() # 全てのグラフを消去する。
rm(list = ls()) # ワークスペースからすべてのファイルを削除する
N = 500 # 観測数
p = 20 # 変数の数
X variable
X = matrix(rnorm(N*p), ncol=p)
# 標準化前
colMeans(X) # 平均値
apply(X,2,sd) # 標準偏差
# スケール : 平均値 = 0, 標準偏差 = 1
X = scale(X)
# 標準化後
colMeans(X) # 平均値
apply(X,2,sd) # 標準偏差
Y variable
beta = c( 0.15, -0.33, 0.25, -0.25, 0.05,rep(0, p/2-5),
-0.25, 0.12, -0.125, rep(0, p/2-3))
# Y 変数、標準化 Y
y = X%*%beta + rnorm(N, sd=0.5)
y = scale(y)
Model
lambda <- 0.01
# 切片(-1)なしの標準線形回帰
li.eq <- lm(y ~ X-1)
# lasso
la.eq <- glmnet(X, y, lambda=lambda,
family="gaussian",
intercept = F, alpha=1)
# Ridge
ri.eq <- glmnet(X, y, lambda=lambda,
family="gaussian",
intercept = F, alpha=0)
Results (lambda=0.01)
df.comp <- data.frame(
beta = beta,
Linear = li.eq$coefficients,
Lasso = la.eq$beta[,1],
Ridge = ri.eq$beta[,1]
)
df.comp
結果。
beta Linear Lasso Ridge X1 0.150 0.203906008 0.196840145 0.202303868 X2 -0.330 -0.423689566 -0.415610815 -0.419470196 X3 0.250 0.327343960 0.318629609 0.324341112 X4 -0.250 -0.277376345 -0.267598251 -0.274612904 X5 0.050 0.016659147 0.008427179 0.016633141 X6 0.000 -0.056507432 -0.046301219 -0.056430887 X7 0.000 -0.007049888 0.000000000 -0.007023013 X8 0.000 -0.034716189 -0.026274040 -0.034599799 X9 0.000 0.003966152 0.000000000 0.003801209 X10 0.000 0.034233995 0.022210885 0.033949470 X11 -0.250 -0.347025589 -0.337845528 -0.343731067 X12 0.120 0.161268813 0.148441700 0.159534076 X13 -0.125 -0.138134657 -0.130234065 -0.136828494 X14 0.000 0.009357876 0.000000000 0.009335303 X15 0.000 -0.029461809 -0.019374633 -0.029580248 X16 0.000 -0.003676459 0.000000000 -0.003517786 X17 0.000 0.016710522 0.005757849 0.016285006 X18 0.000 0.036055451 0.024315152 0.035556900 X19 0.000 -0.071789762 -0.063566916 -0.071340301 X20 0.000 -0.029371234 -0.020993406 -0.029250390
Results (lambda=0.1)
lambda <- 0.1
# lasso
la.eq <- glmnet(X, y, lambda=lambda,
family="gaussian",
intercept = F, alpha=1)
# Ridge
ri.eq <- glmnet(X, y, lambda=lambda,
family="gaussian",
intercept = F, alpha=0)
df.comp <- data.frame(
beta = beta,
Linear = li.eq$coefficients,
Lasso = la.eq$beta[,1],
Ridge = ri.eq$beta[,1]
)
df.comp
係数の収縮 ( λ入力がある場合、またはλ入力がない場合 )
# lasso
la.eq <- glmnet(X, y, family="gaussian",
intercept = F, alpha=1)
# Ridge
ri.eq <- glmnet(X, y, family="gaussian",
intercept = F, alpha=0)
# plot
x11(); par(mfrow=c(2,1))
x11(); matplot(log(la.eq$lambda), t(la.eq$beta),
type="l", main="Lasso", lwd=2)
x11(); matplot(log(ri.eq$lambda), t(ri.eq$beta),
type="l", main="Ridge", lwd=2)
クロスバリデーションの実行とラムダの選択
mod_cv <- cv.glmnet(x=X, y=y, family="gaussian",
intercept = F, alpha=1)
# plot(log(mod_cv$lambda), mod_cv$cvm)
# cvm : クロスバリデーションの平均誤差
# - 長さlength(lambda)のベクトル
# lambda.min : MSEが最小となるλを指定する。
# lambda.min : MSEが最小となるλ。
# lambda.1se : MSEが1標準以内に収まる最大のλ。
# MSEが1標準誤差以内に収まる最大のλ
# 最小の MSE の # 標準誤差の範囲内となる最大の λ。
x11(); plot(mod_cv)
coef(mod_cv, c(mod_cv$lambda.min,
mod_cv$lambda.1se))
print(paste(mod_cv$lambda.min,
log(mod_cv$lambda.min)))
print(paste(mod_cv$lambda.1se,
log(mod_cv$lambda.1se)))
推計結果
### Results (lambda=0.1)
beta Linear Lasso Ridge
X1 0.150 0.203906008 0 2.532565e-37
X2 -0.330 -0.423689566 0 -4.294546e-37
X3 0.250 0.327343960 0 3.634778e-37
X4 -0.250 -0.277376345 0 -2.881883e-37
X5 0.050 0.016659147 0 3.116602e-38
X6 0.000 -0.056507432 0 -1.045405e-37
X7 0.000 -0.007049888 0 -1.477689e-38
X8 0.000 -0.034716189 0 -6.645734e-38
X9 0.000 0.003966152 0 -9.685793e-39
X10 0.000 0.034233995 0 4.278838e-38
X11 -0.250 -0.347025589 0 -3.759453e-37
X12 0.120 0.161268813 0 1.410190e-37
X13 -0.125 -0.138134657 0 -1.505259e-37
X14 0.000 0.009357876 0 1.574418e-38
X15 0.000 -0.029461809 0 -7.327030e-38
X16 0.000 -0.003676459 0 1.173640e-38
X17 0.000 0.016710522 0 -1.187330e-38
X18 0.000 0.036055451 0 2.302968e-38
X19 0.000 -0.071789762 0 -9.991021e-38
X20 0.000 -0.029371234 0 -5.569975e-38
## Results (lambda=0.1)
beta Linear Lasso Ridge
X1 0.150 0.203906008 0.12054182 0.188895963
X2 -0.330 -0.423689566 -0.33066998 -0.385000328
X3 0.250 0.327343960 0.24220099 0.299635890
X4 -0.250 -0.277376345 -0.18160109 -0.252083267
X5 0.050 0.016659147 0.00000000 0.016305312
X6 0.000 -0.056507432 0.00000000 -0.055408297
X7 0.000 -0.007049888 0.00000000 -0.006795353
X8 0.000 -0.034716189 0.00000000 -0.033520493
X9 0.000 0.003966152 0.00000000 0.002549903
X10 0.000 0.034233995 0.00000000 0.031595936
X11 -0.250 -0.347025589 -0.25435556 -0.316732434
X12 0.120 0.161268813 0.04687196 0.145391765
X13 -0.125 -0.138134657 -0.04723004 -0.126120514
X14 0.000 0.009357876 0.00000000 0.009086251
X15 0.000 -0.029461809 0.00000000 -0.030239244
X16 0.000 -0.003676459 0.00000000 -0.002300025
X17 0.000 0.016710522 0.00000000 0.012992814
X18 0.000 0.036055451 0.00000000 0.031594759
X19 0.000 -0.071789762 0.00000000 -0.067468104
X20 0.000 -0.029371234 0.00000000 -0.028188778
推定結果は、データ生成過程の不確実性にもかかわらず、モデル間で同様の結果を示している。特に、Lassoは重要でない、あるいは重要でない変数を係数ゼロとして識別している。変数選択と縮小効果はλに比例して強く、縮小経路である制約パラメータ（log(λ)）の変化に対する推定係数の変化を下図に示す。
モデル選択
Lassoで最も重要なのは、最適なλを選択することである。 これは、クロスバリデーションの過程で決定される。glmnetのcv.glmnet()関数は、λの範囲を適切に設定してクロスバリデーションの結果を提供する。 . この出力を使って、log(λ)とMSE(measuren squared error)のグラフを描くことができる。
上図から、1つ目の候補はMSEが最小となるλだが、このモデルは変数が多いと思われる。これは、やや発見的、経験的なアプローチだが、変数の数を減らすという点ではメリットがある。しかし、縮み過程のパターンを見つけるには、目視検査が非常に重要なツールになる。
以下の結果は、MSE最小化λとMSE最小化1seλの下での推定係数を報告する。をそれぞれ算出する。
21 x 2 sparse Matrix of class "dgCMatrix"
s1 s2
(Intercept) . .
V1 0.2006863980 0.159138784
V2 -0.4015852307 -0.365572144
V3 0.2446988172 0.204828449
V4 -0.3335434947 -0.284715097
V5 . .
V6 -0.0318423289 .
V7 -0.0077604995 .
V8 . .
V9 0.0073594620 .
V10 0.0001035328 .
V11 -0.3393509455 -0.300525072
V12 0.0845500490 0.047706234
V13 -0.1972962269 -0.157268107
V14 . .
V15 -0.0333920845 .
V16 0.0138355335 .
V17 -0.0413405995 -0.003093179
V18 -0.0036850700 .
V19 -0.0379703195 .
V20 . .
> print(paste(mod_cv$lambda.min,
+ log(mod_cv$lambda.min)))
[1] "0.01686454287377 -4.08254191585989"
> print(paste(mod_cv$lambda.1se,
+ log(mod_cv$lambda.1se)))
[1] "0.0565232466025991 -2.87310328115595"
予測
Lasso回帰のパラメータを推定した後、このモデルを用いて予測を行う必要がある。予測演習では、制約項ではなく、推定された係数を使用します。予測方法を見ると、Lasso回帰、Ridge回帰、線形回帰の3つのモデルは、ペナルティ項が推定に使われるだけなので、同じである。
この投稿をもとに、符号制限付きLassoモデルについて説明する。経済理論や実証的な様式化された事実が特定の符号を推奨しているため、係数の符号に制約を設けることは重要である。
Tibshirani, Robert (1996). " Regression Shrinkage and Selection via the lasso," Journal of the Royal Statistical Society 58-1, 267-88.
