２×２クロス表のベイズ分析
２×２のクロス表においてセル内のデータ数が少ないときはフィッシャーの直接確率（Fisher‘s Exact Test）が用いられるが、フィッシャーの直接確率では、行和と列和が固定されている（Conover, 1999）。この条件が不適当なときには、この制約条件を課さない簡単な確率モデルを設定してベイズ分析を用いることができる。また、Fissha-no 直接確率では、分析結果はｐ値のみであるが、ベイズ分析では、確率の差がどれほどであるかという情報が事後分布として得られる。
行を表す変数をX（＝１、２）とおき、列を表す変数をY（＝１、２）とおく。データの４つのセルへの分類は独立であるとする。
行への分類確率を
とおけば、
である。
列への分類確率を、行への分類を条件として、以下のように設定する。
このとき
である。
したがって、４つのセルへの分類確率は以下のようになる。
|
|
|
|
|
|
|
|
|
|
|
|
この確率に基づいて多項分布を適用するStanスクリプトをリスト１のように用意した。
リスト１ 2×２クロス表に多項分布を適用するStanスクリプト（xtable2by2.stan）
data {
int pop1[2];
int pop2[2];
}
transformed data {
int props[4];
vector[4] a;
for (i in 1:2) {
props[i] = pop1[i];
props[2+i] = pop2[i];
}
for (i in 1:4) {
a[i]
= 1.0;
}
}
parameters {
real<lower = 0.0, upper =
1.0> p0;
real<lower = 0.0, upper =
1.0> p1;
real<lower = 0.0, upper =
1.0> p2;
}
transformed parameters {
simplex[4] theta;
theta[1] = p0 * p1;
theta[2] = p0 * (1 - p1);
theta[3] = (1 - p0) * p2;
theta[4] = (1 - p0) * (1 -
p2);
}
model {
p0 ~ uniform(0, 1);
p1 ~ uniform(0, 1);
p2 ~ uniform(0, 1);
props ~ multinomial(theta);
}
generated quantities {
real p1_p2;
p1_p2 = p1 - p2;
}
リスト１では、pop1に（X＝１，Y＝１）の度数と（X＝１、Y=2）の度数を設定し、pop2に（X＝２，Y=1）の度数と（X＝２，Y＝２）の度数を設定する。
リスト１のStanスクリプトを用いるPythonスクリプトをリスト２のように用意した。
リスト２ ２×２クロス表のベイズ分析Pythonスクリプト（XTable2by2Bsn.py）
import numpy as np
import matplotlib.pyplot as plt
import pystan
import seaborn as sb
import pickle
s = input('Pop.1 = ').split()
pop1 = [int(v) for v in s]
s = input('Pop.2 = ').split()
pop2 = [int(v) for v in s]
print('Pop1 =', pop1)
print('Pop2 =', pop2)
sm = pystan.StanModel(file =
'xtable2by2.stan')
fit = sm.sampling(data = {'pop1':pop1,
'pop2':pop2}, n_jobs = 1)
print(fit)
sb.kdeplot(fit['p0'])
plt.title('Posterior Distribution of
p0', fontsize = 16)
plt.savefig('Figp0.png')
plt.show()
sb.kdeplot(fit['p1'], label = 'p1')
sb.kdeplot(fit['p2'], label = 'p2')
plt.title('Posterior Distributions of
p1 and p2', fontsize = 16)
plt.legend(fontsize = 16)
plt.savefig('Figp1p2.png')
plt.show()
ckL = [1 if v > 0 else 0 for v in
fit['p1_p2']]
p_pos = sum(ckL) / len(ckL)
Q1, Med, Q3 =
np.percentile(fit['p1_p2'], [25, 50, 75])
plt.title('Prop(p1-p2>0) =
{0:.2f}'.format(p_pos) \
+ '\nQ1 = {0:.2f}, Med =
{1:.2f}, Q3 = {2:.2f}'.
format(Q1, Med, Q3), fontsize = 16)
sb.kdeplot(fit['p1_p2'], label =
'p1-p2')
plt.xlabel('p1 - p2', fontsize = 14)
plt.legend(fontsize = 14)
plt.tight_layout()
plt.savefig('FigDiffp1p2.png')
plt.show()
リスト１とリスト２のスクリプトファイルは圧縮ファイルScriptFiles.zipにまとめたので、自由に利用されたい。
リスト２のPythonスクリプトは、AnacondaのPython 3.7仮想環境にPyStan 2.17.1.0がインストールされている環境で実行されたものである。
いま、Conover (1999), p. 190の例、
Malesで、account representativeが１人、tellerが9人
Femalesで、account representativeが３人、tellerが１人
の分析を行ってみる。
リスト２のスクリプトを実行すると、クロス表の値の入力が求められる。
次の入力では、Pop.1 =の右辺に行（X＝１）の値として、Malesの１と９が空白文字で区切って設定され、
Pop.2 =の右辺に行（X＝２）の値として、Femalesの3と１が空白文字で区切って設定されている。
(py37) PS
D:\XXX…XX\ScriptFiles> python .\XTable2by2Bsn.py
INFO:numexpr.utils:NumExpr
defaulting to 8 threads.
Pop.1 = 1 9
Pop.2 = 3 1
Pop.2 =の右辺の設定後、StanスクリプトのコンパイルとMCMCサンプリングが始まる。
MCMCサンプリング後、まずp0の事後分布が表示される（図１）。
図１
図１のフォームを閉じると、図２が表示される。
図２
パラメータp1とp2の事後分布が表示されている。
図2のフォームを閉じると図３のグラフが表示される。
図３
確率のパラメータp1とp2の差
の分布は、正の領域は0.01（1％）しかなく、99％が負の領域にある。分布のメディアンも―0.52であり、p1の方がp2より実質的に小さいと判断できる。すなわちFemalesの方がaccount representativeである確率がMalesよりも大きいことが示されている。
Conover (1999)に示されているフィッシャーのthe exact lower-tailed p-valueは、0.041である（p. 190）。p値を示されるより、事後分布（図2、図３）の方が、確率p1とp2の違いに関する情報がより豊かであると言えよう。
なお、Python用Stan（PyStan）の説明を＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞で行っている。
PyStanのインストールについては、このウェブサイトなどで説明している。
参考文献
Conover, W. J. (1999) Practical nonparametric statistics. John Wiley & Sons, Inc.