(cache)Rでスパースモデリング：Adaptive Lasso - データサイエンティスト(仮)

導入

スパース推定の代表的な手法として、Lassoがあります。様々なシーンで活用されているLassoですが、Lassoは変数選択の一致性が保証されないという欠点があります。Adaptive Lassoは、その欠点を補う形で提唱されている手法となっています。こちらは、ある条件のもとで変数選択の一致性が保証*1されており、数理統計学的により好ましい性質を持っています。

このAdaptive Lassoですが、Rでは{glmnet}以外のパッケージを使わないと簡単にできないとかなりの期間勘違いをしてました。そんな中、以下の記事を最近見かけまして、冷静に考えたら{glmnet}でも表現できるよなあと猛省した次第です。

RPubs - Adaptive LASSO Examples

以上の経緯から、挙動を確かめておこうという考えのもと、メモがてらAdaptive Lassoの紹介をしようと思います。

補足

上記の一致性の議論は、特徴量の数を固定して、データ数に対して無限大の極限を取った場合です。特徴量の数をデータ数よりも常に大きくとるようにした場合の極限の振る舞いも議論されており、例えば
http://www.keihirose.com/material/392-399_hirose.pdf
に載っております。

使用するデータセット

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

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

# 目的変数、説明変数
y <- as.matrix(Boston[, 14])
X <- as.matrix(Boston[, -14])

目的変数はmedvで、説明変数はそれ以外です。また、{glmnet}にデータを投入する際は

目的変数と説明変数を分ける
それぞれmatrixにする
変数の中にfactorがあってはならない

といった制約があるので、それらを意識しています。

また、今回は予測に重きを置いているわけではなく、回帰係数の振る舞いに興味があるので、予測のために学習・検証用データに分けることはせず、全データを用いて推定を行います。

Lasso

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

\vec{y} = X \vec{β} + \vec{ϵ}

y, \vec{β}, \vec{ϵ}

はそれぞれ目的変数、回帰係数（推定するパラメータ）、残差で、

X

は計画行列（説明変数をまとめた行列）です。データ数

n

、特徴量の数

p

とすると、

y, \vec{ϵ}

は

n

次元ベクトル、

\vec{β}

は

p + 1

次元ベクトル（切片項も含む）、

X

は

n \times (p + 1)

行列となります。

パラメータ推定をする際の目的関数は、二乗誤差を仮定します。Lasso(Least Absolute Shrinkage and Selection Operator)では、正則化項として $L_{1}$ ノルムを目的関数に加えます。

L (\vec{β}) = {∥ \vec{y} - X \vec{β} ∥}^{2} + λ ∥ \vec{β} ∥_{1}

λ

は正則化項（もしくは罰則項）と呼ばれ、もとの二乗誤差に対して正則化の強さを制御するようなハイパーパラメータです。

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

{\vec{β}}_{L a s s o} = {argmin}_{β} [{∥ \vec{y} - X \vec{β} ∥}^{2} + λ ∥ \vec{β} ∥_{1}]

Lassoの特徴は、パラメータ推定と変数選択が同時にできることです。変数選択は、例えば相関の高い２つの変数があったとき、片方の変数の回帰係数が0と推定され、もう片方は0でない値となる、といったことが起きます。こういった状況のもと、推定結果が0でなかったものが変数として選択された、と考えます。この性質は、LassoがAICやBICといったモデル選択規準を用いてステップワイズ法を考える、といった変数選択のメカニズムが背後にあることに起因しています。こういった状況は、正則化項が離散的な関数になっていると読み取れるのですが、この項をうまく凸関数になるように連続極限をとった場合がLassoになります*2。

他の正則化項として代表的なものとしては、 $L_{2}$ ノルムを加えるRidge、 $L_{1}$ ノルムと $L_{2}$ ノルムの線形結合を加えるElastic netなどがあります。これら正則化やアルゴリズムのより詳しい話は、例えば以下の記事などを参考としてください。
qiita.com

追記

argminの部分、投稿時はargmaxとタイポしてました。ご指摘いただいたynakahashiさんに感謝です！

Lassoの解パス

ハイパーパラメータである正則化項 $λ$ を変化させると、推定結果も変わります。この変化の移り変わりを可視化したものを解パス（Solution Path）といいます。ここでは、ボストン近郊の住宅価格のデータに対してLassoを適用して、解パスがどうなるかを見てみます。

# モデリング
Lasso <- glmnet(x = X, y = y,
                alpha = 1)
# 正則化項を変化させたときの回帰係数の推移（解パス）
plt_res(Lasso)

f:id:tekenuko:20171102005446p:plain

