Rでスパースモデリング：Elastic Net回帰についてまとめてみる

導入

回帰モデル構築の際、汎化性能を向上させるために正則化の手法がたびたび用いられます。これは、考えているデータ数に対して特徴量の数が非常に多い場合や、特徴量間に強い相関（多重共線性）がある場合に有効な方法となっています。このような場合に、通常の回帰モデル構築の際に用いられる2乗誤差などの目的関数に加え、Lp$L_p$ノルム（p$p$は正整数）のような正則化項（もしくは罰則項）加えて最適化をおこなうことで先程の問題を解消することができます。こういった正則化項を加えた上でモデルの最適化をおこなう（ = パラメータを推定する）方法を、正則化法といいます。

代表的な正則化法に、Lasso, Ridge, Elastic Net回帰があります。これらは、解釈性も含めた特徴があり、必ずしも高精度のものだからよいわけではない、というのが私の考えです。しかし一方で、{caret}を使ってこの中で最も精度がよいものを採用しました、ということを行っている状況も（私の周りでは）見受けられます。精度向上だけが目的ならこういった立場でも良いのですが、Lasso回帰を始めとするスパースモデリングの強みの一つである、「人間の事前知識/解釈性」をモデルに取り入れることに関しては、必ずしも活かせないことになります。機械が自動で色々いい感じなものを作ってくれる未来はまだまだ先であり、モデリングの際に、どういったことが目的で、それをどう実現していくかをある程度理論的な背景を踏まえた上でおこなっていくことがまだまだ重要だと私は考えています。

まずは代表的な正則化法に関して整理しておきたいと考えたので、それらに関してまとめたものを備忘録的に載せておきます。

使用するデータセット

ボストン近郊の住宅価格のデータを利用します。Rでは{MASS}パッケージを導入すると付随してこのデータセットを使えるようになるため、{MASS}パッケージを利用します。

# ボストン近郊の住宅価格のデータを利用するため（だけ）に導入
if(!require(MASS)){
  install.packages("MASS", quiet = TRUE)
}
require(MASS)


> head(Boston)
     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

目的変数はmedvで、説明変数候補はそれ以外を想定しています。具体的なモデル構築は、後ほどおこないます。

線形回帰モデル

以下のような線形回帰モデルを仮定します。

y⃗ =Xβ⃗ +ϵ⃗ (1)$$\begin{array}{}\text{(1)}& \overrightarrow{y}=X\overrightarrow{\beta }+\overrightarrow{ϵ}\end{array}$$

y, β⃗ , ϵ⃗ $y,\overrightarrow{\beta },\overrightarrow{ϵ}$はそれぞれ目的変数、回帰係数（推定するパラメータ）、残差で、X$X$は計画行列（説明変数をまとめた行列）です。目的変数と特徴量は、それぞれ中心化・標準化されているとします。データ数n$n$、特徴量の数p$p$とすると、y, ϵ⃗ $y,\overrightarrow{ϵ}$はn$n$次元ベクトル、β⃗ $\overrightarrow{\beta }$はp+1$p+1$次元ベクトル（切片項も含む）、X$X$はn×(p+1)$n\times (p+1)$行列となります。

回帰モデル構築の問題は、回帰係数β⃗ $\overrightarrow{\beta }$を求める最適化問題へと帰着されます。最適化のためのよく用いられる目的関数は、以下のような2乗誤差関数です*1。

S(β⃗ )=(y⃗ −Xβ⃗ )T(y⃗ −Xβ⃗ )(2)$$\begin{array}{}\text{(2)}& S\left(\overrightarrow{\beta }\right)={(\overrightarrow{y}-X\overrightarrow{\beta })}^T(\overrightarrow{y}-X\overrightarrow{\beta })\end{array}$$

この2乗誤差関数を最小にするようなパラメータβ^⃗ $\overrightarrow{\hat{\beta }}$は、実測とモデル式の誤差の2乗和を最小にするという意味で、最適なパラメータだと考えます。このパラメータは、目的関数をパラメータで微分したものを0とおいた方程式の解を求めればよく、簡単な計算で以下のようになることがわかります。

