\[ % 定義 \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \]

1 本講義の目的

  • 正則化
  • 交差検証
  • ハイパーパラメータ

について学ぶ

2 パッケージの読み込み

今回の分析で使うパッケージを読み込む。

pacman::p_load(tidyverse, 
               MLmetrics, # 予測精度の評価指標
               patchwork, # グラフ結合用
               gridExtra, # グラフ結合用
               glmnet) # lassoやリッジ回帰のライブラリ

3 汎化性能と過学習

機械学習で行いたいことは学習に使用したデータのみへの予測精度を高めることではなく、学習に使用していない未知のデータに対しても高い予測精度をもつこと、すなわち汎化性能を高めることである。

逆に、学習用データに対しての予測精度は高いものの、未知のデータに対する予測精度が低い（汎化性能が低い）予測モデルのことを過学習（overfitting）しているという。

数値例を示そう。ある変数\(\{x,y\}\)について50個のデータが得られ、そのうち25個を学習用データにして予測モデルを作るとする。

# seedを固定
set.seed(11111)

# データの生成
n <- 50
x <- 1:n * 1/n
y <- 0.5 + 0.4*sin(2*pi*x) + rnorm(n = n, sd = 0.1)
df <- tibble(y,
             x, 
             x2 = x^2, 
             x3 = x^3, 
             x4 = x^4, 
             x5 = x^5, 
             x6 = x^6, 
             x7 = x^7, 
             x8 = x^8, 
             x9 = x^9)


# データの分割
# ID列を追加
df <- df %>% rownames_to_column("ID")

# 50%を学習用データに
train <- df %>% sample_frac(size = 0.5)

# 学習用データに使っていないIDの行をテスト用データに
test <- anti_join(df, train, by = "ID")

# ID列の削除
df <- df %>% dplyr::select(-ID)
train <- train %>% dplyr::select(-ID)
test <- test %>% dplyr::select(-ID)
train %>% ggplot(aes(x, y)) +
  geom_point()

データは非線形な分布をしているため、以下のような多項式回帰をいくつか作って比較してみるとする。

\[ y = \beta_0 + \beta_1 x + \beta_2 x^2 + \cdots + \beta_p x^p + \varepsilon \]

以下の図が次数をいくつか変えていった多項式回帰の回帰曲線である。

# plot
plots <- list() # グラフを格納するリストを用意
orders <- c(1, 2, 3, 5, 7, 9)  # 試していきたい多項式回帰の次数を用意
for (i in 1:length(orders)) {
  # 訓練データと訓練データにfitさせた多項式回帰のグラフを描き、plotsに格納する
  plots[[i]] <- train %>% ggplot(aes(x, y)) +
    geom_point() +
    geom_smooth(method = "lm", se = F, 
                formula = as.formula(str_c("y ~ poly(x, degree = ", orders[i], ")"))) +
    labs(title = str_c(orders[i], "次多項式回帰"))
}

# gridExtraを使って複数のグラフをまとめる
gridExtra::grid.arrange(grobs = plots, ncol = 3)

# patchworkを使って複数のグラフをまとめる場合
# (plots[[1]] + plots[[2]] + plots[[3]]) / (plots[[4]] + plots[[5]] + plots[[6]])

次数の高いモデルのほうがデータの分布をよく近似する曲線を描いている。

では次数が最も高い多項式回帰が最も良い汎化性能を持つのかというと、そうとは限らない。

次数の異なる多項式回帰モデルごとに学習用データでの誤差とテストデータでの誤差を測ってみよう。

# 次数を変えながら学習・予測
for(p in 1:9){
  reg <- lm(y ~ ., data = train[,1:(p+1)]) # 学習
  # trainへの予測精度
  y_train_pred <- predict(reg, train[,1:(p+1)])   # 予測
  train_error <- RMSE(y_train_pred, train$y) # 誤差
  # testへの予測精度
  y_test_pred <- predict(reg, test[,1:(p+1)])   # 予測
  test_error <- RMSE(y_test_pred, test$y) # 誤差
  if(p == 1) result <- tibble(p, train_error, test_error) # p=1の場合、データフレーム作る
  if(p != 1) result <- rbind(result, tibble(p, train_error, test_error)) # p=2以上の場合はrbindで繋げる
}

