Rを使ってBayesian Lassoで遊んでみた

今回はせっかくBayesian Lassoを使って遊んでみたのでその記録を残しておきたいと思います。(Bayesian Lassoを扱った記事も少ないように思われたので)

# Bayesian Lassoとは

そもそもBayesian Lassoとは何なのでしょうか。
厳密なことであったり、詳しいことは[Park & Casella (2008)](https://www.tandfonline.com/doi/abs/10.1198/016214508000000337)にあるので気になる方はそちらもご参照ください。ここでは簡単に説明するにとどめます。通常の最小二乗法の回帰や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}
$$
の正規分布の尺度混合分布に階層化し、事後分布を推定するために以下のモデルが立式できます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/634709/3c81be55-21b1-28e4-db8a-41017f9ba7a8.png)



# 実際にBayesian Lasso で遊んでみる。

理論的なことをざっくりと把握したので(本当にざっくり)、次は実際に遊んでみましょう。習うより慣れろです。
今回は題材としてはRでBayesian Lassoを行うときに用いるblassoの[公式ドキュメント](https://www.rdocumentation.org/packages/monomvn/versions/1.9-13/topics/blasso)を用います。(公式ドキュメント大事！)
ここにあるサンプルコードを1行ずつ追ってみたいと思います。

```r
library(monomvn); library(lars); library(glmnet);
```
まず、今回必要なライブラリを読み込みたいと思います。インストールしていない場合はインストールしてから実行してください。(install.packageとして直接コマンドに打ち込んでもいいですが、Rstudio)を用いている方だと右下あたりにInstallというところがるので、そちらからインストールしてもよいでしょう。)


```r

data(diabetes)
attach(diabetes)
```
ここでは今回用いるデータを読み込みます。有名なdiabetesのデータセットを用います。
attach()関数を使うと、そのあとは、どのデータフレームを用いるかをかを指定しなくてよくなるため便利です。なので、ここでattachもしておきましょう。

```r
## 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によってでてきたもの点がそのほかの二つの回帰によって得られた推定値となります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/634709/0019a918-abd6-1506-966b-b1241ddd5a22.png)

```r
plot(reg.blas, burnin=200, which="m")
```
ここではblassoによって得られたmの分布を表示します。ここでいうmとはBayesian Lasso回帰を行った際に係数が0でないものの数です。これを見ると7個の係数が0でないと判断された時が一番多いことが分かりますね。下のmchainとは1000回の試行でmがどのように推移していったかを示すものです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/634709/ba329e42-ff88-ad3c-9687-70975326f9ac.png)

```r
## 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だけでなくすべての結果を出力すると以下のようになります。

```r

> 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ならそれは、通常の最小二乗法による回帰が行われていることを表します。

``` r
## 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と同じです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/634709/517a0a76-6382-e29d-8f98-4d315d7140d5.png)

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/634709/4fb1f3ba-7c52-4559-9646-a1f6ac0ca8c1.png)

以上簡単ですが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