◆この記事のテーマ

ベイズファクターを使って比率の差を検定する方法を整理します。この手法は、例えば下図のようなシーンに適用できます。下図では異なる交差点の危険挙動発生割合を比較しています。

交通安全対策実施前後の危険挙動発生割合の変化を判定したいときにも同じ方法論が活用できます。

※一つ目の例は対応のない（two-sampe）グループ間の検定、二つ目の例は対応のある（one-sample）グループ間の検定ですが、どちらの検定も一つのモデルで包括して説明されます。

モデル

参考文献に掲載されたグラフィカルモデルを拝借して以下に示します。またここでは、一つ目の例の対応の無い交差点間の危険挙動割合の比較に着目して考えます。

ここで、$n_1$、$n_2$はそれぞれ交差点1、交差点2における車両毎の危険挙動に関するデータ（0=「危険挙動無し」 OR 1「危険挙動あり」）の標本数、$s_1$、$s_2$はそれぞれ交差点1、交差点2において危険挙動をした車両の台数です。交差点ごとに、独立の二項分布に従い危険挙動をした車両が観測されたと仮定し、発生確率を$\theta_1$、$\theta_2$とおきモデル化します。

比率の差の検定モデル $$ s_1 \sim \mathrm{Binomial}(\theta_1,n_1) \tag{1} $$ $$ \theta_1 \sim \mathrm{Uniform}(0,1) \tag{2} $$ $$ s_2 \sim \mathrm{Binomial}(\theta_2,n_2) \tag{3} $$ $$ \theta_2 \sim \mathrm{Uniform}(0,1) \tag{4} $$

このとき、我々の興味の対象となる確率変数は$\delta=\theta_1 - \theta_2$です。$\theta_1$、$\theta_2$の大小関係についての事前知識がある場合と無い場合に分け、2つの帰無仮説・対立仮説の組を設定します。交差点の比較の例では、交差点2は交差点1と同じ形状であり、信号機が設置されていることから、危険挙動発生率は交差点1より少ない、すなわち$\theta_1 > \theta_2$ではないか、と考えることができます。

$\theta_1$、$\theta_2$には、無情報事前分布$\mathrm{Uniform}(0,1)$を設定することとします。

◆帰無仮説と対立仮説 $\theta_1$と$\theta_2$の大小関係が明らかでないときの仮説は、 $$ \mathrm{Unrestricted~~Model:} \begin{cases} H_0:\delta=0 \\\\ H_1:\delta\neq0 \end{cases} $$ $\theta_1 \geq \theta_2$が明らかな場合の仮説は、 $$ \mathrm{Order-restricted~~Model:} \begin{cases} H_0:\delta=0 \\\\ H_2:\delta \geq 0 \end{cases} $$

ところで、ここまでグラフィカルモデルは対応のないグループ間の比率の差の検定を考えていましたが、冒頭で言ったように、対応のあるグループについても同じモデルで包括出来てしまいます。ここではその理由について説明しようと思います。

対応のあるグループ間の比率の差を考えたい場合のグラフィカルモデルは以下のようになります。

$\theta_1$、$\theta_2$、$\delta$の三角関係が変わっていますが、$\theta_1$と$\theta_2$まわりの尤度計算は対応の無い場合と同じです。さらに$\theta_1$、$\delta$に図示の無情報事前分布を設定してやると、$\theta_2$の事前分布は$\mathrm{Uniform}(0,1)$と、対応の無い場合と同じ事前分布になってしまいます。つまり、$\theta_1$、$\theta_2$の尤度関数と事前分布が対応の場合と同じになってしまうのです。$\delta$は$\theta_2-\theta_1$と、二つの確率変数の差の分布から求まりますから、$\delta$の事後分布は対応の無い場合の事後分布の$y$軸反転になります。よって、比率の差の検定においてはモデル上、対応の有無を問わず一つのモデルで説明出来てしまうのです。（以上、色々考えた上での私見なんですが、どうなんでしょう？）

ベイズファクターを解析的に求める

設定した仮説から分かる通り、本検定はネストされたモデル間の比較を行うので、以前の記事で紹介したSavage-Dickey法を使えばよいと分かります。Savage-Dickey法は、$H_1$に対応するモデルにおける事前分布と事後分布を一点比較をすればベイズファクターを算出できると主張しています。今回の事例では比較する点は$\delta=0$です。つまり

$$ BF_{01} = \cfrac{p(\delta=0 | D, H_1)}{p(\delta=0 |H_1)} \tag{5} $$

ということになります。ここで$D$は全データセットを示します。

この節では本モデルにおける$\delta$の事前分布、事後分布について考え、ベイズファクターを解析的に求めてみます。

まずは$\delta$の事前分布です。