# ggplotで多項式回帰と誤差(RMSE)の関係を描画
ggplot(result %>% gather(key, error, -p),
       aes(x = p, 
           y = error, 
           color = key)) +
  geom_point() +
  geom_line(size = 1) +
  scale_x_continuous(breaks = c(1:9)) +
  labs(color = "", 
       x = "多項式回帰の次数", 
       y = "誤差（RMSE)",
       title = "モデルの複雑さと過学習")

多項式回帰の次数が高くなるほど学習用データへの予測誤差（train_error）は減少しているが、テストデータに対する予測誤差（test_error）は次数が5以上になると増加しており、過学習していることがわかる。

次数が高い多項式回帰モデルは高い表現力をもつ複雑なモデルであるが、データに対して過剰に複雑なモデルを使用してしまうと過学習をおこすことになる。

過学習を抑えたい場合に使われる手法の1つに正則化（regularization）という手法がある。以下では有名な正則化回帰であるリッジ回帰とlassoについて述べる。

4 リッジ回帰

4.1 リッジ回帰

リッジ回帰（ridge regression）は推定するパラメータの二乗和がある値\(R\)を超えない範囲で推定するという制約条件の下で残差二乗和\(\sum^n_{i=1}(y_i - \hat{y}_i)^2\)（\(\hat{y}\)は予測値）を最小化するようなパラメータを推定する。

\[ \min_\beta \sum^n_{i=1}(y_i - \hat{y}_i)^2 ,\hspace{1em} \text{subject to }\sum_{j=1}^p \beta_j^2 \leq R \]

このことについて、パラメータが２つ（\(\beta_1, \beta_2\)）の場合を例に最小二乗法による推定値(\((\hat{\beta}_1^{LS}, \hat{\beta}_2^{LS})\))とリッジ回帰による推定値(\((\hat{\beta}_1^{ridge}, \hat{\beta}_2^{ridge})\))を視覚的に比較したものが下の図である。

最小化したい残差二乗和は放物面のような形になっており（左図）、それを最小化するパラメータが最小二乗推定量\((\hat{\beta}_1^{LS}, \hat{\beta}_2^{LS})\)となる。

左図を真上から見るようにして２次元の図にしたものが中央の図で、楕円は放物面の等高線（誤差二乗和が等しい点を結んだ線）を表している。

最小二乗法は採用するパラメータに制約が無いためパラメータ空間全体を使うことができる（中央図）。

リッジ回帰は制約条件を満たす範囲がパラメータ空間の原点を中心とする円になるため、円の領域の中で残差二乗和を最小化するパラメータを探索することになる（右図）。制約条件の\(R\)は円の大きさを表している。

なお、最小二乗法で得られる推定量がリッジ回帰の円の中にある場合は最小二乗推定量とリッジ回帰の推定量は同じものになる。

