腕組みの性差の2項分布による階層的ロジスティック回帰モデルのベイズ分析
クロス表のベイズ分析が提案されている（Kruscke, 2011；岡本、2019）。クロス表の度数分布に対してポアッソン分布が適用されているが、Spiegelhalter (2019)のデータ（Table 10.2）のようなクロス表の場合は階層的２項分が適当であると考えられる。Spiegelhalterのデータでは、男女別に腕組みが右腕が上か左腕が上かの度数が示されている。これは受講生におけるデータであるので、自然発生的度数と考えるよりは、男女の度数を２項分布で表し、さらに男女別にどちらの腕が上であるかの度数を２項分布で表すことが考えられる。このとき、「男女別」および「右腕が上／左腕が上」の２つの主効果と、「男女別」と「右腕が上／左腕が上」の交互作用を設定したロジスティック回帰モデルによる階層的2項分布モデルを以下のように考える。
一般的に考えるために、男女の要因をA（A1：女性、A2：男性）、腕組みをした時のどちらの腕が上であったかを要因B（B1：左が上、B2：右が上）で表すと、クロス表は表１の形に書くことができる。
|
Table 1. Observed counts of levels of Factor A by levels of Factor B. |
||||
|
|
|
Factor B |
|
|
|
|
|
B1 |
B2 |
Sum |
|
Factor A |
A1 |
|
|
|
|
A2 |
|
|
|
|
|
Sum |
|
|
|
|
Spiegelhalter (2019)のデータ（Table 10.2）における上記Table 1の値は、
5, 
9
7, 
23
である。
要因Aについて、２項分布を適用して
と置く。
要因Aの各水準において、要因Bについて２項分布を適用して
と置く。要因Bの主効果は
で表され、
は交互作用である。
ロジスティック関数
の逆関数
は、図１のようである。
図１
上記モデルによるTable 1のベイズ分析を行うスクリプトを、リスト１およびリスト２に示す。リスト１とリスト２に分割したのは、Windows上では、PyStan 3を用いるリスト１のスクリプトは、WindowsにインストールしたUbuntu上で実行され、リスト１の実行により得られたMCMCサンプルの描画などはWindows上でリスト２のスクリプトにより行われるからである。スクリプトファイルは、まとめて圧縮ファイルhibinfiles.zipとした。
リスト１を実行すると、（D1）と（D2）に示されるデータに対してMCMCサンプリングが行われ、サンプルがファイルに保存される。リスト２を実行すると、保存されたMCMCサンプルが読み込まれて、まず要因Aのパラメータaの事後分布が表示される（図2）。
図２
事後分布は負の領域にシフトしており、
でパラメータaが正である事後確率は0.1％と小さい。統計学的検定の慣例を参考にすれば、パラメータは負であると考えてよさそうである。すなわち、女性の比率は0.5以下と解釈できる。
図２のWindowを閉じると、図３のグラフが表示される。
図３
要因Ｂの主効果のパラメータbの事後分布である。事後分布は負の方向にシフトしているが、
であるので、統計学的検定の慣例を参考にすれば、ｂは負であるとは断言できないが、bが負である事後確率は0.863とかなり高い。女性において左腕が上という傾向があると解釈できる。
図３のWindowを閉じると、図４のグラフが表示される。
図４
交互作用のパラメータcの事後分布である。
であり、原点０の左右にほぼ同じ程度に分布している。交互作用は無いと解釈できる。すなわち、左腕が上であるという傾向は、男女間で同じと解釈できる。
図４のWindowを閉じると、図５のグラフが表示される。
図５
女性における左腕が上のパラメータbと、男性における左腕が上のパラメータb+cの事後分布のグラフである。２つの事後分布はほぼ重なっている。
Kruschke, J. K. (2011) Doing Bayesian data analysis. Elsevier.
岡本安晴（2019） いまさら聞けないPythonでデータ分析．丸善出版
Spiegelhalter, D. (2019) The art of statistics. Penguin Random House
リスト１ MCMCサンプリング： Python3.12/Miniconda/Ubuntu24/WSL/Windows10で実行した。
import stan
import arviz as az
import pickle
"""
D.
Spiegelhalter (2019)
The art of
statistics: Learnig from data.
Penguin
Random House
Table 10.2
"""
# Left arm Right arm
data = [
[ 5,
9 ], # Female
[ 17,
23 ] ] #
Male
stan_code = """
data
{
array[2, 2] int F;
}
transformed data {
int row1_sum;
int row2_sum;
int total;
row1_sum = F[1][1] + F[1][2];
row2_sum = F[2][1] + F[2][2];
total = sum(F[1]) + sum(F[2]);
}
parameters {
real a;
real b;
real c;
}
model {
a ~ normal(0.0, 10.0);
b ~ normal(0.0, 10.0);
c ~ normal(0.0, 10.0);
row1_sum ~ binomial(total, inv_logit(a));
F[1][1] ~ binomial(row1_sum, inv_logit(b));
F[2][1] ~ binomial(row2_sum, inv_logit(b + c));
}
"""
Data = {'F':data}
sm = stan.build(stan_code, data=Data)
fit = sm.sample()
print(fit)
"""
MCMCサンプルの保存
"""
df_fit = fit.to_frame()
with open('logit_fit.pkl', 'wb') as f:
pickle.dump(df_fit, f)
"""
MCMCサンプルの統計量
"""
inf_data =
az.from_pystan(posterior=fit, posterior_model=sm)
smmry = az.summary(inf_data)
print(smmry)
リスト２ MCMCサンプルの描画： Python3.13/Windows10で実行した。
import arviz as az
import pickle
import matplotlib.pyplot as plt
import seaborn as sb
import numpy as np
"""
MCMCサンプルの読み込み
"""
with open('Logit_fit.pkl', 'rb') as f:
df_fit = pickle.load(f)
p_a = np.mean(df_fit['a'] > 0)
plt.title(f'Prob(a > 0) =
{p_a:.3f}', fontsize=16)
sb.kdeplot(df_fit['a'])
plt.xlabel('a', fontsize=16)
plt.ylabel('Density', fontsize=16)
plt.show()
p_b = np.mean(df_fit['b'] > 0)
plt.title(f'Prob(b > 0) =
{p_b:.3f}', fontsize= 16)
sb.kdeplot(df_fit['b'])
plt.xlabel('b', fontsize=16)
plt.ylabel('Density', fontsize=16)
plt.show()
p_c = np.mean(df_fit['c'] > 0)
plt.title(f'Prob(c > 0) =
{p_c:.3f}', fontsize=16)
sb.kdeplot(df_fit['c'])
plt.xlabel('c', fontsize=16)
plt.ylabel('Density', fontsize=16)
plt.show()
sb.kdeplot(df_fit['b'], label='b')
sb.kdeplot(df_fit['b'] + df_fit['c'],
label='b+c')
plt.xlabel('')
plt.ylabel('Density', fontsize=16)
plt.legend(fontsize=18)
plt.show()