◆$\delta$の事前分布 $\delta$の事前分布$p(\delta | H_1)$は重心0の三角分布となる。 $$ p(\delta | H_1) = \begin{cases} 1+\delta ~~~ \mathrm{for} ~ \delta \le 0 \\
1-\delta ~~~ \mathrm{for} ~ \delta > 0 \end{cases} \tag{6} $$

証明は以下を参照のこと。証明では確率変数の変換公式を応用して確率変数の和（差）の分布を求めるテクニックを用います。

◆$\delta$の事前分布の導出 確率密度$f_{\boldsymbol{X}^n}$を持つ$n$次元確率ベクトル$\boldsymbol{X}^n$に一対一の写像$\boldsymbol{\phi}^n$で$\boldsymbol{X}^n = \boldsymbol{\phi}^n(\boldsymbol{Y}^n)$と対応付けされる$n$次元確率ベクトル$\boldsymbol{Y}^n$の確率密度$g_{\boldsymbol{Y}^n}$を求める。 $$ \begin{cases} X_1 = \phi_1(Y_1,Y_2,\ldots,Y_n) \\\\ X_2 = \phi_2(Y_1,Y_2,\ldots,Y_n) \\\\ \vdots \\\\ X_n = \phi_n(Y_1,Y_2,\ldots,Y_n) \end{cases} $$ 領域$B = \prod_{i=1}^{n}(-\inf,y_i] \in \mathbb{R}^n$ での$g_{\boldsymbol{Y^n}}$の積分値と領域$\boldsymbol{\phi}(B)$での$f_{\boldsymbol{X^n}}$の積分値が一致することから、 $$ \int\cdots\int_B g_{\boldsymbol{Y^n}}(y_1,\ldots,y_n) dy_1\ldots dy_n = \left|\int\cdots\int_{\boldsymbol{\phi}(B)} f_{\boldsymbol{X^n}}(x_1,\ldots,x_n) dx_1 \ldots dx_n\right| $$ 右辺の多重積分について、積分変数を$(x_1,\ldots,x_n)$から$(y_1.\ldots,y_n)$に変換する。 $$ \left|\int\cdots\int_{\boldsymbol{\phi}(B)} f_{\boldsymbol{X^n}}(x_1,\ldots,x_n) dx_1\ldots dx_n\right| = \int\cdots\int_{B} f_{\boldsymbol{X^n}}(\phi_1(y_1,\ldots,y_n),\ldots,\phi_n(y_1,\ldots,y_n)) || J_{\boldsymbol{\phi}^n} || dy_1 \ldots dy_n $$ ここで$\left|\left|J_{\boldsymbol{\phi}^n}\right|\right|$はヤコビアンであり変数変換後の単位面積のスケールを調整する役割をもつ。 $$ \left|\left|J_{\boldsymbol{\phi}^n}\right|\right| = \cfrac{dy_1 dy_2\ldots dy_n}{dx_1 dx_2\ldots dx_n} $$ 以上より $$ \int\cdots\int_B g_{\boldsymbol{Y^n}}(y_1,\ldots,y_n) dy_1\ldots dy_n = \int\cdots\int_{B} f_{\boldsymbol{X^n}}(\phi_1(y_1,\ldots,y_n),\ldots,\phi_n(y_1,\ldots,y_n)) || J_{\boldsymbol{\phi}^n} || dy_1 \ldots dy_n $$ すなわち $$ g_{\boldsymbol{Y^n}}(y_1,\ldots,y_n) = f_{\boldsymbol{X^n}}(\phi_1(y_1,\ldots,y_n),\ldots,\phi_n(y_1,\ldots,y_n)) || J_{\boldsymbol{\phi}^n} || \tag{確率変数の変換公式} $$ 次に、確率変数$\Theta_1$、$\Theta_2$に対して、$\Delta = \Theta_1 - \Theta_2$の分布を求める。以下の一対一写像 $$ {\boldsymbol{\phi}^{2}}^{-1} = \begin{cases} \Delta = \Theta_1 - \Theta_2 \\\\ Ε = \Theta_2 \end{cases} $$ を考える。このとき、 $$ \boldsymbol{\phi}^{2} = \begin{cases} \Theta_1 = \Delta + Ε \\\\ \Theta_2 = Ε \end{cases} $$ ヤコビアン$\left|J_{\boldsymbol{\phi}^n}\right|$は、 $$ \left|J_{\boldsymbol{\phi}^n}\right| = \cfrac{d\delta d\epsilon}{d\theta_1 d\theta_2 } = \left| \begin{array}{ccc} \cfrac{\partial \theta_1}{\partial \epsilon} & \cfrac{\partial \theta_1}{\partial \delta} \\ \cfrac{\partial \theta_2}{\partial \epsilon} & \cfrac{\partial \theta_2}{\partial \delta} \end{array} \right| = \left|\begin{array}{ccc} 1 & 1 \\\\ 1 & 0 \end{array}\right| = -1 $$ $\Theta_1$、$\Theta_2$が独立の一様分布に従うと仮定したとき、確率変数の変数変換の公式より、$\Delta$、$Ε$の結合確率密度は、 $$ \begin{split} g_{\Delta,Ε}(\delta, ϵ) &=& f_{\Theta_1,\Theta_2}(\delta+\epsilon,\epsilon)\left|\left|J_{\boldsymbol{\phi}^2}\right|\right| = f_{\Theta_1,\Theta_2}(\delta+\epsilon,\epsilon) = f_{\Theta_1}(\delta+\epsilon)f_{\Theta_2}\epsilon\\\\ &=& \begin{cases} 1~~~\mathrm{for}~~~ 0 \leq \delta + \epsilon \leq 1~~\cap~~0 \leq \epsilon \leq 1 \\\\ 0~~~\mathrm{for~~~otherwise} \end{cases} \end{split} $$ これを$Ε$について積分し、$\Delta$の周辺確率密度を求める。 $-1\leq \delta < 0$のとき、$-\delta\leq \epsilon \leq 1$で$g_{\Delta,Ε}(\delta, ϵ)$は$1$をとるから、 $$ g_{\Delta}(\delta) = \int_{-\delta}^{1}1d\epsilon = 1 + \delta $$ $0 \leq \delta \leq 1$のとき、$0 \leq 1-\delta$で$g_{\Delta,Ε}(\delta, ϵ)$は$1$をとるから、 $$ g_{\Delta}(\delta) = \int_{0}^{1-\delta}1d\epsilon = 1 - \delta $$ 以上より、$\Delta$の確率密度は重心0の三角分布となる。 $$ g_{\Delta}(\delta) = \begin{cases} 1+\delta ~~~ \mathrm{for} ~ \delta \le 0 \\\\ 1-\delta ~~~ \mathrm{for} ~ \delta > 0 \end{cases} $$