書き方を変えると、リッジ回帰のパラメータの推定量（のベクトル）\(\hat{\boldsymbol{\beta}}^{ridge}\)は \[ \hat{\boldsymbol{\beta}}^{ridge} = \argmin_{\boldsymbol{\beta}} \left\{\sum_{i=1}^n(y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^p \beta_j^2\right\} \] というものになる。

最小二乗法の推定量は残差二乗和\(\sum_{i=1}^n(y_i - \hat{y}_i)^2 \)を最小化するパラメータ \[ \hat{\boldsymbol{\beta}}^{OLS} = \argmin_{\boldsymbol{\beta}} \sum_{i=1}^n(y_i - \hat{y}_i)^2 \]

であった。リッジ回帰は残差二乗和に\(\lambda \sum_{j=1}^p \beta_j^2\)の項を加えたものを誤差関数とし、それを最小化している。

\(\lambda\)は正則化パラメータ（regularization parameter）と呼ばれ、正則化項（regularization term）\(\sum_{j=1}^p \beta_j^2\)の影響の程度を調整するパラメータである。

\(\lambda\)は予測モデルによってデータから推定されるパラメータ\(\beta_j \)とは異なり、分析者が適切な値を指定する必要がある。こうした予測モデル外から設定する必要があるパラメータをハイパーパラメータ（hyperparameter）という。

4.2 データの標準化

変数の測定単位が異なる場合、回帰分析のパラメータ\(\beta_j\)も異なってくる。

リッジ回帰は正則化項\(\sum_{j=1}^p \beta_j^2\)でパラメータの値に基づいて推定を行うため、正則化法を使用する際には各変数の測定単位の違いやパラメータ\(\beta\)の大きさの差をなくすために説明変数を標準化する必要がある。

標準化とは、変数の平均を0、分散を1にすることであり、ある変数\(x\)の平均\(\mu\)と標準偏差\(\sigma\)を用いて

\[ z = \frac{x - \mu}{\sigma} \]

のように変換することである。

適当に乱数を生成して比較したものが次のグラフである。標準化することによって分布の中心が原点\((0,0)\)になる。これを中心化という。

Rではscale()関数で標準化を行うことができる。

4.3 Rでの実践

冒頭で作成したデータを用いて検証していく。

リッジ回帰を実行する関数glmnet()では説明変数と目的変数を分ける必要があるので、あらかじめ分けておく

# glmnetは説明変数Xと目的変数yを別々に入れる必要があるので分ける

# train
X_train <- train %>% select(-y) %>% as.matrix() # as.matrix()でmatrix型のデータにする
y_train <- train$y

# test
X_test <- test %>% select(-y) %>% as.matrix()
y_test <- test$y

4.3.1 glmnet関数

{glmnet}パッケージのglmnet()関数で正則化回帰を行うことができる。

glmnet()関数は以下のように書く。

glmnet(x = X, y = y, lambda = 1, alpha = 0)
  • x：説明変数。データフレームではなくmatrix型で指定する必要がある。
    • 説明変数のみを含むデータフレームをas.matrix()関数に通せばよい
  • y：被説明変数（目的変数）。ベクトル型で指定する。
  • lambda：正則化パラメータ\(\lambda\)の値
  • alpha：モデルの指定。alpha = 0でリッジ回帰になる。

なお、事前に説明変数の標準化を行うかどうかを指定するstandardizeという引数もあり、デフォルトではTRUEになっているため自動的に標準化を行ってくれる。

# 9次多項式でリッジ回帰
reg_ridge <- glmnet(x = X_train, # 特徴量
                    y = y_train, # 目的変数
                    lambda = 0.1, # λ：任意の値
                    alpha = 0) # alpha = 0でリッジ回帰

# trainに対する予測精度
y_train_pred <- predict(reg_ridge, X_train) # 予測
train_error <- RMSE(y_train_pred, y_train) # 誤差
train_error %>% round(3) # 結果を表示
## [1] 0.157
# testに対する予測精度
y_test_pred <- predict(reg_ridge, X_test) # 予測
test_error <- RMSE(y_test_pred, y_test) # 誤差
test_error %>% round(3) # 結果を表示
## [1] 0.183

学習用データに対する予測誤差（RMSE）は0.157、テストデータに対しては0.183となった。

この結果は正則化パラメータ\(\lambda\)の値によって大きく変化する。次節では適切な\(\lambda\)の推定方法について述べる。

5 交差検証とハイパーパラメータのチューニング

5.1 交差検証

k-分割交差検証（k-fold cross validation）とは予測精度の測り方の一つであり、

  1. データセットを\(k\)個（例えば5個や10個）のサブセットへ無作為に分割する
  2. \(k-1\)個のサブセットのデータで予測モデルを学習させる
  3. 残りの\(1\)個のサブセットのデータをテスト用に回して予測誤差を測る
  4. これを\(k\)回繰り返して予測誤差の平均値を使って予測性能を評価する

という手順で予測精度を測定するものである。

これまではデータを学習用データ・テスト用データに2分割して予測性能を評価してきたが（この方法をホールドアウト法という）、 交差検証は手元のデータすべてを学習とテストに使用できるため、データが少ない場合においてはホールドアウト法よりもよい予測性能の評価方法になる。

\(k\)をサンプルサイズと同じにした場合、つまり1レコードだけをテストデータに回して残りの\(n-1\)レコードを訓練データにする交差検証を1つ抜き交差検証（leave-one-out cross validation: LOOCV）と呼ぶ。

5.2 cv.glmnet()関数によるハイパーパラメータの探索

機械学習の分野においては基本的にパラメータ（parameter）という言葉はデータによって学習（推定）される対象を指す。一方で正則化パラメータ\(\lambda\)のような値はデータから推定されるものではなく、学習に先立って分析者が設定していくことになる。このような値のことをハイパーパラメータ（hyperparameter）と呼ぶ。

ハイパーパラメータの調整（tuning）は、ハイパーパラメータを変えて学習と予測精度の評価を行う実験を繰り返して最適値を探して調整していくことになる。

{glmnet}パッケージには正則化パラメータ\(\lambda\)を少しずつ変えながら交差検証で予測性能を評価することにより良い予測性能を与える\(\lambda\)を推定するための関数cv.glmnet()があるため、簡単に調整を行うことができる（なお、cv=cross validation=交差検証である）。

# ラムダを変えながら交差検証で誤差を評価
ridge_cv <- cv.glmnet(x = X_train, # 説明変数
                      y = y_train, # 目的変数
                      nfolds = nrow(X_train), # 1つ抜き交差検証
                      grouped = FALSE,
                      lambda = seq(0.000001, 0.001, length.out = 200), # 探索する範囲を指定（任意）
                      alpha = 0) # alpha = 0でリッジ回帰

# ラムダと交差検証誤差のグラフ
plot(ridge_cv)

plot()で\(\lambda\)を変えながら交差検証で誤差（MSE)を計算したときのグラフを見ることができる。縦の破線は最適と考えられる\(\lambda\)の値を示している。左側の破線は交差検証でMSEが最小だった\(\lambda\)であり、右側の破線は「最小のMSE+その標準誤差未満」の範囲で最大の\(\lambda\)である。