β^⃗ =(XTX)−1XTy⃗ (3)$$\begin{array}{}\text{(3)}& \overrightarrow{\hat{\beta }}={\left(X^{T}X\right)}^{-1}{X}^T\overrightarrow{y}\end{array}$$

式(3)の右辺は観測された量だけで構成されているため、それらを代入すれば回帰係数を求めることができます。

正則化法

先程の線形回帰モデルの議論は、特徴量とデータの関係によってはうまく機能しないことがあります。まず、どういった場合にどのような問題が生じるかを紹介し、それらをどのように解決していくか、という視点で正則化法の説明をしていきます。

多重共線性

式(3)で問題になるのは、多重共線性がある（特徴量間の相関が強い）場合です。

(XTX)$(X^{T}X)$の逆行列を求めるためには、(XTX)$(X^{T}X)$の行列式det(XTX)$\text{det}(X^{T}X)$を知る必要があります。しかし、多重共線性がある場合、(XTX)$(X^{T}X)$のなかにほとんど同じ値を持つような列が複数出現することになります。線形代数などで行列式を求める計算を思い出すと、全く同じ値をもった列があると行列式が0になり、行列が特異（正則行列ではない）になっている状態となります。つまり、多重共線性がある場合は、det(XTX)$\text{det}(X^{T}X)$が非常に0に近くなります。

(XTX)$(X^{T}X)$の逆行列は1/det(XTX)$1/\text{det}(X^{T}X)$に比例しているため、det(XTX)$\text{det}(X^{T}X)$が0に近づくと逆行列(XTX)−1$(X^{T}X{)}^{-1}$は無限大へ発散します。これは、 推定量の共分散β^⃗ =σ2(XTX)−1$\overrightarrow{\hat{\beta }}={\sigma }^2(X^{T}X{)}^{-1}$(σ$\sigma$は残差が正規分布に従うとした場合の標準偏差)ということから、推定量が非常に不安定な値を持つことを意味します。

Ridge回帰 (Hoerl and Kennard, 1970)

多重共線性の問題を回避する方法の一つが、正則化法です。Ridge回帰の場合は、L2$L_2$ノルムを加えたものを目的関数とします。

S(β⃗ )=(y⃗ −Xβ⃗ )T(y⃗ −Xβ⃗ )+λβ⃗ Tβ⃗ (4)$$\begin{array}{}\text{(4)}& S\left(\overrightarrow{\beta }\right)={(\overrightarrow{y}-X\overrightarrow{\beta })}^T(\overrightarrow{y}-X\overrightarrow{\beta })+\lambda {\overrightarrow{\beta }}^T\overrightarrow{\beta }\end{array}$$

λ (λ≥0)$\lambda (\lambda \ge 0)$は正則化の強さを制御する（ハイパー）パラメータです。通常の線形回帰の同様に、解を求めると以下のようになります。

β^⃗ Ridge=(XTX+λ1⃗ )−1XTy⃗ (5)$$\begin{array}{}\text{(5)}& {\overrightarrow{\hat{\beta }}}_{\text{Ridge}}={(X^{T}X+\lambda \overrightarrow{1})}^{-1}{X}^T\overrightarrow{y}\end{array}$$

