【R言語】リッジ回帰・ラッソ回帰　glmnetを用いたL1, L2正則化

R言語でリッジ回帰・ラッソ回帰を実行する関数やその実行例を紹介します。

この記事では、モデルの過適合を回避するための正則化（regularization）をRで行い、実際にテストデータに対する予測精度が良くなっているか確認します。

この記事で紹介するRファイルは以下からダウンロードできます。

R言語　L1, L2正則化

1 ファイル 2.09 KB

ダウンロード

通常の重回帰分析や一般化線形モデル、多項式回帰については以下の記事を参照。

: 【R言語】関数lmによる線形回帰　単回帰分析・重回帰分析　回帰分析で予測する
R言語の回帰分析の実行方法を紹介していきます。回帰分析の実行方法だけでなく、回帰分析により作ったモデルでテストデータの予測も行っていきます。実行方法では、回帰分析の結果の見方やその抽出方法を詳しく解 ...
続きを見る

: 【R言語】一般化線形モデル　関数glmの使いかた
一般化線形モデルを行う方法を解説していきます。 R言語で回帰分析を行う方法については、以下の記事で紹介しています。一般化線形モデルを用いることで、母分布が正規分布でなくても線形モデルを構築することが ...
続きを見る

正則化

物理学や機械学習などでモデルの過適合や過学習を回避する手法の1つとして正則化があります。

次のように、誤差にペナルティー項を加えることで回帰係数が過剰に大きくならないように最適化を行うことができます。

\begin{align} \| \boldsymbol{y} - f(\boldsymbol{\beta}) \|_2^2+ \lambda\cfrac{1}{p} \| \boldsymbol{\beta} \|_p^p,\end{align}

ここに、$\|\boldsymbol{\beta} \|_p^p$は$\boldsymbol{\beta}$の$p$-ノルム。

$p = 2$ときL2正則化（リッジ回帰）

\begin{align}\label{eq1} \| \boldsymbol{y} - f(\boldsymbol{\beta}) \|_2^2 + \lambda\cfrac{1}{2} \| \boldsymbol{\beta} \|_2^2\tag{1}\end{align}

となり、$p=1$のときにL1正則化（ラッソ回帰）

\begin{align}\label{eq2} \| \boldsymbol{y} - f(\boldsymbol{\beta}) \|_2^2 + \lambda \| \boldsymbol{\beta} \|_1\tag{2}\end{align}

に対応します。

また、機械学習でよく用いられるエラスティックネットは、次のL2正則化とL1正則化のペナルティー項の和からなるモデルで与えられます。

\begin{align}\label{eq3} \| \boldsymbol{y} - f(\boldsymbol{\beta}) \|_2^2 + \lambda\left[\alpha \|\boldsymbol{\beta}\|_1+ \cfrac{(1 - \alpha)}{2} \| \boldsymbol{\beta} \|_2^2\right], \quad 0\leq \alpha \leq 1. \tag{3}\end{align}

$\alpha = 0$のとき\eqref{eq1}のL1正則化、$\alpha = 1$のとき\eqref{eq2}のL2正則化に一致します。

関数glmnet

以下関数glmnetとその引数の説明です。引数が多いのでよく用いるものだけ載せときます。

