こんにちは，株式会社Nospareリサーチャー・千葉大学の小林です．２項ロジットモデルは，ある試験に合格・不合格になる，消費者がある製品カテゴリから購入する・しない，ある疾患にかかる・かからない，といった結果に対してそれらが起こる確率が説明変数の値によってどのように変化するかを調べることができます．また２項ロジットモデルの一般化である多項ロジットモデルはある製品カテゴリ内にある複数ブランドからあるブランドをひとつ選択するといったような行動が観測されたデータを分析するのに用いられます．さらにロジットの構造は被説明変数が単純に何らかの選択行動の結果である場合に対する回帰モデルだけでなく，mixture of experts model（混合エキスパートモデル）やzero inflation（カウントデータにおけるゼロ過剰）などといったようなデータに内在するカテゴリーの分類確率を記述するのにも用いられ，とても幅広い応用の場面があります．本記事は前半部分にあたり，ロジットモデルをベイズ推定するための便利なモデル定式化とその応用例を紹介します．
２項ロジットモデル
まず２項ロジットモデルを紹介します．被説明変数を$y_i$, $i=1,\dots,N$とします．$y_i$は$0$か$1$の値を取り，例えば来店時にある製品カテゴリからの購入があった場合に$y_i=1$，なかった場合に$y_i=0$が観測されます．ロジットモデルでは$y_i=1$が起こる確率（ここでは選択確率と呼ぶことにします）を説明変数のベクトル$x_i$を用いて次のようにモデル化します．
\Pr(y_i=1)=\frac{\exp(x_i'\beta)}{1+\exp(x_i'\beta)},\quad i=1,\dots,N.
ここで$\beta$は$x_i$に対する係数パラメータです．簡単に言えば$y_i$成功確率が上のようなベルヌーイ分布に従います．この確率をもとに，尤度関数に対する$i$番目のデータについての部分は
p(y_i|\beta)=\left[\frac{\exp(x_i'\beta)}{1+\exp(x_i'\beta)}\right]^{y_i} \left[\frac{1}{1+\exp(x_i'\beta)}\right]^{1-y_i},\quad i=1,\dots,N,
と書けます．尤度関数は$p(y|\beta)=\prod_{i=1}^Np(y_i|\beta)$で，最尤法の場合にはこのまま尤度関数を$\beta$に関して最大化しますが，解析的に最大化は行えないので数値的な手法を利用します．ベイズでは$\beta$の事前分布$p(\beta)$を特定し，事後分布$p(\beta|y)\propto p(y|\beta)p(\beta)$に対する推測を行います．尤度関数が$\beta$についてきれいな形ではないので，きれいなギブスサンプラーは構築できず，メトロポリス・ヘイスティングス（MH）アルゴリズムなどによって事後分布からのサンプリングを行うのが普通でした．
ベイズではプロビットのほうが人気だった（読み飛ばしても大丈夫）
二項選択問題に対して最尤法による分析をする場合には二項ロジットモデルはよく用いられてきました．一方でベイズ推定においては二項選択問題においてもうひとつの定番的な手法であるプロビットモデルが定番で1990年代以降用いられてきました．これにはプロビットモデルに対するMCMCはデータ拡大法を使うことでギブスサンプリングだけで済ますことができるということが大きな理由だと考えられます．プロビットモデルは選択確率を$\Pr(y_i=1)=\Phi(x_i'\beta)$（$\Phi$は標準正規分布の分布関数）と書けます．ここで潜在変数$y_i^{*}$を導入し，$y_i^{*}=x_i'\beta+e_i$という回帰モデルを考えます（$e_i\sim N(0,1)$）．この潜在変数$y_i^{*}$と実際のデータ$y_i$を，$y_i^{*}\geq 0$のときに$y_i=1$が観測される，$y_i^*<0$のときに$y_i=0$が観測される，というように結びつけます（$y_i^{*}$を積分消去するともとのプロビットの選択確率が得られます）．プロビットモデルに対するギブスランプラーは$y_i^{*}$を完全条件付分布（切断正規分布）からサンプリングすることで，あたかも普通の回帰モデルで被説明変数が観測されているかのように$\beta$の完全条件付分布からサンプリングを行うことができます．
二項ロジットモデルに対しても同じようなデータ拡大法がHeld and Holmes (2006) によって提案されました．再び$y_i^{*}=x_i'\beta+e_i$という潜在変数に対する回帰モデルを考えます．プロビットモデルと異なり，ロジットモデルでは$e_i$はロジスティック分布に従います．この論文はロジスティック分布に対して次のような表現のもとにギブスサンプラーを構築しました：$e_i\sim N(0,\lambda_i)$, $\lambda=(2\psi_i)^2$, そして$\psi_i$はKolmogorov-Smirnov（KS）分布という分布に従います．このアプローチでのギブスサンプリングは，$y_i^{*}$，$\lambda_i$，$\beta$のサンプリングステップからなりますが，$\lambda_i$の完全条件付分布が標準的なものでなく，サンプリングがやや面倒というのもあって今まであまり使われてこなかったのかと思います．また予めパラメータとウェイトの値を定めておいた正規分布の有限混合モデルによってロジスティック分布を近似するという方法も提案されましたが，あくまでも近似モデルからのサンプルになってしまうということで定番的な方法として定着することはなかったと思います．
これは便利かも！Polya-Gamma mixtureによる定式化
さて遠回りになりましたが，本題のロジットモデルに対する画期的な定式化がPolson etal.(2013)によって提案されました．この論文の結果（Theorem 1）はロジスティック型の関数を次の積分で表現できるというものです．
\frac{(e^\psi)^a}{(1+e^\psi)^b}=2^{-b}e^{\kappa\psi}\int_0^\infty e^{-\frac{w \psi^2}{2}}p(w|b,0)dw,\quad \kappa=a-b/2, \quad (1)
ここで$w$はPolya-Gamma分布$PG(b,0)$という分布に従う変数で$p(w|b,0)$という密度関数をもつとします（これについては後で書きます）．
$(1)$でまず重要なのは$w$を所与にしてみると，積分の中に$\psi$について正規分布の密度関数の核が現れるということです．$\psi$は$x'\beta$のように回帰式で表される部分なので，この表現を利用するとデータと$w$が与えられたもとで$\beta$の完全条件付分布が正規分布になりそうな予感がします．実際に二項ロジットモデルの場合には，
p(y_i|\beta)=\frac{\exp(x_i'\beta)^{y_i}}{1+\exp(x_i'\beta)}
と書けるので，上の結果$(1)$を使って
p(y_i|\beta)\propto \exp(\kappa_i x_i'\beta)\int_0^\infty \exp\left\{-\frac{w_i}{2}(x_i'\beta)^2\right\}p(w_i|b,0)dw_i
となります（$a=y_i$, $b=1$, $\kappa_i=y_i-1/2$）．
$\beta$の事前分布を$p(\beta)$とすると，データと$w=(w_1,\dots,w_N)$を所与としたときの$\beta$の完全条件付分布は
\begin{align}
p(\beta|w,y)&\propto p(\beta)\prod_{i=1}^N p(y_i|\beta) \\
&=p(\beta)\prod_{i=1}^N \exp\left\{-\frac{w_i}{2}((x_i'\beta)^2-2\kappa_i x_i'\beta/w_i)\right\}\\
&\propto p(\beta)\prod_{i=1}^N\exp\left\{-\frac{w_i}{2}(\kappa_i/w_i-x_i'\beta)^2\right\}\\
&\propto p(\beta)\exp\left\{-\frac{1}{2}(z-X\beta)'W(z-X\beta)\right\}
\end{align}
という形になります（$z=(z_i,\dots,z_N)'$, $z_i=\kappa_i/w_i$, $X=(x_1',\dots,x_N')$, $W=diag(w_1,\dots,w_N)$）．このことから例えば$\beta$に対して正規事前分布を仮定すると完全条件付分布も正規分布になることがわかります．
さて$(1)$において$w$が従うPolya-Gamma分布ですが，一般的に$X\sim PG(b,c)$の場合には
X= \frac{1}{2\pi}\sum_{k=1}^\infty \frac{g_k}{(k-1/2)^2+c^2/(4\pi^2)}, \quad g_k\sim Ga(b,1),
という独立なガンマ変数の無限和で表現されます．$PG(b,c)$の密度関数は$PG(b,0)$の密度関数のexponential tilting（指数傾斜）で表されます：
p(x|b,c)=\frac{\exp\left\{-\frac{c^2x}{2}\right\}}{E\left[\exp \left(-\frac{c^2}{2}w\right)\right]}p(x|b,0)
ここで$p(x|b,0)$は$PG(b,0)$の密度関数で，分母の期待値は$w\sim PG(b,0)$についてとったものです．では，$(1)$において$w|\psi$の条件付分布を考えてみると，
p(w|\psi)=\frac{ \exp\left\{-\frac{w \psi^2}{2}\right\}p(w|b,0)}{\int\exp\left\{-\frac{w \psi^2}{2}\right\}p(w|b,0)dw}
となることから，$w|\psi\sim PG(b,\psi)$となることがわかります．
二項ロジットモデルの場合には，$w_i$の完全条件付分布は単純に$PG(1,x_i'\beta)$となります．よって$\beta\sim N(b_0,B_0)$の場合，次のふたつのステップから成るギブスサンプラーを回すことになります．
- $i=1,\dots,N$まで，$w_i$を$PG(1,x_i'\beta)$からサンプリングする
- $\beta$を$N(b_1,B_1)$からサンプリングする，ただし$B_1=\left[X'W X+ B_0^{-1}\right]^{-1}$, $b_1=B_1\left[X'W z+ B_0^{-1}b_0\right]$
この論文では$PG(1,c)$からの効率的なサンプリング方法を提案しています．指数分布と逆ガウシアン分布に基づく棄却サンプリングによるものですが採択確率が99.9%以上という極めて効率的なものとなっています．$n$が自然数のときには$n$個の独立な$PG(1,c)$変数を足し合わせるだけで$PG(n,c)$からのサンプリングができます．一般的な$PG(r,c)$の場合にはWindle et al.(2014)で提案されています．Rでの実装においてはBayesLogitパッケージ 内のrpg()関数が利用できます．このパッケージを利用した実装は後編で紹介したいと思います．
Polya-Gammaはいろいろなモデルで使える！
多項ロジットモデル
この定式化を利用することで，二項ロジットモデルの一般化である多項ロジットモデルに対するギブスサンプラーも構築することができます．多項ロジットモデルでは$J$個のカテゴリ（ブランドなど）からひとつのカテゴリが選択・実現されるという状態を考えます．被説明変数を$y_{ij}$とし，$j$番目のカテゴリが選択された場合に$y_{ij}=1$，$y_{ih}=0$, $h\neq j$が観測されます．多項ロジットモデルの選択確率$p_{ij}=\Pr(y_{ij}=1)$は
p_{ij}=\frac{\exp(\psi_{ij})}{\sum_{h=1}^J\exp(\psi_{ih})}
で与えられ，$\psi_{ij}=x_i'\beta_j$とモデル化します．尤度関数は
f(y|\beta)=\prod_{i=1}^N\prod_{j=1}^J\left(\frac{\exp(\psi_{ij})}{\sum_{h=1}^J\exp(\psi_{ih})}\right)^{y_{ij}}
となります．一見Polya-Gammaが使えそうにないですが，$\beta_j$以外のパラメータ$\beta_{-j}$を条件付けると
f(y|\beta_j,\beta_{-j})=\prod_{i=1}^N\left(\frac{e^{\eta_{ij}}}{1+e^{\eta_{ij}}}\right)^{y_{ij}}\left(\frac{1}{1+e^{\eta_{ij}}}\right)^{1-y_{ij}}
と書けます（$\eta_{ij}=x_i'\beta_j-C_{ij}$, $C_{ij}=\log \sum_{h\neq j}e^{x_i'\beta_h}$）．こうすると冒頭の二項ロジットモデルの尤度のような形になっていて，$\kappa_{ij}=y_{ij}-1/2$, $b=1$でPolya-Gammaが使えることがわかります：
\prod_{i=1}^N\exp\left\{\kappa_{ij}\eta_{ij}-\frac{w_{ij}}{2}\eta_{ij}^2\right\}p(w_{ij}|b,0).
$\beta_j$の事前分布を$N(b_{j0}, B_{j0})$とすると次のギブスサンプラーを構築できます．
- $i=1,\dots,N$, $j=1,\dots,J$について$w_{ij}$を$PG(1,\eta_{ij})$からサンプリングする．
- $j=1,\dots,J$について$\beta_j$を$N(b_{1j},B_{1j})$からサンプリングする，ただし$B_{1j}=\left[X'W_j X+ B_{j0}^{-1}\right]^{-1}$, $b_{1j}=B_{1j}\left[X'W_j (C_j+z_j) + B_0^{-1}b_0\right]$, $C_j=(C_{1j},\dots,C_{Nj})'$, $z_j=(z_{1j},\dots,z_{Nj})'$, $z_{ij}=\kappa_{ij}/w_{ij}$.
負の二項分布（Negative binomial）
負の二項分布はポアソン分布に代わるカウントデータに対する確率分布として広く使われてきました．確率関数は
p(y_i)=\frac{\Gamma(\xi+y_i)}{\Gamma(\xi) y_i! } (1-\psi_i)^\xi\psi_i^{y_i},\quad y_i=0,1,2,\ldots
と書け（$\xi>0$はoverdispersionのパラメータ），負の二項分布における回帰モデルは$\psi_i=e^{x_i'\beta}/(1+e^{x_i'\beta})$と与えられます．ここでまたPolya-Gammaが使える形になっていることに気づきます．実際，
(1-\psi_i)^\xi \psi_i^{y_i} =\frac{\exp(x_i'\beta)^{y_i}}{[1+\exp(x_i'\beta)]^{\xi+y_i}}
なので， $a_i=y_i$, $b_i=\xi+y_i$とすればあとは二項ロジットの場合と同じようなギブスサンプラーが構築できます．事前分布$\beta\sim N(b_0,B_0)$の場合，
1. $i=1,\dots,N$まで，$w_i$を$PG(y_i+\xi,x_i'\beta)$からサンプリングする．
2. $\beta$を$N(b_1,B_1)$からサンプリングする，ただし$B_1=\left[X'W X+ B_0^{-1}\right]^{-1}$, $b_1=B_1\left[X'W z+ B_0^{-1}b_0\right]$, $z=(z_1,\dots,z_N)'$, $z_i=\frac{y_i-\xi}{2w_i}$.
3. $\xi$をMHアルゴリズムなどでサンプリングする．
となります．
おわりに
前篇となる本記事ではPolya-Gammaを使った２項ロジットモデルの再定式化から始まり，Polya-Gammaを使うと幅広いクラスのモデルに対してギブスサンプラーを構築できるようになることを紹介しました．ここでは$\psi$や$\beta$について単純な回帰モデルや事前分布を考えていましたが，時間相関や空間相関が含まれるより複雑な階層モデルや，説明変数の選択のための馬蹄事前分布（horseshoe prior）などといったより高度な$\beta$の階層事前分布のもとでの推定も簡単な拡張で実行できます．後編ではRによる実装を紹介したいと思います．
一緒にお仕事しませんか！
株式会社Nospareではベイズ統計学に限らず統計学の様々な分野を専門とする研究者が所属しており，新たな知見を日々追求しています．統計アドバイザリーやビジネスデータの分析につきましては弊社までお問い合わせください．
