LASSO and Ridge regression

今回はLASSOとリッジ回帰についてです。

パッケージは「glmnet」、「lars」、「lasso2」で実行できます。

glmnetとlarsの作者はFriedman、Hastie、Efron、Tibshiraniと有名な先生ですが、lasso2の作者は知らないです。。

内容もほぼ一緒なので、LASSOをするときはglmnet一択で良いと思います。



まずは使用例から。。。


データはLARSパッケージにあるdiabetesを使います。

このデータである結果変数y(中性脂肪?)をx(性別や血圧など)によって予測するモデルを作ります。

まずは単純な線形回帰をします。


library(lars)
library(glmnet)

data(diabetes)
Linear     <- lm(diabetes$y ~ diabetes$x)
Linear$coefficients

これが推定結果です。

  (Intercept) diabetes$xage diabetes$xsex diabetes$xbmi diabetes$xmap  diabetes$xtc diabetes$xldl diabetes$xhdl 
    152.13348     -10.01220    -239.81909     519.83979     324.39043    -792.18416     476.74584     101.04457 
diabetes$xtch diabetes$xltg diabetes$xglu 
    177.06418     751.27932      67.62539 


この線形回帰をLASSOで推定します。

ここで注意しなければいけないのが、xはすべて標準化しなくてはなりません(オプションデフォルトでちゃんと標準化されるようになってます)。

LASSO.net  <- glmnet(diabetes$x, diabetes$y, alpha = 1)
plot(LASSO.net)

f:id:isseing333:20110309203644p:image

LASSOは難しそうな手法に見えますが、実はそんなことありません。

古典的な回帰手法を、「ペナルティを付けて」推定しているだけです。

この結果のグラフで、横軸がペナルティの強さ、縦軸が回帰係数の推定結果を示しています。

それぞれの折れ線が説明変数xを示してまして、例えば2本の緑の線はBMIとLTGです。

ペナルティは左が強く、右に行くに従って弱くなっています。



まずは直感的な説明ですが、最もペナルティが強い左端ではどの変数の回帰係数も0になってます。

そして最もペナルティの弱い右端では変数の回帰係数が単純な線形回帰と一致するのです。

つまり最初強くしていたペナルティを弱めていく事によって、結果変数と関連の強い変数から推定していき、最終的には線形回帰と等しくなるのです。




次にLASSOがどのようにペナルティをかけているのかを数式で説明します。

参考文献はこちら→http://www.stanford.edu/~hastie/Papers/glmnet.pdf

"Regularization Paths for Generalized Linear Models via Coordinate Descent"




LASSOの数式部分を抜粋します。



f:id:isseing333:20110309210911p:image



数式(1)の1の部分だけだと通常の線形回帰(最小二乗法)で、LASSOには2のペナルティ項が入っています。

ペナルティの中身は数式(3)になっていて、α=1だとLASSO、α=0だとリッジ回帰、その他だとElastic Netになります。

このペナルティ項はλの大きさで変わります。


λが大きいとペナルティがきつく、λが小さいとペナルティは緩いです。

λが0だと最小二乗法と一致します。


さっきの図で横軸に示されていたL1 normは、(3)の式の値です。

これは推定値βの絶対値を足したものです。




さて、このペナルティがどうなっているのかを図で表現してみます。

乱数でデータを発生させて、自作関数で線形回帰とLASSOを作ってみます。


正しいモデルは「y = 2x + 誤差」で、x1はyとは無相関のダミーです。

まずは通常の線形回帰。

#---simu data
set.seed(1)
x   <- rnorm(1000)
set.seed(2)
x1 <- rnorm(1000)
y   <- 2 * x + rnorm(1000)

plot(x, y)


LinearRegression <- function(beta){
	(y - beta[1] - beta[2]*x)%*%(y - beta[1] - beta[2]*x)
}
optim(c(0, 0), LinearRegression)$par
lm(y ~ x)

自作関数とlm関数の結果が微妙に違いますが、許容範囲ということで。。。



次に、これを利用してLASSOを作ります。

LassoRegression <- function(beta){
	LassoRes <- sum(abs(beta))
	RSS <- (y - beta[1] - beta[2]*x - beta[3]*x1)%*%
				(y - beta[1] - beta[2]*x - beta[3]*x1)/(2*length(x))
	RSS + lambda * LassoRes
}

lambda <- 1
optim(c(0, 0, 0), LassoRegression)$par

線形回帰にペナルティを入れてます。



この関数を使って、以下のように回帰残差、ペナルティ項、目的関数と推定値のプロットを描いてみます。

(1)の数式では、1が回帰残差、2がペナルティ項、1と2を足したものが目的関数です。

#---relationship RSS and beta
BetaX      <- seq(-5, 10, by = 0.1)
BetaInt    <- mean(y) - BetaX*mean(x)
lambda   <- 1
LassoRes <- lambda * (abs(BetaInt) + abs(BetaX))


RSS       <- rep(NA, length(BetaX))
for(i in 1:length(BetaX)){
	RSS[i]    <- (y - BetaInt[i] - BetaX[i]*x)%*%(y - BetaInt[i] - BetaX[i]*x)/(2*length(x))
}
LassoBeta <- BetaX[(RSS + LassoRes) == min(RSS + LassoRes)]

plot(BetaX, RSS, type = "n", ylab = "objective function", main = paste("Lambda = ", lambda, sep = ""))
lines(BetaX, RSS)
lines(BetaX, LassoRes, col = "blue")
lines(BetaX, RSS + LassoRes, col = "red")

abline(v = LassoBeta, lty = 2)
text(LassoBeta + 0.1, 30, paste("Beta = ", LassoBeta, sep = ""), adj = 0)
legend("bottomright", c("RSS", "LassoPanalty", "ObjectiveFunction"), col = c("black", "blue", "red"), lty = 1)


λを3、2、1、0と変えて図を描いてみます。

黒のRSSが回帰残差、青がペナルティ項、赤が2つを足した目的関数です。

この目的関数が最小になる点が、回帰係数(β)の推定値になります。

f:id:isseing333:20110309212634p:image

f:id:isseing333:20110309212635p:image

f:id:isseing333:20110309212636p:image

f:id:isseing333:20110309212637p:image


λが小さくなるとペナルティ項の傾きが緩くなっています。

それによって目的関数の形が変わり、最終的には回帰残差だけが残ります。

回帰残差だけということは最小二乗法と等しいということです。

また、回帰係数の推定値がが0から2へ大きくなっていることも確認できます。



LASSOではペナルティ項が折れ線になっていますが、リッジ回帰だと2次関数になるわけです。

つまりLASSOもリッジ回帰も「λを変化させて回帰係数の推定値を変化させている」のです。

実用的には、cv.glmnetでクロスバリデーションをしながら最適なλを求めることになります。



そこで求まった回帰モデルは、単純な線形モデルです。

LASSOはover-fittingしないようにパラメータ推定をする手順であり、結果自体が複雑なわけではないのです。



さらにfamilyを指定することで、ロジスティック回帰や生存時間解析にもLASSOを応用できます。

単純に回帰をするより、LASSOを使う方が良さそうですね。

ページTOPへ