今回はせっかくBayesian Lassoを使って遊んでみたのでその記録を残しておきたいと思います。(Bayesian Lassoを扱った記事も少ないように思われたので)
Bayesian Lassoとは
そもそもBayesian Lassoとは何なのでしょうか。
厳密なことであったり、詳しいことはPark & Casella (2008)にあるので気になる方はそちらもご参照ください。ここでは簡単に説明するにとどめます。通常の最小二乗法の回帰やLasso,Ridge回帰の結果に対するinferenceの手法としてbootstrapと並ぶ2大手法といったイメージでしょうか。(違ったら訂正お願いします。)
Park & Casella (2008)によると条件付きラプラス分布
$$
\pi\left(\boldsymbol{\beta} | \sigma^{2}\right)=\prod_{j=1}^{p} \frac{\lambda}{2 \sqrt{\sigma^{2}}} \exp \left(-\frac{\lambda\left|\beta_{j}\right|}{\sqrt{\sigma^{2}}}\right)
$$
を事前分布とし、以下の式
$$
\begin{aligned}
\frac{a}{2} \exp (-a|z|) & \
=\int_{0}^{\infty} \frac{1}{\sqrt{2 \pi s}} \exp \left(-\frac{z^{2}}{2 s}\right) \cdot \frac{a^{2}}{2} \exp \left(-\frac{a^{2} s}{2}\right) d s
\end{aligned}
$$
の正規分布の尺度混合分布に階層化し、事後分布を推定するために以下のモデルが立式できます。
実際にBayesian Lasso で遊んでみる。
理論的なことをざっくりと把握したので(本当にざっくり)、次は実際に遊んでみましょう。習うより慣れろです。
今回は題材としてはRでBayesian Lassoを行うときに用いるblassoの公式ドキュメントを用います。(公式ドキュメント大事！)
ここにあるサンプルコードを1行ずつ追ってみたいと思います。
library(monomvn); library(lars); library(glmnet);
まず、今回必要なライブラリを読み込みたいと思います。インストールしていない場合はインストールしてから実行してください。(install.packageとして直接コマンドに打ち込んでもいいですが、Rstudio)を用いている方だと右下あたりにInstallというところがるので、そちらからインストールしてもよいでしょう。)
data(diabetes)
attach(diabetes)
ここでは今回用いるデータを読み込みます。有名なdiabetesのデータセットを用います。
attach()関数を使うと、そのあとは、どのデータフレームを用いるかをかを指定しなくてよくなるため便利です。なので、ここでattachもしておきましょう。
## Ordinary Least Squares regression
reg.ols <- regress(x, y)
## Lasso regression
reg.las <- regress(x, y, method="lasso")
## Bayesian Lasso regression
reg.blas <- blasso(x, y)
plot(reg.blas, burnin=200)
points(drop(reg.las$b), col=2, pch=20)
points(drop(reg.ols$b), col=3, pch=18)
legend("topleft", c("blasso-map", "lasso", "lsr"),
col=c(2,2,3), pch=c(21,20,18))
次に実際に回帰を行います。そして、その結果得られた各説明変数の係数について図示します。ここでは最小二乗法を用いた回帰、通常のLasso回帰、Bayesian Lasso による回帰を行っております。これを出力した結果が以下です。コード中のburninに見覚えがない人が多いとは思われますが、これはバーニンを表します。これはMCMC法において良く用いられる用語で、MCMC法において初期値の影響が出来ると考えられる最初の方のデータを取り除くことを表します。今回の場合は最初から200番目までのデータは捨てるということを意味しています。boxplotであらわされているのがBayesian Lassoによってでてきたもの点がそのほかの二つの回帰によって得られた推定値となります。
plot(reg.blas, burnin=200, which="m")
ここではblassoによって得られたmの分布を表示します。ここでいうmとはBayesian Lasso回帰を行った際に係数が0でないものの数です。これを見ると7個の係数が0でないと判断された時が一番多いことが分かりますね。下のmchainとは1000回の試行でmがどのように推移していったかを示すものです。
## get the summary
s <- summary(reg.blas, burnin=200)
## calculate the probability that each beta coef != zero
s$bn0
b.1 b.2 b.3 b.4 b.5 b.6
0.25375 0.98875 1.00000 1.00000 0.63125 0.44000
b.7 b.8 b.9 b.10
0.85625 0.47375 1.00000 0.36750
次にBayesian Lassoの結果のsummaryを得たものを用意し、そのうちのbn0を表示します。ここでいうbn0とはコードにおいてコメントアウトされている部分にもありますがそれぞれの係数が0でない確率を表しています。この結果と先ほどのboxplotを見比べてみると(当然ですが)うまく対応していて面白いですね。
ちなみに、bn0だけでなくすべての結果を出力すると以下のようになります。
> s
$call
blasso(X = x, y = y)
$B
[1] 200
$T
[1] 1000
$thin
[1] 10
$coef
mu b.1
Min. :142.3 Min. :-137.1054
1st Qu.:150.3 1st Qu.: 0.0000
Median :152.1 Median : 0.0000
Mean :152.0 Mean : -0.1623
3rd Qu.:153.8 3rd Qu.: 0.0000
Max. :160.4 Max. : 173.6658
b.2 b.3 b.4
Min. :-363.0 Min. :324.7 Min. : 42.98
1st Qu.:-241.9 1st Qu.:486.9 1st Qu.:263.40
Median :-201.1 Median :529.0 Median :306.93
Mean :-198.7 Mean :529.1 Mean :306.84
3rd Qu.:-157.8 3rd Qu.:572.1 3rd Qu.:349.41
Max. : 0.0 Max. :705.8 Max. :514.01
b.5 b.6
Min. :-892.02 Min. :-400.18
1st Qu.:-181.38 1st Qu.: -24.49
Median : -61.41 Median : 0.00
Mean :-106.02 Mean : -10.47
3rd Qu.: 0.00 3rd Qu.: 0.00
Max. : 251.09 Max. : 606.66
b.7 b.8 b.9
Min. :-485.8 Min. :-249.20 Min. :255.9
1st Qu.:-277.5 1st Qu.: 0.00 1st Qu.:451.2
Median :-214.4 Median : 0.00 Median :508.8
Mean :-194.3 Mean : 48.05 Mean :513.4
3rd Qu.:-120.1 3rd Qu.: 77.96 3rd Qu.:568.5
Max. : 109.3 Max. : 516.95 Max. :871.6
b.10
Min. :-97.78
1st Qu.: 0.00
Median : 0.00
Mean : 23.93
3rd Qu.: 35.57
Max. :225.64
$s2
Min. 1st Qu. Median Mean 3rd Qu. Max.
2348 2834 2969 2981 3109 3586
$lambda2
Min. 1st Qu. Median Mean 3rd Qu.
0.008017 0.050188 0.083472 0.100220 0.132015
Max.
0.621306
$tau2i
tau2i.1 tau2i.2
Min. : 0.0030 Min. :0.003606
1st Qu.: 0.0610 1st Qu.:0.027874
Median : 0.1126 Median :0.054310
Mean : 0.6519 Mean :0.096423
3rd Qu.: 0.3412 3rd Qu.:0.107230
Max. :35.3498 Max. :1.630862
NA's :597 NA's :9
tau2i.3 tau2i.4
Min. :0.002158 Min. :0.002433
1st Qu.:0.015728 1st Qu.:0.022765
Median :0.025703 Median :0.040232
Mean :0.031993 Mean :0.057687
3rd Qu.:0.040386 3rd Qu.:0.071626
Max. :0.190404 Max. :0.407075
tau2i.5 tau2i.6
Min. : 0.00286 Min. :0.0032
1st Qu.: 0.02884 1st Qu.:0.0363
Median : 0.06145 Median :0.0854
Mean : 0.24860 Mean :0.2770
3rd Qu.: 0.14471 3rd Qu.:0.2260
Max. :26.01915 Max. :6.4826
NA's :295 NA's :448
tau2i.7 tau2i.8
Min. : 0.00297 Min. : 0.0035
1st Qu.: 0.02769 1st Qu.: 0.0376
Median : 0.05071 Median : 0.0721
Mean : 0.11148 Mean : 0.4795
3rd Qu.: 0.09714 3rd Qu.: 0.1745
Max. :10.73615 Max. :27.7849
NA's :115 NA's :421
tau2i.9 tau2i.10
Min. :0.0005692 Min. : 0.0046
1st Qu.:0.0155462 1st Qu.: 0.0452
Median :0.0254085 Median : 0.1076
Mean :0.0330458 Mean : 1.0034
3rd Qu.:0.0425784 3rd Qu.: 0.2784
Max. :0.1597780 Max. :106.2569
NA's :506
$bn0
b.1 b.2 b.3 b.4 b.5 b.6
0.25375 0.98875 1.00000 1.00000 0.63125 0.44000
b.7 b.8 b.9 b.10
0.85625 0.47375 1.00000 0.36750
$m
Min. 1st Qu. Median Mean 3rd Qu. Max.
5.000 6.000 7.000 7.011 8.000 10.000
attr(,"class")
[1] "summary.blasso"
ここで得られたmと先ほど得られたmの分布を見比べても面白いですね。
せっかくなのでここで出てきた値の意味について重要なものをピックアップして解説していきたいと思います。
TとはMCMC法によって得られたサンプル数の合計を表します。coefはそれぞれの係数について(今回の場合1000回得られる)の代表値を示しています。lambda2は今回はLasso回帰ということもあり、この罰則項の係数を示しています。この値が0ならそれは、通常の最小二乗法による回帰が行われていることを表します。
## summarize s2
plot(reg.blas, burnin=200, which="s2")
s$s2
Min. 1st Qu. Median Mean 3rd Qu. Max.
2348 2834 2969 2981 3109 3586
## summarize lambda2
plot(reg.blas, burnin=200, which="lambda2")
s$lambda2
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.008017 0.050188 0.083472 0.100220 0.132015 0.621306
ここでは先ほど解説したs2,lambda2についてその分布を図にして出力しております。見方は基本的には先ほどのmと同じです。
以上簡単ですがblassoの公式ドキュメントのサンプルコードを1行ずつ追ってみました。
個人的見解ですが、こういった統計で用いられるような概念というものはもちろん理論を一から理解しようとする姿勢は重要ですが、ざっくりと理論を理解してしまったら一度RやPythonを用いて、実習してみるとより理解が深まると思います。
参考文献
https://www.tandfonline.com/doi/abs/10.1198/016214508000000337
https://www.rdocumentation.org/packages/monomvn/versions/1.9-13/topics/blasso