ここで、plt_resはglmnet関数の出力をうまくggplot2で可視化する自作関数です。

# Library
if(!require(ggplot2)){
  install.packages("ggplot2", quiet = TRUE)
}
require(ggplot2)
if(!require(tidyr)){
  install.packages("tidyr", quiet = TRUE)
}
require(tidyr)

plt_res <- function(dgMatrix){

  # Convert to data.frame
  df_X <- as.data.frame(t(as.matrix(dgMatrix$beta)))
  df_loglambda <- dgMatrix$lambda
  df <- cbind(df_loglambda, df_X)
  # melt
  df_melt <- df %>%
    gather(key = df_loglambda, value = value)
  colnames(df_melt) <- c("lambda", "coefficient", "beta")
  # plot
  plt <- ggplot(df_melt) + 
    geom_line(aes(x = lambda, y = beta, group = coefficient, col = coefficient)) + 
    scale_x_log10()
  
  return(plt)
  
}

図を見ると、正則化項の寄与が小さい（ $λ$ が小さい）ときは、0でない値に推定される回帰係数が多く、正則化項の寄与が大きくなると0と推定される係数が増えていくといった振る舞いになっています。このように、正則化項の寄与を変化させると、推定される係数も変化していきます。

回帰係数

上で解パスを見ましたが、ハイパーパラメータはどの値を取るのがよいでしょうか。上記の解パスの各 $λ$ ごとにRMSEなどが評価できるはずなので、その評価指標が最も良い点を選ぶ、というのも一つです。ここでは、より汎用的に評価するために、cross-validationを考えます。{glmnet}では、cv.glmnet関数を用いると簡潔に計算できます。

fit.Lasso.cv <- cv.glmnet(x = X, y = y,
                      nfolds = 10,
                      alpha = 1, 
                      standardize = TRUE)

nfoldはデータを10分割し、そのうちの一つをテストデータとする、といったサンプルを10パターン生成することを意味しています。alphaは $L_{1}$ ノルムと $L_{2}$ ノルムの相対的な寄与の大きさをコントロールしているパラメータでalpha=1のときがLasso、alpha=0のときがRidgeとなります。standardizeはモデリングの際にデータを標準化するかどうかのオプションです。

この結果、plot関数を用いると以下のように可視化されます。

plot(fit.Lasso.cv)

f:id:tekenuko:20171102010846p:plain
これは、 $λ$ とRMSEの関係をプロットしたもので、エラーバーはデータ分割の仕方が起源として発生しています。この図を見ると、 $λ$ は非常に小さい値が好まれるようです。

最もRMSEが小さい $λ$ の場合の推定結果は、以下になります。

coef(fit.Lasso.cv, s = fit.Lasso.cv$lambda.min)

# 出力
14 x 1 sparse Matrix of class "dgCMatrix"
                        1
(Intercept)  34.880894915
crim         -0.100714832
zn            0.042486737
indus         .          
chas          2.693903097
nox         -16.562664331
rm            3.851646315
age           .          
dis          -1.419168850
rad           0.263725830
tax          -0.010286456
ptratio      -0.933927773
black         0.009089735
lstat        -0.522521473

indusとageの係数が0と推定されています。

Adaptive Lasso

Adaptive Lassoは、 $L_{1}$ ノルムを目的関数に加える点はLassoと同じですが、重みが付く形になります。

L (\vec{β}) = {∥ \vec{y} - X \vec{β} ∥}^{2} + λ \sum_{i = 1}^{p} ω_{i} | β_{i} |

ここで、重みベクトル

\vec{ω} = (ω_{1}, \dots, ω_{p})

は以下になります。

\vec{ω} = 1 / {| \vec{\hat{β}} |}^{γ}

\vec{\hat{β}}

はある一致推定量で、一致推定量の各成分の絶対値の

γ

の逆数が

\vec{ω}

の各成分となっています。また

γ > 0

です。

一致推定量は、真の値との差が $O (1 / \sqrt{n})$ 以下程度の精度が担保されるものなら何でもよいです。多くの場合、OLSの結果やRidgeの結果を使うようです。今回は、Ridge回帰の結果を用いることにします。

Ridge回帰により何らかの一致推定量を得る

必ずしもcross-validationまでやる必要はありませんが、10-fold cross-validationにより、最もRMSEが小さくなる場合の回帰係数を取得します。

# Ridge回帰
fit.ridge.cv <- cv.glmnet(x = X, y = y,
                          type.measure = "mse",
                          alpha = 0, 
                          nfolds = 10, 
                          standardize = TRUE)
