Fitbitデータでスパース推定に入門してみる② 〜lasso, elastic net, SCAD, adaptive lassoの概要とRでの数値例〜
はじめに
前回の記事では、スパース推定の概要と代表的な手法であるlassoについてまとめ、fitbitデータに適用することで、睡眠効率に影響を与える要因を探索しました。
今回のエントリでは、lasso以外の様々なスパース正則化法を紹介し、前回と同様にFitbitデータに適用することで、それぞれの手法を比較したいと思います。なお、今回のエントリは、以下の書籍の第１・２章を自分なりに簡単にまとめた*1だけなので、丁寧に学びたい方は、書籍を読んでいただければと思います。
スパース推定法による統計モデリング (統計学One Point)
- 作者: 川野秀一,松井秀俊,廣瀬慧
- 出版社/メーカー: 共立出版
- 発売日: 2018/03/08
- メディア: 単行本
- この商品を含むブログを見る
lassoの概要と弱点
概要
lassoは以下のように、最小二乗推定量を算出する式の後ろに、正則化項にL1ノルムを考えたものです。
Lassoはパラメータを0と推定することで、変数選択をすることが可能です。例えば相関が高い２つの説明変数があったとき、片方の説明変数のパラメータが0と推定されて、その変数はモデルから除外されます。
弱点
しかし、lassoは２つの問題があることが知られています。
- 相関の高い２つの説明変数が目的変数に関係している場合、一方の説明変数しかモデルに含まれない
- サンプルサイズより説明変数の数の方が大きい場合（所謂 n < p の場合）、高々サンプルサイズ分の変数までしか選択されない
これらの問題を解決するという観点でelastic netが提案されました。さらにオラクル性質（oracle property）*2を持たせるという観点で、正則化項として非凸正則化項を使うSCADやadaptive lassoが提案されています。以降ではこれらの概要を説明してから、Rで数値計算してみます。
lasso正則化項の拡張
elastic net
elastic netは上記２つの問題を解決させるために提案された手法です。最小２乗推定量を算出する式の後ろに、下記のような正則化項を加えたものを最小化することで、パラメータを推定します。
これはL1ノルムとL2ノルムを合わせた形になっており、α=1の時はlasso, α=0の時はリッジとなります。このような正則化項を取ることで、上記の問題が解決されます。つまり、相関が高い変数を両方採用したい場合や、サンプルサイズ以上の変数を採用したい場合はelastic netを使用するのが良いということになります。
SCAD
以下のような正則化項をもつスパース推定法をSCADと言います。
λは正則化パラメータ、a(＞2)は調整パラメータです*3。SCADはオラクル性質を満たしており、漸近正規性や変数選択の一致性を有します。
adaptive lasso
SCADはオラクル性質を満たしますが、複雑な正則化項の形をしており、lassoの形とは大きく異なります。そこで、lassoのL1ノルムの部分に適当な重みwを置くadaptive lassoが提案されています。
これはある正則条件の下でオラクル性質を満たします。重みwは既知でなければなりませんが、例えば最小二乗推定量βを使って、w=1/βと推定して用います。 wが大きい（βが小さい）とモデルに含む必要がないため回帰係数を0に推定する方向に持っていき、逆にwが小さい（βが大きい）と最小２乗推定量になるべく近い推定量を採用することを意味します。上記推定式は適切な式変形をすればlassoと同じアルゴリズムで推定値を得ることができますが、これについては以下のRでの数値計算例の時に簡単に説明します。
数値例
前回と同様にfitbitデータと自作の生活行動データを使用して、睡眠効率を目的変数としました。最小２乗推定、lasso、elastic net、SCAD、adaptive lassoそれぞれのアプローチで数値計算し結果を比較します。
データ*4
> glimpse(df) Observations: 43 Variables: 9 $ totalTimeInBed <int> 423, 357, 590, 426, 463, 426, 440, 417, 381, 588... $ efficiency <dbl> 0.9290780, 0.9607843, 0.9067797, 0.9436620, 0.95... $ Work <int> 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, ... $ Drink <int> 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, ... $ Glass <int> 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, ... $ Bath <int> 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, ... $ Stretch <int> 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, ... $ weekday <int> 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, ... $ step8000 <int> 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, ...
標準化します。
y <- unlist(select(df, efficiency)) y <- y - mean(y) x <- scale(as.matrix(select(df, -efficiency)))
本来であれば、しっかりと可視化をしてデータの分布などを確認すべきですが割愛します。
lasso
lassoは前回の記事で紹介した通りです。
library("glmnet") fit <- glmnet(x, y, standardize = FALSE) autoplot(fit) + theme_pander() + scale_color_pander() # solution path set.seed(0) cvfit <- cv.glmnet(x, y, standardize = FALSE, nfolds = 10) coef(cvfit, s = "lambda.min") # 回帰係数
elastic net
glmnetパッケージのαパラメータに0~1の数値を代入すればよく（ここではα=0.5とした）、それ以外はlassoと全く同じです。
fit <- glmnet(x, y, alpha = 0.5, standardize = FALSE) autoplot(fit) + theme_pander() + scale_color_pander() # solution path set.seed(0) cvfit <- cv.glmnet(x, y, alpha = 0.5, standardize = FALSE, nfolds = 10) coef(cvfit, s = "lambda.min") # 回帰係数
lassoとほぼ同じ（と言うか全く同じ？）ですね。正直よくわかりません（やっぱりデータがよくないのかな…）。
SCAD
ncvregパッケージを使用します。
library("ncvreg") set.seed(0) fit <- ncvreg(x, y, family="gaussian", penalty = "SCAD") plot(fit) # solution path
影の部分はnonconvex region（非凸領域）です*5。 デフォルトのグラフでも綺麗ですが、各線がどの変数に紐づいているか不明なので、自作します。
lambda_df <- as.data.frame(fit$beta) # パラメータの抽出、data frame化 lambda_df$variable <- rownames(lambda_df) # tidy化 df <- lambda_df %>% gather(key="lambda", value="beta", -variable) %>% filter(variable != "(Intercept)") %>% # 切片変数を除外 mutate(lambda = as.double(lambda)) # 可視化 ggplot(df, aes(lambda, beta, color=variable)) + geom_line() + theme_pander() + scale_color_pander() + scale_x_reverse()
lassoやelastic netのsolution pathとは大きく異なりますが、出てくる変数（パラメータが非0となる変数）の順番は同じっぽいです。
cross validationによりパラメータλを選択し、回帰係数を算出します。
set.seed(0) cvfit <- cv.ncvreg(x, y ,family="gaussian", penalty = "SCAD") opt_lambda <- cvfit$lambda.min # 方程式を最小にするλ fit_lambda.min <- ncvreg(x, y, family="gaussian", penalty = "SCAD", lambda = opt_lambda) fit_lambda.min$beta # 回帰回数
adaptive lasso
adaptive lassoもlassoと同じglmnetパッケージを使用することができますが、少し工夫が必要です。以下のように式変形すると、
lassoのアルゴリズムが使用できます。以下、ざっとした計算の流れをまとめます。
- 重みwを最小２乗推定値により算出*6
を算出・・・①
- glmnetによりパラメータ推定
によりパラメータを得ることができます。
W <- diag(1/abs(coef(lm(y~x))[-1])) xstar <- x %*% solve(W) # ①を算出 colnames(xstar) <- colnames(x) fit <- glmnet(xstar, y, standardize = FALSE) autoplot(fit) + theme_pander() + scale_color_pander() # solution path set.seed(0) cvfit <- cv.glmnet(xstar, y, standardize = FALSE, nfolds = 10) beta <- coef(cvfit, s = "lambda.min") solve(W) %*% beta[-1] # 回帰係数
※ を変換した
に対するsolution pathです。参考までにやってみたものの、どのように解釈していいか不明です。
結果の解釈
推定結果をまとめた表が以下になります。
書籍には、
非凸正則化によるスパース推定法では、lassoに比べてより多くの回帰係数が0と推定されている。非0と推定されている部分に注目すると、これらの手法はlassoに比べて、最小２乗推定値に近い値が推定されていることがわかる。
と書いているが、ここではデータが微妙なせいか、よく分かりませんでした。
さいごに
今回のエントリでは、様々なスパース正則化法についての概要をまとめ数値例を示しました。*7
最後に個人的に気になった点、学びたい点をMEMOしておきます。
分からないことだらけですね。まだまだINPUTが足らないようなので、引き続き勉強頑張ります。
なお、本エントリは雰囲気で書いた部分も多いので、間違いがあれば遠慮なくご指摘くださいませ。
mikutaifuku
*1:時間の都合上（難しかったので）、計算アルゴリズムは飛ばしています。
*2:漸近正規性、変数選択の一致性を有する性質のことを言うらしいが、時間の都合上（難しかったので）詳細は書籍に譲ります
*3:aはリスク最小化の観点からa=3.7が用いられているらしいです。cross validationなどで選択することも可能です
*4:データ数少ないですね。最近サボり気味なのでしっかりデータ集めます。
*5:よくわからなかったので余力があればしっかりと調べます
*6:リッジでも問題ないが、リッジ推定を用いた場合はオラクル性質を満たさないらしい
*7:疲れて最後の方は少し適当になってしまった
*8:adaptive lassoでもある条件の下だけらしいが