（後者を選択するというのは、交差検証のMSEが小さい中でもできるだけ\(\lambda\)が大きい（正則化項の影響が大きい）モデルを採用して過学習を回避しようとする\(\lambda\)の選択方法である。これは1標準誤差ルール（one standard error rule）と呼ばれる。）

5.3 推定値

coef()関数で推定されたパラメータを見ることができる。その際にs = "lambda.min"と指定することで交差検証でMSEが最小だった\(\lambda\)にアクセスできる。

# 誤差(MSE)を最小にしたラムダを採択したときの係数推定値
coef(ridge_cv, s = "lambda.min") %>% round(3)
## 10 x 1 sparse Matrix of class "dgCMatrix"
##                 s1
## (Intercept)  0.609
## x            1.795
## x2          -2.579
## x3          -4.694
## x4          -0.956
## x5           3.608
## x6           5.732
## x7           4.398
## x8          -0.202
## x9          -7.429
# (参考)MSE+1SEの中で最大のラムダを採用したときの係数推定値
coef(ridge_cv, s = "lambda.1se") %>% round(3)
## 10 x 1 sparse Matrix of class "dgCMatrix"
##                 s1
## (Intercept)  0.732
## x            1.024
## x2          -2.676
## x3          -1.858
## x4           0.071
## x5           1.686
## x6           2.349
## x7           1.840
## x8           0.163
## x9          -2.576

5.4 予測精度

続いて予測精度を評価してみる。

predict()でも同様にs = "lambda.min"と指定することで、誤差を最小化していた\(\lambda\)を採用したときの予測値を算出できる。