# 最もRMSEが小さくなる正則化項の場合の回帰係数
best_ridge_coef <- as.numeric(coef(fit.ridge.cv, s = fit.ridge.cv$lambda.min))[-1]

Adaptive Lassoの解パス

上で取得した係数を反映します。Adaptive Lassoは、glmnet関数のオプションpenalty.factorに1 / abs(best_ridge_coef))を代入することで実現できます。このオプションを認識していなかったので、glmnetでは簡単にAdaptive Lassoを書けないと思っていました。。解パスを表示すると、以下のようになります。

# モデリング
AdaLasso <- glmnet(x = X, y = y,
                  alpha = 1,
                  penalty.factor = 1 / abs(best_ridge_coef))
# 解パスを表示
plt_res(AdaLasso)

f:id:tekenuko:20171102013024p:plain

線が上下している係数もあるので、実はあまり推定がうまくいっていないのかも、という懸念はありますが、Lassoの場合と比較して、だいぶ振る舞いが変わりました。どのタイミングでどの係数が0になるか、線の形状、色々変化しています。動く $λ$ の範囲は、実質的にRidgeの結果で割りこんだものを考えているので変化するのは理解できます。他の変化の仕方は明確に理由を説明するのは難しそうです。

回帰係数

Adaptive Lassoの場合も係数を表示しておきましょう。cross-validationした場合の図は省略します。

fit.AdaLasso.CV <- cv.glmnet(x = X, y = y,
                        type.measure = "mse",
                        nfolds = 10,
                        alpha = 1,
                        penalty.factor = 1 / abs(best_ridge_coef))
coef(fit.AdaLasso.CV, s = fit.AdaLasso.CV$lambda.min)

# 出力
14 x 1 sparse Matrix of class "dgCMatrix"
                        1
(Intercept)  38.501580556
crim         -0.098114529
zn            0.022165389
indus        -0.016969903
chas          3.129489627
nox         -20.469852957
rm            3.948669580
age           .          
dis          -1.345793859
rad           0.097441871
tax           .          
ptratio      -1.020675606
black         0.001601036
lstat        -0.546279678

0となる係数の組み合わせ、および係数の値が変化しています。このように、回帰係数のパターンもLassoと比較して変わるようです。

精度比較

精度が全てというわけではありませんが、LassoとAdaptive Lassoの両者で精度比較を行っておきます。評価指標はRMSEと決定係数 $R^{2}$ ですが、上述の通り、今回は全データを推定に使ったので、学習データに関するものの評価になっています*3。RMSEは{Metrics}のメソッドで計算し、決定係数 $R^{2}$ は予め関数を定義しておきます。

# 決定係数
r_squared <- function(y, ypred) {
  ybar <- mean(y)
  ## Total SS
  ss_tot <- sum((y - ybar)^2)
  ## Residual SS
  ss_res <- sum((y - ypred)^2)
  ## R^2 = 1 - ss_res/ ss_tot
  1 - (ss_res / ss_tot)
}

# LassoとAdaptive LassoそれぞれRMSEとR^2を計算
c("RMSE(Lasso) : ", round(rmse(y, pred_Lasso), 4))
c("R^2(Lasso) : ", round(r_squared(y, pred_Lasso), 4))
c("RMSE(Lasso) : ", round(rmse(y, pred_AdaLasso), 4))
c("R^2(Adaptive Lasso) : ", round(r_squared(y, pred_AdaLasso), 4))

# 出力
"RMSE(Lasso) : " "5.0299"   
"R^2(Lasso) : " "0.7003" 

"RMSE(Adaptive Lasso) : " "5.0435"
"R^2(Adaptive Lasso) : " "0.6987"

Adaptive Lassoの方がRMSEが大きく、決定係数が小さいという結果になってしまいました。変数選択をより正確にしたことの代償なのか…。今回はテストデータも含めて検証をしておりませんし、何よりデータに依存する話かもしれないので、ここの比較結果はあくまで参考程度にお願いします。

まとめ

今回は、数理統計学的に性質のよいとされるAdaptive Lassoのglmnetでの実装例を紹介しました。今まで壮大に勘違いをしていたのですが、案外シンプルにglmnetで書けることがわかったので、変数選択の一致性を気にする場合は積極的に活用していこうと思います。

*1:他にも推定量が漸近正規性をもつという性質も合わせて、オラクルプロパティを持つ、といいます。

*2:このあたりの議論は、統計数理研究所の公開講座「スパース推定」でされております。私はこの話を聴いたとき非常に感銘をうけました。

*3:実務接続を意識して比較する場合は、予測に関する評価もきちんと行わなければなりません。