ベイズ分位点回帰

こんにちは，株式会社Nospareリサーチャー・千葉大学の小林です．
前回の記事では，分位点回帰モデルの紹介と頻度論の立場からの適用について書きました．ベイズ分位点回帰モデルはその拡張性の高さから大きな関心を集めてきました．今回はベイズ統計学の枠組みで分位点回帰を取り上げたいと思います．

ベイズ統計学においては尤度関数を定義する必要があります．ベイズ分位点回帰モデルでは尤度の作り方がいくつかあるのですが，ここでは最も簡単なアプローチとして非対称ラプラス分布を用いたものを取り上げます．非対称ラプラス分布の密度関数は次のように与えられます．

f(x;\mu,\tau,\sigma)=\frac{\tau(1-\tau)}{\sigma}\exp\left\{-\rho_\tau\left(\frac{x-\mu}{\sigma}\right)\right\},\quad -\infty<x<\infty,

ここで$\mu$は位置パラメータ，$\sigma>0$は尺度パラメータ，$\alpha\in(0,1)$は形状パラメータで，$\mu$がちょうど分布の$\tau$分位点に位置しています．
指数の中の$\rho_\tau$は通常の分位点回帰の目的関数であるチェックロスです

\rho_\tau(u)=u(\tau-I(u<0)).

非対称ラプラス分布$\mu=0$，$\sigma=$の密度分布は下図のような形になります．

ベイズ分位点回帰モデルでは回帰の誤差項がこの非対称ラプラス分布に従うと仮定します．この密度関数はそのまま扱うこともできなくないですが，以下のような正規分布と指数分布の混合表現を使うとより効率的にMCMCを適用することができます：

\epsilon = v\theta_1 + \sqrt{\sigma v \theta_2^2}z

ここで$v\sim Exp(\sigma)$，$z\sim N(0,1)$，$\theta_1=(1-2\tau)/\tau(1-\tau)$，$\theta_2^2=2/\tau(1-\tau)$．指数分布に従う$v$を条件つけることで，$\epsilon|v$は正規分布に従うことがわかります．データ拡大法によって$v$を完全条件付分布からサンプリングするだけで，他のパラメータのサンプリングは正規分布に基づいた回帰モデルの枠組みでギブスサンプラーのが構築できます．この定式化は打ち切りのある分位点回帰を扱うのにも非常に便利になります．

bayesQRパッケージを使ってみる

bayesQRパッケージでは非対称ラプラス分布に基づいた線形分位点回帰モデルをMCMCによってベイズ推定します．bayesQRパッケージで中心となるのはbayesQR()関数です．この関数では連続と２値（0か1）の被説明変数に対する分位点回帰モデルを扱います（個人的には２値変数に対する分位点回帰はあまりオススメしません）．係数パラメータ$\beta_\tau$について通常の事前分布およびadaptive Lassoに基づいた事前分布を指定できます．

• formula：分位点回帰モデルの記述
• quantile：分位点のレベル（デフォルトで0.5），ベクトル指定可能
• alasso：adaptive Lassoを使うかを表す論理値（デフォルトではFALSE）
• prior：$\beta_\tau$と$\sigma$の事前分布についての設定，prior()関数で設定したものを使う
• ndraw：MCMCの繰り返し回数
• keep：thinningの数，ndraw個のMCMCサンプルのうちkeep番目のサンプルだけをとっておく
• normal.approx：事後分布を正規近似するかどうかの論理値（デフォルトではTRUE）．TRUEにしておくと$\sigma$の事前分布の設定が不要になる

prior()関数においてパラメータの事前分布を設定します．この人工データの例ではadaptive LASSOを使わずに$\beta_\tau\sim N(0,1000I)$，$\sigma\sim IG(3,0.1)$（逆ガンマ分布）とするので，beta0，V0，shape0，scale0の値をそれぞれ行列またはスカラーで指定します（これらの事前分布のデフォルトは$N(0,100I)$，$IG(0.01,0.01)$です)．

library(bayesQR)
# データを発生させる
set.seed(1)
n <- 200
x <- rnorm(n)
y <- rnorm(n, mean=1+2.0*x)
# 事前分布を設定する．
prior <- prior(y~x, beta0=matrix(0,nrow=2,ncol=1), V0=1000*diag(2), shape0=3.0, scale0=0.1)

事前分布の指定ができたらモデルの推定を$\tau=0.01,0.02,\dots,0.98,0.99$で行います（quantregパッケージよりはちょっと時間はかかります）．

tau <- seq(0.05,0.95,0.05)
res <- bayesQR(y~x, prior=prior, quantile=tau, ndraw=10000, keep=1, normal.approx=FALSE)