1⃗ $\overrightarrow{1}$は単位行列です。この単位行列（のλ$\lambda$倍）の行列の寄与のお陰で、(XTX)$(X^{T}X)$ではほとんど同じ列だったものがそうではなくなり、特異ではなくなります。そのため、多重共線性がある状況でも推定が可能になります。線形代数の言葉で言うと、正則化 = 正則行列化((正則行列は、逆行列が存在するような正方行列のことです。行列式が0の場合、逆行列が存在しないのですが、正則化項を導入することで、逆行列が存在するような行列へと補正をかけています。)というわけです。

ただし、λ$\lambda$の値が大きすぎると、推定量の共分散β^⃗ =σ2(XTX)−1$\overrightarrow{\hat{\beta }}={\sigma }^2(X^{T}X{)}^{-1}$(σ$\sigma$は小さくなりますが、バイアス((XTX)/λ+1⃗ )−1β⃗ $((X^{T}X)/\lambda +\overrightarrow{1}{)}^{-1}\overrightarrow{\beta }$が大きくなるため、適切なλ$\lambda$の調整は必要です。

Lasso回帰 (Tibshirani, 1996)

Ridge回帰によって、多重共線性の問題が回避できることがわかりました。しかし、Ridge回帰には、変数選択ができないという課題があります。式(5)から、大きなλ$\lambda$をとることにより回帰係数全体が縮小される傾向があることがわかります。しかしながら、どれかの係数が完全に0になるといった性質は持っておらず、特徴量の数が多いときに解釈のし易い（特定の変数のみが寄与する）モデルが得ることは難しいという問題があります。

一方、変数選択はAICなどのモデル選択規準を用いて変数の追加・削除を行い最適な変数のセットを見つける問題です。こちらは、解釈性の高いモデル構築をできる可能性がありますが、多重共線性の問題がある場合は必ずしも安定な推定量が得られるわけではありません。

Lasso回帰は、Ridge回帰と変数選択の良いとこ取りをしたようなモデルになっています。この場合、L1$L_1$ノルムを加えたものを目的関数とします。

S(β⃗ )=(y⃗ −Xβ⃗ )T(y⃗ −Xβ⃗ )+λ∑j=1p|βj|(6)$$\begin{array}{}\text{(6)}& S\left(\overrightarrow{\beta }\right)={(\overrightarrow{y}-X\overrightarrow{\beta })}^T(\overrightarrow{y}-X\overrightarrow{\beta })+\lambda \sum _{j=1}^p\left|{\beta }_j\right|\end{array}$$

式(6)は原点で不連続なため、通常の線形回帰のように厳密に解を求めることはできません。数値的に解く必要があります。数値計算の方法はいくつか種類があり、よく用いられる手法に座標降下法と呼ばれるものがあります。こういった何らかのアルゴリズムを用いて推定した結果をβ^⃗ LASSO${\overrightarrow{\hat{\beta }}}_{\text{LASSO}}$とします。

β^⃗ LASSO=argminβ[(y⃗ −Xβ⃗ )T(y⃗ −Xβ⃗ )+λ∑j=1p|βj|](7)$$\begin{array}{}\text{(7)}& {\overrightarrow{\hat{\beta }}}_{\text{LASSO}}={\text{argmin}}_{\beta }[{(\overrightarrow{y}-X\overrightarrow{\beta })}^T(\overrightarrow{y}-X\overrightarrow{\beta })+\lambda \sum _{j=1}^p\left|{\beta }_j\right|]\end{array}$$

Lasso回帰は、パラメータ推定と変数選択を同時に行うことができ、スパースなモデルを構築することが可能です。特徴量が2つの場合の直感的なイメージ図は以下になります。正則化法は、制約条件付きの2乗誤差の最小化問題に置き換えることができます。ピンクの領域は、正則化項の大きさを固定した場合に許されるパラメータの領域となっています。最小2乗法による解がこのピンクの領域よりも外側にある場合、制約条件込みで目的関数を最小にするのは、最小2乗法による解を中心とした等高線とピンクの領域が交わる部分（図の星印）となります。Ridge回帰では、制約条件が円状であるため、回帰係数が全体的に小さくなる領域が解になりやすく、Lasso回帰では、尖っている部分（微分不可能な点で、特定の回帰係数の大きさが0となっている）が解になりやすくなります。そのため、Lasso回帰ではいくつかの係数が0となる解が推定されやすくなっています。

Lasso回帰の欠点

Lasso回帰は、スパースなモデルを構築するためには非常に有用ですが、一方で以下のような欠点もあります。

• データ数n$n$, 特徴量の数をp$p$とした場合に、p>n$p>n$のときは高々n$n$個の変数までしか選択できない
• 相関の高い変数群がある場合、Lassoはその中の1つしか変数として選択できない

一つ目は、本当に重要な変数（0でないと推定したい変数）がn$n$よりも多い場合に、重要な変数を選択しきれない、いう意味で問題となります。2つめは、Grouping効果を考慮できないということを意味しています。遺伝子データのような場合、事前知識から似ている特徴量が複数あり、しかもそれらを変数として選択したい状況もあるかと思います。このような場合、Lasso回帰では選択される変数が似ている特徴量の中のどれか一つとなってしまいます。そうすると、モデル作成者が達成したい目的はLasso回帰では達成できなくなります。

Elastic Net回帰 (Zou and Hastie, 2005)

上記のようなLasso回帰の欠点を解決すべく考案されたものが、Elastic Net回帰です。Elastic Net回帰は、正則化項として、L1$L_1$ノルムとL2$L_2$ノルムの和を考えます。

S(β⃗ )=(y⃗ −Xβ⃗ )T(y⃗ −Xβ⃗ )+λ∑j=1p{α|βj|+(1−α)β2j}(8)$$\begin{array}{}\text{(8)}& S\left(\overrightarrow{\beta }\right)={(\overrightarrow{y}-X\overrightarrow{\beta })}^T(\overrightarrow{y}-X\overrightarrow{\beta })+\lambda \sum _{j=1}^p\{\alpha \left|{\beta }_j\right|+(1-\alpha ){\beta }_j^2\}\end{array}$$

λ (λ≥0)$\lambda (\lambda \ge 0)$は正則化項の大きさを制御するハイパーパラメータで、α (0≤α≤1)$\alpha (0\le \alpha \le 1)$はL1$L_1$ノルムとL2$L_2$ノルムの相対的な寄与を調整するハイパーパラメータになります。α=0$\alpha =0$の場合がRidge、α=1$\alpha =1$の場合がLasso回帰に帰着します。

以下、Lasso回帰の欠点をElastic Net回帰がどのように解決しているかを見ていきます。

データ数n$n$, 特徴量の数をp$p$とした場合に、p>n$p>n$のときは高々n$n$個の変数までしか選択できない

こちらは、Elastic Net回帰の場合の推定量の計算を具体的に導出することでわかってきます。今、式(8)に登場する量たちを

X∗=(1+(1−α)λ)(X(1−α)λ−−−−−−−√1⃗ p),   y⃗ ∗=(y⃗ 0⃗ p), γ=αλ/1+(1−α)λ−−−−−−−−−−√,   β⃗ ∗=1+(1−α)λ−−−−−−−−−−√β⃗ (9)$$\begin{gathered} \begin{array}{}\text{(9)}& X^{\ast }=(1+(1-\alpha )\lambda )\left(\begin{array}{c}X\\ \sqrt{(1-\alpha )\lambda }{\overrightarrow{1}}_p\end{array}\right),{\overrightarrow{y}}^{\ast }=\left(\begin{array}{c}\overrightarrow{y}\\ {\overrightarrow{0}}_p\end{array}\right), \\ \gamma =\alpha \lambda /\sqrt{1+(1-\alpha )\lambda },{\overrightarrow{\beta }}^{\ast }=\sqrt{1+(1-\alpha )\lambda }\overrightarrow{\beta }\end{array} \end{gathered}$$

とおくと、式(8)は

S(β⃗ ∗)=(y⃗ ∗−X∗β⃗ ∗)T(y⃗ ∗−X∗β⃗ ∗)+γ∑j=1p|βj|(10)$$\begin{array}{}\text{(10)}& S\left({\overrightarrow{\beta }}^{\ast }\right)={({\overrightarrow{y}}^{\ast }-X^{\ast }{\overrightarrow{\beta }}^{\ast })}^T({\overrightarrow{y}}^{\ast }-X^{\ast }{\overrightarrow{\beta }}^{\ast })+\gamma \sum _{j=1}^p|{\beta }_j|\end{array}$$

とLasso回帰の形で書き直せます。そのため、式(10)の問題はLasso回帰の問題として解くことができ、その解をβ^⃗ ∗${\overrightarrow{\hat{\beta }}}^{\ast }$とすると、Elasitc Net回帰に関する推定量は

β^⃗ Elastic Net=11+(1−α)λ−−−−−−−−−−√β^⃗ ∗(11)$$\begin{array}{}\text{(11)}& {\overrightarrow{\hat{\beta }}}_{\text{Elastic Net}}=\frac{1}{\sqrt{1+(1-\alpha )\lambda }}{\overrightarrow{\hat{\beta }}}^{\ast }\end{array}$$

と求めることができます。今、X∗$X^{\ast }$は(n+p)×p$(n+p)\times p$行列となっているため、この行列のランクはp$p$です。そのため、式(11)の推定量として、最高でp$p$の成分が0でないと推定されます。これは、Lasso回帰のp>n$p>n$のときは高々n$n$個の変数までしか選択できない（0でない回帰係数の最大数はn$n$）の欠点を解消するものになっています。

ここで、X∗$X^{\ast }$というX$X$よりもランクが大きくなりうる行列を定義していますが、これはRidge回帰の性質である疑似データ生成がもとになっています。つまるところ、Lasso回帰の欠点をRidge回帰の性質で補うという構造になっています。

相関の高い変数群がある場合、Lassoはその中の1つしか変数として選択できない

こちらは、原論文に証明があります。今β^i${\hat{\beta }}_i$をi$i$番目のElastic Net回帰の推定量とし、ある(i,j)$(i,j)$に関してβ^i×β^j>0${\hat{\beta }}_i\times {\hat{\beta }}_j>0$が成立しているとします。これは、これらの係数の符号が一致している状況になります。このとき、以下のような量を定義します。

D(α,λ)=1∑ni=1|yi|∣∣β^i−β^j∣∣(12)$$\begin{array}{}\text{(12)}& D(\alpha ,\lambda )=\frac{1}{\sum _{i=1}^n|y_i|}|{\hat{\beta }}_i-{\hat{\beta }}_j|\end{array}$$

詳しい導出は省略しますが、今考えている回帰係数に対応する特徴量間の相関係数をρ$\rho$とすると、式(12)は以下のように評価することができます。

D(α,λ)≤1(1−α)λ2(1−ρ)−−−−−−−√(13)$$\begin{array}{}\text{(13)}& D(\alpha ,\lambda )\le \frac{1}{(1-\alpha )\lambda }\sqrt{2(1-\rho )}\end{array}$$

これは、ρ=1$\rho =1$（相関係数が1）の場合にβ^i=β^j${\hat{\beta }}_i={\hat{\beta }}_j$となることを意味しています。よって、Lasso回帰の場合と異なり、Elastic Net回帰では、相関の高い係数はどちらも0になるか、0にならないなら同じ値になる、といったGrouping効果を持つことがわかりました。

glmnetによる実験

以上のことから、それぞれの正則化法のおおざっぱな特徴は以下のようになります。

• Ridge回帰：回帰係数全体が縮小推定される
• Lasso回帰：相関の高い特徴量群がある場合、そのうちの一つが選択される（対応する回帰係数以外が0と推定される）
• Elasitic Net回帰：相関の高い特徴量は、回帰係数が両方0になる、もしくは両方同じ値になる

これらの特徴を、ある程度恣意的な状況でモデルを構築して確かめてみることにします。使用データはボストン近郊の住宅価格のデータで、議論を簡単にするために変数をある程度絞って議論します。絞る材料とするために、いったんglmでモデルを構築し、stepwise法で変数選択をしてみます。

# glmでモデル作成
mod.full <- glm(medv ~., data = Boston)
# AICを基準にstepwise
mod.step <- stepAIC(mod.full, direction = "both", trace = 0)
# summaryを出力
summary(mod.step)

# 係数の部分だけ記載
Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  36.341145   5.067492   7.171 2.73e-12 ***
crim         -0.108413   0.032779  -3.307 0.001010 ** 
zn            0.045845   0.013523   3.390 0.000754 ***
chas          2.718716   0.854240   3.183 0.001551 ** 
nox         -17.376023   3.535243  -4.915 1.21e-06 ***
rm            3.801579   0.406316   9.356  < 2e-16 ***
dis          -1.492711   0.185731  -8.037 6.84e-15 ***
rad           0.299608   0.063402   4.726 3.00e-06 ***
tax          -0.011778   0.003372  -3.493 0.000521 ***
ptratio      -0.946525   0.129066  -7.334 9.24e-13 ***
black         0.009291   0.002674   3.475 0.000557 ***
lstat        -0.522553   0.047424 -11.019  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

結果を見ると、相対的にはrm, dis, ptratio, lstatの係数が有意であるようです。というわけで、これら4つの変数に絞り、さらにrmに対しては全く同じ変数のコピーを新たに加えたデータを作成します。

# rmのコピー
rm.dummy <- data.frame(rm_dummy = Boston$rm)
# new data.frame
Boston.new <- cbind(Boston[, c(6, 8, 11, 13, 14)], rm.dummy)

以下では、rmのコピー変数の回帰係数の解パスをモデルごとを見ていくことで、それぞれの正則化法の違いを見ていくことにします。その際に、glmnetに関する便利な補助ツールが幾つかあることが最近わかったので、活用していきます。具体的には、glmnetでもformula的に書くことができるようにするglmnetUtilsと、解パスをggplot2の形式で可視化してくれるggfortulyというものです。こちらに関しては、hoxo_mさんが最近紹介してくださりました。
qiita.com

Ridge回帰

先程恣意的に作成したデータで、Ridge回帰をし、解パスを見てみましょう。

# 必要なライブラリ
if(!require(glmnet)){
  install.packages("glmnet", quiet = TRUE)
}
require(glmnet)
require(glmnetUtils)
if(!require(githubinstall)){
  install.packages("githubinstall", quiet = TRUE) 
}
require(githubinstall)
if(!require("ggfortify")){
  githubinstall("ggfortify")
}
require(ggfortify)

# Ridge回帰（glmnetUtilsを併用）
Ridge <- glmnet(medv ~ ., data = Boston.new, alpha = 0)
# ggfortifyで可視化
autoplot(Ridge, xvar = "lambda")

（途中0を横切っているのもありますが…）正則化パラメータλ$\lambda$を大きくしていくと、係数の大きさが全体的に0に縮小していく傾向が見えます。また、rmとそのコピーrm_dummyはほぼほぼ同じような動きをしています。

Lasso回帰

Lasso回帰の場合にも、同じように解パスを見てみます。

# Lasso回帰（glmnetUtilsを併用）
Lasso <- glmnet(medv ~ ., data = Boston.new)
# ggfortifyで可視化
autoplot(Lasso, xvar = "lambda")

Ridge回帰の場合と比較し、rmとrm_dummyの振る舞いが明らかに違います。rm_dummyの回帰係数が常に0と推定されています。これは、Lasso回帰の「相関の高い特徴量群がある場合、そのうちの一つが選択される」という性質に対応しています。

Elastic Net回帰

# Elastic Net回帰（glmnetUtilsを併用）
ElasticNet <- glmnet(medv ~ ., data = Boston.new, alpha = 0.5)
# ggfortifyで可視化
autoplot(ElasticNet, xvar = "lambda")

Lasso回帰の欠点であった、Grouping効果が反映されています。少し見づらいですが、rmとrm_dummyの係数が0でなくなるタイミングが同じになっており、いったん0でないと係数された場合は、両者の回帰係数が似たような値になっています。

精度について

今回は主題ではないのであえて触れませんでしたが、実際に精度検証をしてみると、今回のデータではglm w/ stepwise > Elastic Net, Lasso > Ridgeの順で決定係数が大きかったです。例によってテストデータを含めた検証をしっかりやっているわけではないので、ここは参考程度で。実務の際はちゃんと汎化性能・解釈性などを踏まえて議論します。

まとめ

今回は、ある程度理論的な性質をもとに代表的な正則化法をまとめてみました。かなり恣意的な状況のデータではありましたが、相関が高い特徴量に対する実際の挙動の違いも確認しました。精度追求が要件の場合には今回の議論はあまり意味をなさないと思いますが、一方で、業務知識などである程度取り入れたい考えがあるときには、それに応じてモデルを構築するヒントの一助となるのではないかと思います。

*1:最尤法の枠組みで、正規分布に従う残差を仮定した場合でも、同じ目的関数となります。