# trainに対する予測精度
y_train_pred <- predict(ridge_cv, 
                        X_train, 
                        s = "lambda.min") # 予測

train_error <- RMSE(y_train_pred, y_train) # 誤差
train_error %>% round(3)
## [1] 0.091
# testに対する予測精度
y_test_pred <- predict(ridge_cv, 
                       X_test, 
                       s = "lambda.min") # 予測

test_error <- RMSE(y_test_pred, y_test) # 誤差
test_error %>% round(3)
## [1] 0.117

学習用データに対する予測誤差（RMSE）は0.091、テストデータに対しては0.117となった。

リッジ回帰でなく最小二乗法で9次多項式回帰を推定した際のテストデータに対する誤差が0.596であったのと比べると、過学習が抑制されている。

テスト用データに対する予測値のグラフをみても、最小二乗法（OLS）は（学習用データに対してぴったり予測するように過学習しているため）極端に外している値があるのに対し、リッジ回帰（Ridge）は9次多項式回帰であっても穏当な誤差になっている様子がみられる。

# testへの9次多項式OLS
reg_ols <- lm(y ~ ., train)
y_test_pred_ols <- predict(reg_ols, test)

# data
data_for_plot <- tibble(x = test$x, 
                        y = test$y,
                        OLS = y_test_pred_ols,
                        Ridge = c(y_test_pred)) 

#long形式に変換
data_for_plot <- data_for_plot %>%
  tidyr::pivot_longer(OLS:Ridge,
                      names_to = "model",
                      values_to = "value")

# plot
data_for_plot %>% ggplot(aes(x = x,
                             y = value,
                             color = model)) +
  geom_point(aes(x = x,
                 y = y),
             size = 1.5, 
             color = "black") +
  geom_line(size = 1.5, alpha = 0.5) +
  geom_point(size = 2, alpha = 0.5) +
  scale_color_brewer(palette = "Set1")+
  labs(title = "9次多項式回帰によるtestデータに対する予測",
       y = "y and predicted y")

6 lasso

6.1 特徴量選択

本資料の序盤で述べた多項式回帰の例のように、無駄な説明変数が多く過度に複雑なモデルは過学習を起こしやすい。

より良い予測モデルを構築するには、適切な説明変数（機械学習の分野では特徴量 featureとも呼ばれる）を投入する必要がある。

そのためのプロセスは特徴量エンジニアリング（feature engineering）と呼ばれ、次の作業から構成される。

  • 特徴量選択（feature selection）：既存の特徴量から予測モデルの構築に役立つ特徴量を選択する。
  • 特徴量抽出（feature extraction）：既存の特徴量を組み合わせて、より役に立つ１つの特徴量を生み出す。
  • 新しいデータの収集による特徴量の追加

次節で紹介するlassoは特徴量エンジニアリングのうち特徴量選択をデータから自動的に行うことができる手法である。

6.2 lasso

lasso（least absolute shrinkage and selection operator）は \[ \min_\beta \sum^n_{i=1}(y_i - \hat{y}_i)^2 ,\hspace{1em} \text{subject to }\sum_{j=1}^p |\beta_j| \leq R \]

の制約付き最小化問題の解によってパラメータを推定するもので、そのパラメータは