bayesQRではquantregと同様に回帰係数の事後分布を$\tau$の関数としてプロットするメソッドがついています．plot()関数でplottype="quantile"と指定すると，var番目の係数パラメータの各分位点における事後平均と95%信用区間をつなげて描画します（変数ごとに描画させるためループを使っています）．一方でplottype="trace"や"hist"と指定すると，各quantile分位点におけるパラメータのMCMCのトレースプロットやヒストグラムを描画します．

par(mfrow=c(1,2))
for (i in 1:2){
    plot(res, plottype="quantile", var=i)
}

最後に$\tau=0.1, 0.5, 0.9$について条件付分位点の事後平均をプロットしてみましょう．

ind <- c(2,10,18)
plot(y=y, x=x)
for (i in 1:3){
    # 事後平均を計算
    beta <- apply(res[[ind[i]]]$betadraw, 2, mean)
    abline(a=beta[1], b=beta[2])    
}

Brqパッケージを使ってみる

Brqパッケージでも分位点回帰のベイズ推定を行うことができます．このパッケージのいいところは打ち切りのある被説明変数を扱えたり，LASSOや適応的LASSOによる変数選択も行えたりすることろにあります．一方で事前分布のハイパーパラメータの指定がほぼできないという欠点もあります．これによりあと共変量の関数をノンパラメトリックに推定することも事前分布による罰則を導入できないということでちゃんと推定することができないです．（あとマニュアルが不親切です）．

基本となるBrq()関数はbayesQR()と使い方はほぼ同じですが引数methodに指定する内容によっていろいろなモデルを推定することができます．例えば

• Bqr：普通のベイズ分位点回帰
• BALqr：適応的LASSOによるベイズ分位点回帰
• BLqr：LASSOによるベイズ分位点回帰
• Btqr：左側打ち切りのあるベイズ分位点回帰モデル

などです．

先程の生成したデータをそのまま使ってBrq()関数を試してみましょう．

library(Brq)
tau <- c(0.1,0.5,0.9)
plot(y=y, x=x)
res <- Brq(y~x, tau=tau, runs=11000, burn=1000)
for (i in 1:3)
    abline(a=res$coefficients[1,i], b=res$coefficients[2,i])
}
`summary(res)`

直線は事後平均を表していますが，推定結果はbayesQRパッケージと同じですね．事後分布の要約を見ると，各分位点において事後平均は真値と近く，95%信用区間も真値を含んでいることがわかります．

Call:
Brq.formula(formula = y ~ x, tau = 0.5, runs = 11000, burn = 1000)

tau:[1] 0.5

          Estimate L.CredIntv U.CredIntv
Intercept 1.039527  0.9006036   1.194918
x         1.951499  1.7992500   2.124598
> res <- Brq(y~x, tau=tau, runs=11000, burn=1000)
> summary(res)
$tau
[1] 0.1

$coefficients
       Estimate L.CredIntv U.CredIntv
[1,] -0.3190949 -0.4637631 -0.1682399
[2,]  2.1343379  1.9954108  2.2722162

$tau
[1] 0.5

$coefficients
     Estimate L.CredIntv U.CredIntv
[1,] 1.045066  0.9028536   1.201596
[2,] 1.955265  1.8022694   2.126541

$tau
[1] 0.9

$coefficients
     Estimate L.CredIntv U.CredIntv
[1,] 2.314270   2.167803   2.474030
[2,] 1.950449   1.822644   2.103496

$\tau=0.5$の場合をもう一度推定して，各係数パラメータの事後分布からのサンプルのヒストグラムをプロットしてみましょう．

res <- Brq(y~x, tau=0.5, runs=11000, burn=1000)
par(mfrow=c(1,2))
# 引数Coefficientsにループでインデックスを渡す
for (j in 1:2){
    plot(res, plottype="hist", Coefficients=j)
}

結果は以下のとおりです．ラベルが変ですが左が定数項，右が$x$の係数パラメータの事後分布です．真値の周辺に事後分布が集中しているのがわかります．

最後に

ベイズ分位点回帰モデルのRパッケージを使った推定について取り上げました．本記事で取り上げた方法についてもっと詳しく知りたい方はKozumi and Kobayashi (2011)を参考にしてください．ベイズ分位点回帰モデルについてはまた今後もまた記事を書いていきたいと思います．

参考文献

1. Kozumi, H. and Kobayashi, G. (2011). Gibbs sampling methods for Bayesian quantile regression. Journal of Statistical Computation and Simulation, 81, 1565-1578.

株式会社Nospareでは分位点回帰モデルに限らず，統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては弊社 までお問い合わせください．