医療から宇宙へトップ分析会社

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)

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の数式部分を抜粋します。

数式(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)