glmnet(x, y, family=c("gaussian","binomial","multinomial"), alpha = 1, nlambda = 100, lambda.min = ifelse(nobs

関数glmnetの引数
x	説明変数に対応するベクトルまたは行列。
y	目的変数に対応するベクトルまたは行列。
family	character型。目的変数の母集団分布。"gaussian"で正規分布、"binomial"で二項分布、"multinomial"多項分布。
alpha	numeric型。\eqref{eq3}の$\alpha$。
nlambda	numeric型。ハイパーパラメータ$\lambda$の数。
lambda.min	numeric型。ハイパーパラメータ$\lambda$の最小値。

glmnetでリッジ回帰、ラッソ回帰、エラスティックネットを実行するには、次のように引数alphaを与える必要があります。

L2正則化（リッジ回帰）：alpha = 0
L1正則化（ラッソ回帰）：alpha = 1
エラスティックネット：0 < alpha < 1

実行例

上で紹介した関数glmnetを用いて、実際にL1正則化（ラッソ回帰）とL2正則化（リッジ回帰）を行っていきます。

多項式回帰にL2正則化を適用することで過適合を抑えたり、スパ―ス推定の観点でのL1正則化とL2正則化の違いについてみたりしていきます。

リッジ回帰を用いて過適合を回避する（L2正則化）

まずは、L2正則化（リッジ回帰）の実行例について見ていきます。多項式回帰にL2正則化を適用していきます。

データセットとして、次の多項式にホワイトノイズを加えた乱数を使用します。

library(glmnet)
library(ggplot2)
##過適合を回避する例##
set.seed(7)

polynomial <- function(x) {
  return(x * (x + 2) * (x + 1))
}

n <- 30
x <- runif(n, -2, 1)
whiteNoise <- rnorm(n, sd = var(x))

samples <- polynomial(x) + whiteNoise
dataset <- data.frame(x = x, f = samples)ぱ

library(glmnet)

library(ggplot2)

##過適合を回避する例##

set.seed(7)

polynomial <- function(x) {

return(x * (x + 2) * (x + 1))

}

n <- 30

x <- runif(n, -2, 1)

whiteNoise <- rnorm(n, sd = var(x))

samples <- polynomial(x) + whiteNoise

dataset <- data.frame(x = x, f = samples)ぱ

パッケージggplot2とglmnetを使うので事前にinstall.packagesを使ってインストールしておいてください。

datasetをプロットすると分かるように、下の画像のように多項式x * (x + 2) * (x + 1)の周辺にデータがサンプリングされています。

ggplot(dataset, aes(x = x, y = f)) + 
  geom_point() + 
  geom_line(data = data.frame(x = seq(-2, 1, 0.001), f = polynomial(seq(-2, 1, 0.001))),
            aes(x = x, y = f), size = 1)

ggplot(dataset, aes(x = x, y = f)) +

geom_point() +

geom_line(data = data.frame(x = seq(-2, 1, 0.001), f = polynomial(seq(-2, 1, 0.001))),

aes(x = x, y = f), size = 1)

散布図

次に乱数を用いてトレインデータとテストデータを4:1の割合で作成します。

trainID <- sample(seq_len(nrow(dataset)), (4 / 5) * nrow(dataset))
testID <- setdiff(seq_len(nrow(dataset)), trainID)

dataset$validation <- NA
dataset$validation[trainID] <- "train"
dataset$validation[testID] <- "test"
dataset$validation <- as.factor(dataset$validation)                  
 
trainData <- dataset[trainID, ] 
testData <- dataset[testID, ]

trainID <- sample(seq_len(nrow(dataset)), (4 / 5) * nrow(dataset))

testID <- setdiff(seq_len(nrow(dataset)), trainID)

dataset$validation <- NA

dataset$validation[trainID] <- "train"

dataset$validation[testID] <- "test"

dataset$validation <- as.factor(dataset$validation)

trainData <- dataset[trainID, ]

testData <- dataset[testID, ]

多項式回帰を用いてdatasetのfをxで予測していきます。次のように、関数lmを用いることで多項式回帰を実行できます。

#多項式回帰
regression <- lm(f ~ poly(x, 10, raw = TRUE), data = trainData)

1 2	#多項式回帰 regression <- lm(f ~ poly(x, 10, raw = TRUE), data = trainData)

事前に用意しておいたtestDataに対するモデルの予測精度を検証してみます。関数predictを用いて予測値を計算し平均二乗誤差とテストデータの散布図と回帰曲線は次のようになります。

predictedData_curve <- data.frame(x = seq(-2, 1, 0.001), 
                                  predict = predict(regression, data.frame(x = seq(-2, 1, 0.001))))
MSE_train <- mean((predict(regression, trainData) - trainData$x)^2)
MSE_test <- mean((predict(regression, testData) - testData$x)^2)

ggplot(dataset) + 
  geom_point(aes(x = x, y = f, color = validation)) +
  geom_line(data = predictedData_curve, aes(x = x, y = predict), size = 1) +
  annotate("text", x = -1, y = 5.5, label = paste0("MSE for trainData: ", MSE_train)) +
  annotate("text", x = -1, y = 5, label = paste0("MSE for testData: ", MSE_test))

predictedData_curve <- data.frame(x = seq(-2, 1, 0.001),

predict = predict(regression, data.frame(x = seq(-2, 1, 0.001))))

MSE_train <- mean((predict(regression, trainData) - trainData$x)^2)

MSE_test <- mean((predict(regression, testData) - testData$x)^2)

ggplot(dataset) +

geom_point(aes(x = x, y = f, color = validation)) +

geom_line(data = predictedData_curve, aes(x = x, y = predict), size = 1) +

annotate("text", x = -1, y = 5.5, label = paste0("MSE for trainData: ", MSE_train)) +

annotate("text", x = -1, y = 5, label = paste0("MSE for testData: ", MSE_test))

多項式回帰

グラフから分かること

トレインデータに対する平均二乗誤差は小さいのですが、反対にテストデータに対する誤差は著しく大きいことが分かります。これがいわゆる過適合です。

これを回避する手法の一つして正則化があり、今回はL2正則化（リッジ回帰）を適用したいと思います。リッジ回帰はパッケージglmnetの関数glmnetの引数alphaを0にすることで実行できます。

#L2正則化
X_train_modelMatrix <- model.matrix(f ~ poly(x, 10, raw = TRUE), data = trainData)[, -1]

regression_ridge_cv <- cv.glmnet(X_train_modelMatrix, trainData$f, alpha = 0, nlambda = 100, nfold = 5)
regression_ridge <- glmnet(X_train_modelMatrix, trainData$f, alpha = 0, lambda = regression_ridge_cv$lambda.min)

#L2正則化

X_train_modelMatrix <- model.matrix(f ~ poly(x, 10, raw = TRUE), data = trainData)[, -1]

regression_ridge_cv <- cv.glmnet(X_train_modelMatrix, trainData$f, alpha = 0, nlambda = 100, nfold = 5)

regression_ridge <- glmnet(X_train_modelMatrix, trainData$f, alpha = 0, lambda = regression_ridge_cv$lambda.min)

通常の多項式回帰はpolyを用いて多項式を表現しましたが、glmnetの引数xには多項式に対応する計画行列を代入する必要があるため、2行目でmodel.matrixを用いて計画行列を作っています。また、4行目でcv.glmnetを用いてますが、これは最適なハイパーパラメータを決定するためです。クロスバリデーションにより最適なλを決定し、この値をもとに5行目でリッジ回帰を実行しています。今回はサンプル数が少ないのでフォールド数を5にしています（nfold = 5）。

また、plotを用いることでハイパーパラメータが変化したときの各回帰係数の挙動を確認することもできます。

plot(glmnet(X_train_modelMatrix, trainData$f, alpha = 0), xvar = "lambda")
legend("topright", legend = paste0("beta", seq(10)), lty = 1, col = seq(10))

1 2	plot(glmnet(X_train_modelMatrix, trainData$f, alpha = 0), xvar = "lambda") legend("topright", legend = paste0("beta", seq(10)), lty = 1, col = seq(10))

L2正則化　回帰係数

関数cv.glmnetの結果をプロットすると、次のハイパーパラメータと誤差の関係がプロットされます。

L2正則化　誤差

リッジ回帰の結果は次のようにregression_ridgeの後ろに$を付けることで取得できます。

a0 <- regression_ridge$a0
betas <- regression_ridge$beta

resultTable <- rbind(a0, betas)
write.csv(resultTable, "L2正則化_リッジ回帰.csv")

a0 <- regression_ridge$a0

betas <- regression_ridge$beta

resultTable <- rbind(a0, betas)

write.csv(resultTable, "L2正則化_リッジ回帰.csv")

上を実行すると作業ディレクトリに次の画像のcsvファイルが作成されます。

リッジ回帰結果

最後に、リッジ回帰の予測精度をみていきます。次を実行するとリッジ回帰の平均二乗誤差と回帰曲線をプロットすることができます。

X_test_modelMatrix <- model.matrix(f ~ poly(x, 10, raw = TRUE), data = testData)[, -1]
modelMatrix_curve <- model.matrix(f ~ poly(x, 10, raw = TRUE), data = data.frame(x = seq(-2, 1, 0.001), f = polynomial(seq(-2, 1, 0.001))))[, -1]

predictedData_curve <- data.frame(x = seq(-2, 1, 0.001), 
                                  predict = c(predict(regression_ridge, modelMatrix_curve)))
MSE_train_ridge <- mean((predict(regression_ridge, X_train_modelMatrix) - trainData$x)^2)
MSE_test_ridge <- mean((predict(regression_ridge, X_test_modelMatrix) - testData$x)^2)

ggplot(dataset) +
  geom_point(aes(x = x, y = f, color = validation)) +
  geom_line(data = predictedData_curve, aes(x = x, y = predict), size = 1) +
  annotate("text", x = -1, y = 5.5, label = paste0("MSE for trainData: ", MSE_train_ridge)) +
  annotate("text", x = -1, y = 5, label = paste0("MSE for testData: ", MSE_test_ridge))

X_test_modelMatrix <- model.matrix(f ~ poly(x, 10, raw = TRUE), data = testData)[, -1]

modelMatrix_curve <- model.matrix(f ~ poly(x, 10, raw = TRUE), data = data.frame(x = seq(-2, 1, 0.001), f = polynomial(seq(-2, 1, 0.001))))[, -1]

predictedData_curve <- data.frame(x = seq(-2, 1, 0.001),

predict = c(predict(regression_ridge, modelMatrix_curve)))

MSE_train_ridge <- mean((predict(regression_ridge, X_train_modelMatrix) - trainData$x)^2)

MSE_test_ridge <- mean((predict(regression_ridge, X_test_modelMatrix) - testData$x)^2)

ggplot(dataset) +

geom_point(aes(x = x, y = f, color = validation)) +

geom_line(data = predictedData_curve, aes(x = x, y = predict), size = 1) +

annotate("text", x = -1, y = 5.5, label = paste0("MSE for trainData: ", MSE_train_ridge)) +

annotate("text", x = -1, y = 5, label = paste0("MSE for testData: ", MSE_test_ridge))

リッジ回帰

グラフから分かること

画像から分かるように、多項式回帰のときと比べてテストデータに対する平均二乗誤差が大幅に減少していることが分かります。

回帰係数が小さくなるような最適化の条件を与えることで、トレインデータへの過適合を回避しています。

ラッソ回帰とリッジ回帰の比較（L1, L2正則化の比較）

次にラッソ回帰の実行例を紹介します。ラッソ回帰とリッジ回帰の比較を行いたいと思います。

使用するデータセットは次のBostonです。BostonはパッケージMASSに入っているので事前にinstall.packagesでインストールしてください。

##L1,L2正則化の比較##
library(MASS)
set.seed(0)

dataset <- Boston

##L1,L2正則化の比較##

library(MASS)

set.seed(0)

dataset <- Boston

Bostonは次の14個の変数から成ります。一番最後の住宅価格medvを他の変数で予測したいと思います。

> head(dataset)
     crim zn indus chas   nox    rm  age    dis rad tax ptratio  black lstat medv
1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90  4.98 24.0
2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90  9.14 21.6
3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83  4.03 34.7
4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63  2.94 33.4
5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90  5.33 36.2
6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12  5.21 28.7

> head(dataset)

crim zn indus chas nox rm age dis rad tax ptratio black lstat medv

1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 15.3 396.90 4.98 24.0

2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90 9.14 21.6

3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83 4.03 34.7

4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63 2.94 33.4

5 0.06905 0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 396.90 5.33 36.2

6 0.02985 0 2.18 0 0.458 6.430 58.7 6.0622 3 222 18.7 394.12 5.21 28.7

まず、リッジ回帰を実行したいと思います。次を実行するとリッジ回帰を実行することができます。

#L2正則化 Ridge
X_modelMatrix <- model.matrix(medv ~., data = dataset)
regression_ridge_cv <- cv.glmnet(X_modelMatrix, dataset$medv, alpha = 0)
regression_ridge <- glmnet(X_modelMatrix, dataset$medv, alpha = 0, lambda = regression_ridge_cv$lambda.min)

a0 <- regression_ridge$a0
betas <- regression_ridge$beta

resultTable_ridge <- rbind(a0, as.matrix(betas))

#L2正則化 Ridge

X_modelMatrix <- model.matrix(medv ~., data = dataset)

regression_ridge_cv <- cv.glmnet(X_modelMatrix, dataset$medv, alpha = 0)

regression_ridge <- glmnet(X_modelMatrix, dataset$medv, alpha = 0, lambda = regression_ridge_cv$lambda.min)

a0 <- regression_ridge$a0

betas <- regression_ridge$beta

resultTable_ridge <- rbind(a0, as.matrix(betas))

次にラッソ回帰を行います。ラッソ回帰は、関数glmnetの引数alphaを1に指定することで行うことができます。

#L1正則化 Lasso
regression_lasso_cv <- cv.glmnet(X_modelMatrix, dataset$medv, alpha = 1)
regression_lasso <- glmnet(X_modelMatrix, dataset$medv, alpha = 1, lambda = regression_lasso_cv$lambda.min)

#L1正則化 Lasso

regression_lasso_cv <- cv.glmnet(X_modelMatrix, dataset$medv, alpha = 1)

regression_lasso <- glmnet(X_modelMatrix, dataset$medv, alpha = 1, lambda = regression_lasso_cv$lambda.min)

ラッソ回帰の回帰係数とハイパーパラメータの関係をプロットすると次の画像のようになります。

plot(glmnet(X_modelMatrix, dataset$medv, alpha = 1,), xvar = "lambda")
legend("bottomright", legend = rownames(regression_lasso$beta), lty = 1, col = seq(ncol(dataset) - 1))

1 2	plot(glmnet(X_modelMatrix, dataset$medv, alpha = 1,), xvar = "lambda") legend("bottomright", legend = rownames(regression_lasso$beta), lty = 1, col = seq(ncol(dataset) - 1))

L1正則化　回帰係数

リッジ回帰と同様に回帰係数をデータフレームとして保存します。

a0 <- regression_lasso$a0
betas <- regression_lasso$beta

resultTable_lasso <- rbind(a0, as.matrix(betas))

a0 <- regression_lasso$a0

betas <- regression_lasso$beta

resultTable_lasso <- rbind(a0, as.matrix(betas))

リッジ回帰とラッソ回帰の回帰係数を比較してみます。

> resultTable <- data.frame(resultTable_ridge, resultTable_lasso)
> resultTable
                    Ridge         Lasso
a0           28.051686825  34.905311445
(Intercept)   0.000000000   0.000000000
crim         -0.087904300  -0.100851620
zn            0.032606201   0.042632690
indus        -0.038328191   0.000000000
chas          2.902980774   2.692311675
nox         -12.005369287 -16.565500129
rm            4.014163735   3.849597177
age          -0.003862644   0.000000000
dis          -1.120903830  -1.420950722
rad           0.154161048   0.264420345
tax          -0.005729860  -0.010322154
ptratio      -0.855862908  -0.933656448
black         0.009068108   0.009088836
lstat        -0.471596371  -0.522600774