２×２クロス表のベイズ分析
2025.11改訂
Spiegelhalter(2019)は、腕組み（arm crossing）において左腕が上になる確率に性差があるかどうかの検定例を取り上げている。
女性14名において5名が左腕が上であり、男性40名において17名が左腕が上であった。Spiegelhalter(2019)では、カイ２乗検定が説明されているが、ベイズ分析を試みる。
いま、2つの条件、列1と列２、のそれぞれにおいて2つの属性、行１と行２、のデータが表1に示す個数得られたとする。Spiegelhalterの例では、列１が女性、列２が男性であり、行１が左腕が上、行２が右腕が上である。
|
表１ ２×２分割表 |
||
|
|
列１ |
列２ |
|
行１ |
|
|
|
行２ |
|
|
すなわち、表２のようになる。
|
表２ ２×２分割表（表１の改変） |
||
|
|
女性 |
男性 |
|
左腕が上 |
|
|
|
右腕が上 |
|
|
女性の条件（列１）において左腕が上（行１）である確率を
、男性（列２）において左腕が上（行１）である確率を
とおき、次の２項分布を設定する。
いま、表２のデータにおいて女性である（列１である）確率を
で表し、２項分布を適用する。すなわち、
とおく。
上記のモデルによるベイズ分析を行うStanスクリプトをリスト１のように用意した。
リスト１のStanスクリプトを用いてベイズ分析を行うスクリプトをリスト２のように用意した。ファイルは、c10files.zipにまとめた。
リスト１のスクリプトファイル（xtable.stan）とリスト２のスクリプトファイルを同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移す。CmdStanPyのインストールされている環境において、次のコマンドを実行する。CmdStanPyについては、簡単な説明をこのウェブサイトに用意した。
(stan) *****/c10files$
python ba_xtable.py
上のコマンドを実行すると、Stanスクリプトのコンパイルが始まる。このコンパイルは、多少時間が掛かる。コンパイルが終了すると、MCMCサンプリングが始まり、終了後、トレース図が表示される（図１）。
図１
左側に事後分布のグラフ、右側にサンプリングのトレースが描かれている。サンプリングに特に問題は無いようである。
図１のWindowを閉じると、図２のグラフが表示される。
図２
女性の条件（列１）において左腕が上（行１）である確率、男性（列２）において左腕が上（行１）である確率の事後分布のグラフである。ほぼ重なっており、差はなさそうである。
図２のWindowを閉じると、図３のグラフが表示される。
図３
確率
と
の差
の事後分布のグラフである。0より大きい事後確率が65.8%、0より小さい事後確率が34.2％であるので、は０である、とに差はないと考えられる。すなわち、どちらの腕が上であるかの確率に、男女差は認められないと言える。
図３のWindowを閉じると、図４のグラフが表示される。
図４
パラメータの事後分布のグラフである。パラメータは、表２のデータにおいて、サンプルが女性である確率を表している。が0.5より小さい事後確率が100.0％であるので、女性の方が少数であると言える。
図４のWindowを閉じると、スクリプトの実行終了である。
リスト１ 表１のデータに対するStanスクリプト（xtable.stan）
data {
/*
o11 o12
o21 o22
*/
int o11;
int o12;
int o21;
int o22;
}
parameters {
real<lower=0, upper=1>
pfm;
//
(o11+o21)/(o11+o21+o12+o22)
real<lower=0, upper=1>
p1; // o11/(o11 + o21)
real<lower=0, upper=1>
p2; // o12/(o12 + o22)
}
model {
o11 + o21 ~
binomial(o11+o21+o12+o22, pfm);
o11 ~ binomial(o11+o21, p1);
o12 ~ binomial(o12+o22, p2);
}
generated quantities {
real diff_p;
diff_p
= p2 - p1;
}
リスト２ リスト１のStanスクリプトによりベイズ分析を行るPythonスクリプト
from cmdstanpy
import CmdStanModel
import numpy
as np
import pandas as pd
import arviz
as az
import matplotlib.pyplot
as plt
import seaborn as sb
#
# Spiegelhalter(2019), Table
10.2
#
Female Male
data = [[5, 17], # Left arm on top
[9,
23]] # Right arm on top
Data = {'o11':data[0][0],
'o12':data[0][1],
'o21':data[1][0], 'o22':data[1][1]}
model = CmdStanModel(stan_file='xtable.stan')
fit = model.sample(data=Data)
print(fit.summary())
infdata = az.from_cmdstanpy(fit) # Arviz InfereceData型へ変換
summary = az.summary(infdata)
print('\nSummary...\n',
summary)
az.plot_trace(infdata)
plt.tight_layout()
plt.show()
dframe = fit.draws_pd()
# Pandas DataFrame型へ変換
sb.kdeplot(dframe['p1'],
color='b', label='p1')
sb.kdeplot(dframe['p2'],
color='g', label='p2')
plt.xlabel('p1, p2', fontsize=14)
plt.legend(fontsize=16)
plt.tight_layout()
plt.show()
diff_p = np.array(dframe['diff_p'])
L_diff_p = diff_p <=
0.0
print('L_diff_p
=',np.mean(L_diff_p))
R_diff_p = diff_p
>= 0.0
print('R_diff_p
=', np.mean(R_diff_p))
az.plot_kde(dframe['diff_p'])
plt.vlines([0.0], 0.0, 3.5)
plt.text(0.01, 3.3, 'p2$-$p1$=0$', fontsize=16)
plt.hlines([0.0], -0.4, 0.4)
plt.text(-0.152, 0.21, f'{np.mean(L_diff_p)*100:.1f}%',
va='bottom', fontsize=16)
plt.text(0.05, 0.21, f'{np.mean(R_diff_p)*100:.1f}%',
va='bottom', fontsize=16)
plt.xlabel('p2$-$p1', fontsize=16)
plt.tight_layout()
plt.show()
pfm = np.array(dframe['pfm'])
L_pfm = pfm <=
0.5
print('L_pfm
=', np.mean(L_pfm))
R_pfm = pfm >=
0.5
print('R_pfm
=', np.mean(R_pfm))
az.plot_kde(dframe['pfm'])
plt.vlines([0.5], 0.0, 3.5)
plt.hlines([0.0], 0.0, 1.0)
plt.text(0.2, 0.21, f'{np.mean(L_pfm)*100:.1f}%',
va='bottom', fontsize=16)
plt.text(0.51, 0.21, f'{np.mean(R_pfm)*100:.1f}%',
va='bottom', fontsize=16)
plt.text(0.5, 3.5, 'pfm=0.5',
fontsize=16)
plt.xticks(np.arange(0,
1.01, 0.1))
plt.xlabel('pfm', fontsize=16)
plt.tight_layout()
plt.show()
参考文献
岡本安晴（2019） いまさら聞けないPythonでデータ分析．丸善出版
Spiegelhalter, D. (2019). The art of statistics: Learning from data. Pelican