\[ \hat{\boldsymbol{\beta}}^{lasso} = \argmin_{\boldsymbol{\beta}} \left\{\sum_{i=1}^n(y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^p |\beta_j|\right\} \]

となる。上記の式は、リッジ回帰のものと非常に似ているが、正則化項が\(\sum_{j=1}^p \beta_j^2\)ではなく\(\sum_{j=1}^p |\beta_j|\)である点が異なっている。なお、前者を\(L_2\)正則化、後者を\(L_1\)正則化と呼ぶ。

最小二乗法による推定値(\((\hat{\beta}_1^{LS}, \hat{\beta}_2^{LS})\))とlassoによる推定値(\((\hat{\beta}_1^{lasso}, \hat{\beta}_2^{lasso})\))を視覚的に比較したものが下の図である。

lassoの制約条件を満たす推定量はパラメータ空間の原点を中心にしたひし形の領域になる（右図）。

もし最小二乗推定量\((\hat{\beta}_1^{LS}, \hat{\beta}_2^{LS})\)がlassoのひし形の制約領域の外にある場合、制約領域内で残差二乗和を最小化するパラメータがlassoの推定量となるため、残差二乗和の等高線の楕円とlassoのひし形の領域が接する点がlassoの推定量となる。残差二乗和の等高線の形状によってはlassoの推定量が\((0, \hat{\beta}_2^{lasso})\)のように一部のパラメータが0と推定される。

0と推定されたパラメータに対応する説明変数は目的変数に寄与しないと解釈することができ、特徴量の選択が行われたモデルを得ることができる。リッジ回帰はパラメータを縮小するものの、ぴったり0に縮小することはない。パラメータをぴったり0に推定することで特徴量選択を行うことができるのはlassoの長所である。

なお、正則化パラメータ\(\lambda\)を大きくすると0と推定されるパラメータの数が多くなり、\(\lambda \)を小さくすると0と推定されるパラメータの数は少なくなる。

6.3 Rでの実践

リッジ回帰のときと同じデータ（xの9次多項式のデータ）でlassoを実行してみよう。

glmnet()関数で引数にalpha = 1と指定すればlassoを実行できる。

まずcv.glmnet()関数で適切な\(\lambda\)を推定する

# ラムダを変えながら交差検証で誤差を評価
lasso_cv <- cv.glmnet(x = X_train, # 特徴量
                      y = y_train, # 目的変数
                      nfolds = nrow(X_train), # 1つ抜き交差検証
                      grouped = FALSE,
                      lambda = seq(-20, 0, 0.5) %>% exp(), # 探索する範囲を指定（任意）
                      alpha = 1) # alpha = 1でlasso

# ラムダと交差検証MSEのグラフ
plot(lasso_cv)

推定されたパラメータを確認する

# 誤差(MSE)を最小にしたラムダを採択したときの係数推定値
coef(lasso_cv, s = "lambda.min") %>% round(3)
## 10 x 1 sparse Matrix of class "dgCMatrix"
##                 s1
## (Intercept)  0.638
## x            1.402
## x2          -1.139
## x3          -6.361
## x4           .    
## x5           0.000
## x6          11.714
## x7           1.821
## x8           .    
## x9          -7.827

いくつかのパラメータはぴったり0と推定されたため説明変数として採用されなかったことが.で示されている。

# (参考) MSE + 1SEの中で最大のラムダを採用したときの係数推定値
coef(lasso_cv, s = "lambda.1se") %>% round(3)
## 10 x 1 sparse Matrix of class "dgCMatrix"
##                 s1
## (Intercept)  0.632
## x            2.030
## x2          -5.223
## x3           .    
## x4           .    
## x5           .    
## x6           6.453
## x7           .    
## x8           .    
## x9          -3.289
# trainに対する予測精度
y_train_pred <- predict(lasso_cv, 
                        X_train, 
                        s = "lambda.min") # 予測

train_error <- RMSE(y_train_pred, y_train) # 誤差

str_c("train_error: ", round(train_error, 3)) %>% print() # 結果を表示
## [1] "train_error: 0.09"
# testに対する予測精度
y_test_pred <- predict(lasso_cv, 
                      X_test, 
                      s = "lambda.min") # 予測

test_error <- RMSE(y_test_pred, y_test) # 誤差(MSE)

str_c("test_error: ", round(test_error, 3)) %>% print() # 結果を表示
## [1] "test_error: 0.12"

テスト用データに対する予測の誤差は0.12となった。

最小二乗法で9次多項式回帰を推定した際のテストデータに対する誤差が0.596であったのと比べると、過学習が抑制されており、よりよい予測モデルになっている。

7 まとめ

7.1 推定量

最小二乗法は残差二乗和\(\sum_{i=1}^n(y_i - \hat{y}_i)^2\)という誤差関数を最小化するパラメータ\(\hat{\boldsymbol{\beta}}^{LS}\)を採用するというものだった。

\[ \hat{\boldsymbol{\beta}}^{LS} = \argmin_{\boldsymbol{\beta}} \sum_{i=1}^n(y_i - \hat{y}_i)^2 \]

リッジ回帰は残差二乗和に\(\lambda \sum_{j=1}^p \beta_j^2\)を加えた誤差関数を最小化するパラメータ\(\hat{\boldsymbol{\beta}}^{ridge}\)を採用するというものだった。

\[ \hat{\boldsymbol{\beta}}^{ridge} = \argmin_{\boldsymbol{\beta}} \left\{\sum_{i=1}^n(y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^p \beta_j^2\right\} \]

lassoは残差二乗和に\(\lambda \sum_{j=1}^p |\beta_j|\)を加えた誤差関数を最小化するパラメータ\(\hat{\boldsymbol{\beta}}^{lasso}\)を採用するというものだった。

\[ \hat{\boldsymbol{\beta}}^{lasso} = \argmin_{\boldsymbol{\beta}} \left\{\sum_{i=1}^n(y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^p |\beta_j|\right\} \]

7.2 推定結果

今回、例として使用した9次多項式回帰モデルでは、次のような推定結果になった。

まず、最小二乗法、リッジ回帰、lassoで推定されたパラメータは次のようになった。

OLS Ridge Lasso
(Intercept) 6.315 0.609 0.638
x -198.461 1.795 1.402
x2 2551.378 -2.579 -1.139
x3 -16400.664 -4.694 -6.361
x4 60332.666 -0.956 0.000
x5 -134969.809 3.608 0.000
x6 186701.721 5.732 11.714
x7 -155990.432 4.398 1.821
x8 72186.258 -0.202 0.000
x9 -14219.122 -7.429 -7.827

テスト用データに対する予測値は次のようになった。

その予測誤差（RMSEで測ったもの）は次のようになった。

OLS Ridge Lasso
RMSE 0.596 0.117 0.12

8 （参考）エラスティックネット

リッジ回帰とlassoを結合したエラスティックネット（elastic net）というアルゴリズムも存在する。

\[ \hat{\boldsymbol{\beta}}^{EN} = \argmin_{\boldsymbol{\beta}} \left\{ \sum_{i=1}^n(y_i - \hat{y}_i)^2 + \lambda \sum_{j=1}^p \left ( \alpha |\beta_j| + (1-\alpha) \beta_j^2 \right ) \right\} \]

\(\alpha\)はlassoとリッジ回帰のどちらに重きを置くかのウェイトで、これもハイパーパラメータである。

実はglmnet()はエラスティックネットの関数であり、glmnet()を実行するときに、「リッジ回帰ならalpha = 0」と指定していたのがこの\(\alpha\)にあたる。

# Elastic Net

# ラムダを変えながら交差検証で誤差を評価
EN_cv <- cv.glmnet(x = X_train, # 特徴量
                   y = y_train, # 目的変数
                   alpha = 0.35) # 任意のα

# testに対する予測精度
y_test_pred <- predict(EN_cv, X_test, s = "lambda.min") # 予測
test_error <- RMSE(y_test_pred, y_test) # 誤差
test_error %>% round(3)
## [1] 0.12

9 参考文献

Jamesほか（2018）『Rによる統計的学習入門』、朝倉書店。 （原著ウェブサイト、原著PDF）

Hastieほか（2014）『統計的学習の基礎-データマイニング・推論・予測』、共立出版。 （原著ウェブサイト、原著PDF）

岩波データサイエンス刊行委員会編（2017）『岩波データサイエンス Vol.5 スパ―スモデリングと多変量データ解析』

川野ほか（2018） 『スパース推定法による統計モデリング』、共立出版。