事前分布について$p(\delta=0 |H_1)=1$であることがわかりました。よって、$BF_{01} = p(\delta=0 | D, H_1)$となります。これを求めた結果を示します。

◆比率の差の検定におけるベイズファクター 比率の差を検定するためのベイズファクター$BF_{01}$は $$ BF_{01} = \cfrac{{}_{n_{1}} C _ {s_{1}} {}_ {n_{2}} C _ {s_{2}}}{{}_ {n_1+n_2} C _ {s_1+s_2}}\cfrac{(n_1+1)(n_2+1)}{n_1+n_2+1} \tag{7} $$

証明は以下を参照のこと。周辺尤度を無視した事後分布の計算と、確率変数の変換公式を応用して確率変数の和（差）の分布を求めるテクニックを用います。

◆$p(\delta=0 | D, H_1)$の導出 $\delta = \theta_1 - \theta_2$であるから、$p(\theta_1 | D^{n_1},H_1)$と$p(\theta_2 | D^{n_2},H_1)$を求め、これらの差の分布を求めることを考える$(D = D^{n_1}+D^{n_2})$。一様分布$\mathrm{Uniform}(0,1)$はベータ分布$\mathrm{Beta}(1,1)$に等しいから、 $$ \begin{cases} p(\theta_1 | D^{n_1}, H_1) ∝ \mathrm{Binomial}(s_1|\theta_1, n_1) \mathrm{Beta}(\theta_1 | 1,1)　\\\\ p(\theta_2 | D^{n_2}, H_1) ∝ \mathrm{Binomial}(s_2|\theta_2, n_2) \mathrm{Beta}(\theta_2 | 1,1) \end{cases} $$ ここで、 $$ \begin{split} p(\theta_1 | D^{n_1}, H_1) &=& \cfrac{{}_{n_{1}} C _ {s_{1}}}{\mathrm{B}(1,1)}\theta_1^{s_1}(1-\theta_1)^{n_1 - s_1} \\\\ &∝& \mathrm{Beta}(\theta_1 | s_1+1, n_1-s_1+1) \end{split} $$ よって、 $$ \begin{cases} p(\theta_1 | D^{n_1}, H_1) = \mathrm{Beta}(\theta_1 | s_1+1, n_1-s_1+1)　\\\\ p(\theta_2 | D^{n_2}, H_1) = \mathrm{Beta}(\theta_2 | s_2+1, n_2-s_2+1) \end{cases} $$ ここで、確率変数$\Theta_1$、$\Theta_2$に対して、$\Delta = \Theta_1 - \Theta_2$の分布を求めると、 $$ g_{\Delta,Ε}(\delta, ϵ) = f_{\Theta_1,\Theta_2}(\delta+\epsilon,\epsilon)\left|\left|J_{\boldsymbol{\phi}^2}\right|\right| = f_{\Theta_1,\Theta_2}(\delta+\epsilon,\epsilon) = f_{\Theta_1}(\delta+\epsilon)f_{\Theta_2}\epsilon $$ これを$\epsilon$について積分して、 $$ g_{\Delta}(\delta) = \int_{-1}^{1}f_{\Theta_1}(\delta+\epsilon)f_{\Theta_2}\epsilon d\epsilon $$ $f_{\Theta_1}(\theta_1)$、$f_{\Theta_1}(\theta_1)$について $$ \begin{cases} f_{\Theta_1}(\theta_1) = \mathrm{Beta}(\theta_1 | s_1+1, n_1-s_1+1)　\\\\ f_{\Theta_1}(\theta_2) = \mathrm{Beta}(\theta_2 | s_2+1, n_2-s_2+1) \end{cases} $$ を代入すると、 $$ \begin{split} g_{\Delta}(\delta) &=& \int_{-1}^{1}f_{\Theta_1}(\delta+\epsilon)f_{\Theta_2}\epsilon d\epsilon \\\\ &=& \int_{-1}^{1} \mathrm{Beta}( \delta+\epsilon | s_1+1, n_1-s_1+1)\mathrm{Beta}(\delta | s_2+1, n_2-s_2+1) d\epsilon \\\\ &=& \int_{-1}^{1}\cfrac{(\delta+\epsilon)^{s_1}(1-\delta-\epsilon)^{n_1-s_1}\epsilon^{s_2}(1-\epsilon)^{n_2-s_2}}{\mathrm{B}(s_1+1,n_1-s_1+1)\mathrm{B}(s_2+1,n_2-s_2+1)}d\epsilon \end{split} $$ 求めたいのは $\left.g_{\Delta}(\delta)\right|_{\delta=0}$なので、 $$ \begin{split} \left.g_{\Delta}(\delta)\right|_{\delta=0} &=& \cfrac{1}{\mathrm{B}(s_1+1,n_1-s_1+1)\mathrm{B}(s_2+1,n_2-s_2+1)} \int_{-1}^{1}\epsilon^{s_1}(1-\epsilon)^{n_1-s_1}\epsilon^{s_2}(1-\epsilon)^{n_2-s_2}d\epsilon \\\\ &=& \cfrac{\mathrm{B}(s_1+s_2+1,n_1+n_2-s_1-s_2+1)}{\mathrm{B}(s_1+1,n_1-s_1+1)\mathrm{B}(s_2+1,n_2-s_2+1)} \int_{-1}^{1} \mathrm{Beta}(\epsilon | s_1+s_2+1,n_1+n_2-s_1-s_2+1)d \epsilon \\\\ &=&\cfrac{\mathrm{B}(s_1+s_2+1,n_1+n_2-s_1-s_2+1)}{\mathrm{B}(s_1+1,n_1-s_1+1)\mathrm{B}(s_2+1,n_2-s_2+1)} \\\\ &=& \cfrac{{}_{n_{1}} C _ {s_{1}} {}_ {n_{2}} C _ {s_{2}}}{{}_ {n_1+n_2} C _ {s_1+s_2}}\cfrac{(n_1+1)(n_2+1)}{n_1+n_2+1}~~~(\mathrm{B}(a,b)=\cfrac{(a-1)!(b-1)!}{(a+b-1)!}より) \end{split} $$ よって、 $$ \begin{split} p(\delta=0 | D, H_1) &=& \left( \left. \mathrm{Binomial}(s_1|\theta_1, n_1) \mathrm{Beta}(\theta_1 | 1,1) - \mathrm{Binomial}(s_2|\theta_2, n_2) \mathrm{Beta}(\theta_2 | 1,1) \right) \right | _{\theta_1 - \theta_2=0} \\\\ &=& \cfrac{{}_ {n_{1}} C _ {s_{1}} {}_ {n_{2}} C _ {s_{2}}}{{}_ {n_1+n_2} C _ {s_1+s_2}}\cfrac{(n_1+1)(n_2+1)}{n_1+n_2+1} \end{split} $$

おわりに

今回は、比率の差をベイズファクターを使って検定する手法を整理しました。本内容の実践版はこちらの姉妹記事を参照してください。姉妹記事では解析的にベイズファクターを計算するだけではなく、MCMCを用いた近似手法でも計算しています。なぜ精度の劣る方法で…と思うかもしれませんが、従来通りの検定に対応する帰無仮説・対立仮設の設定に留まらず、片側検定に対応する仮説の比較をすることもできることにそのメリットがあります。詳しくは姉妹記事を参照のこと。

2021.04.18 00:00

ベイズファクターを用いた仮説検定～比率の差の検定～

モデル

ベイズファクターを解析的に求める

おわりに