ポアッソン分布によるクロス表のベイズ分析
主効果と交互作用
岡本安晴
2026.07.07改訂
クロス表の標準的分析法であるカイ二乗検定では、交互作用の情報のみが得られる。クロス表には主効果という情報も含まれているが、これを無視することはクロス表に含まれる重要な情報の一部を無視することになる。主効果と交互作用の情報をともに扱う分析の枠組みとして分散分析がよく知られている。
表１の形式のデータ（岡本、2019、表8.2.2より）に対して、分散分析のデザインでポアッソン分布を用いたベイズ分析を行ってみた。
|
表１ |
|
|
|
|
|
|
心理学≒臨床心理学 |
その他 |
心理学＞臨床心理学 |
計 |
|
小学生 |
1 |
1 |
10 |
12 |
|
中学生 |
3 |
2 |
10 |
15 |
|
高校生 |
6 |
2 |
4 |
12 |
|
計 |
10 |
5 |
24 |
39 |
いま、一般的に考えて、行数がNR、列数がNCのクロス表を考える。
第
行、第
列の度数を
で表す。
度数の分布をポアッソン分布で表し（Spiegelhalter, 2019）、
とおく。
上のモデルのStanスクリプトをリスト１のように用意した。ファイルは、poissonfiles.zipにまとめた。
リスト１ ポアッソン分布のStanスクリプト（xtable_poisson.stan）
data {
int nr;
int nc;
array[nr, nc] int F;
}
parameters {
array[nr, nc]
real<lower=0> lmbd;
}
model {
for (ir in 1:nr) {
for
(ic in 1:nc) {
lmbd[ir][ic] ~ uniform(0.0001, 100);
F[ir][ic] ~ poisson(lmbd[ir][ic]);
}
}
}
パラメータ
に対して、主効果、交互作用を以下のように与える。
パラメータ全体の平均を
とおく。
行主効果を
とおき、列主効果を
とおく。
交互作用を次式で与える。
したがって、
である。
上のモデルによるクロス表の分析を行うPythonスクリプトをリスト２のように作成した。ファイルは、poissonfiles.zipにまとめた。
リスト１のStanスクリプトファイル（xtable_poisson.stan）とリスト２のPythonスクリプトファイル（xtablePoisson.py）、および入力データファイルを同じフォルダに置き、次のPythonコマンドを実行する。
(stan) ****/poissonfiles$
python xtablePoisson.py
Data
file(*.xlsx) = sample.xlsx
入力データファイル名を聞いてくるので、上の例ではsample.xlsxを設定している。
ファイルsample.xlsは、表１のデータのExcelファイルである（図０）。
図０ 表１のExcelファイル
１行目と１列目に行カテゴリ名と列カテゴリ名をおき、２行２列目以降に度数データを置いている。
ファイル名が設定されると、次のスクリプトによりデータが読み込まれる。
flnm = input('Data file(*.xlsx) = ')
pd_data = pd.read_excel(flnm)
F = pd_data.values[:,1:]
print(F)
NR, NC = np.shape(F)
print(NR,NC)
ファイルが読み込まれと、次のスクリプトでStanスクリプトのコンパイルとMCMCサンプリングが実行される。
model =
CmdStanModel(stan_file="xtable_poisson.stan") # コンパイル
fit = model.sample(data={'nr':NR,
'nc':NC, 'F':F}) # サンプリング
MCMCサンプリングが終了すると、次のスクリプトでトレース図（図１）が表示される。
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
# トレースプロット
az.plot_trace_dist(inf_data)
plt.tight_layout()
plt.show()
図１
図１のWindowの右上のX印をクリックして閉じると、計算が次に進む。
次のスクリプトにより、MCMCサンプルがPandasのDataFrame型として取り出され、配列Lmbdに格納される。
df_sample = fit.draws_pd() # PandasのDataFrame型に変換
Lmbd = []
for ir in range(NR):
t_lmbd = []
for ic in range(NC):
t_lmbd.append(df_sample[f"lmbd[{ir+1},{ic+1}]"])
Lmbd.append(t_lmbd)
Lmbd = np.array(Lmbd)
この配列から主効果と交互作用のサンプル値が次のスクリプトにより算出される。
g = []
a = [] # row effects
b = [] # column effects
ab = [] # interaction
for i in range(len(Lmbd[0,0,:])):
v = Lmbd[:,:,i]
g.append(np.mean(v))
va = np.mean(v, axis=1)
a.append(va - g[-1])
vb = np.mean(v, axis=0)
b.append(vb - g[-1])
t_v0 = []
for ir in range(NR):
t_v
= []
for
ic in range(NC):
t_v.append(v[ir,ic] - va[ir] - vb[ic] + g[-1])
t_v0.append(t_v)
ab.append(t_v0)
g = np.array(g)
a = np.array(a)
b = np.array(b)
vab = np.array(ab)
a = a.T
b = b.T
ab = np.empty((NR, NC, len(g)))
for ir in range(NR):
for ic in range(NC):
for
i in range(len(g)):
ab[ir][ic][i] = vab[i][ir][ic]
主効果の事後分布が、次のスクリプトにより表示される（図２）。
plt.figure(figsize=(12,5))
plt.subplot(121)
for ir in range(NR):
sb.kdeplot(a[ir],
label=f"a{ir+1}")
plt.xlabel('a', fontsize=14)
plt.legend()
plt.tight_layout()
plt.subplot(122)
for ic in range(NC):
sb.kdeplot(b[ic],
label=f"b{ic+1}")
plt.xlabel('b', fontsize=14)
plt.legend()
plt.tight_layout()
plt.show()
図２
左のグラフに行主効果aの事後分布が描かれている。３つの分布はかなり重なっていて実質的差は認められない。行は、心理学という分野があることを知った時期であるが、表１に記されているデータでは、3つの時期、小学生、中学生、高校生、に統計学的差は認められないということになる。なお、その他に回答した回答者は０人であった。
右のグラフは、列主効果bの事後分布である。b3の事後分布が他の２つより右にある。列カテゴリは、大学受験のときに、心理学≒臨床心理学と思っていた、心理学とは何をする分野なのか知らなかった、心理学には臨床心理学以外の分野があることを知っていたであるが、心理学には臨床心理学以外の分野があることを知っていたと答えた回答者が他の回答者より統計学的に多いことが示されている。ちなみに、回答者は某女子大学心理学科の１年生である（2015年頃のアンケート調査）。
交互作用の事後分布は、次のスクリプトで箱図として描いた（図３）。
xdata = {}
for ir in range(NR):
for ic in range(NC):
xdata[f"a{ir+1}b{ic+1}"] = ab[ir][ic]
plt.plot([-1, len(xdata)], [0, 0],
ls='--')
sb.boxplot(xdata)
plt.title('Boxplots of Posterior
Distributions', fontsize=16)
plt.tight_layout()
plt.show()
図３
箱が０より離れている交互作用は、a1b1、a1b3、a3b1、a3b3の４つである。交互作用
は、主効果
で説明されないものを表すものである（式（1））。
交互作用a1b1の事後分布が負の領域にあることは、「小学生」×「心理学≒臨床心理学」の度数が主効果で説明されるよりも少ないことを表している。
交互作用a1b3が正の領域にあることは、「小学生」×「心理学＞臨床心理学」の度数が、主効果で説明されるよりも多いことを示している。
交互作用a3b1が正の領域にあることは、「高校生」×「心理学≒臨床心理学」の度数が、主効果で説明されるよりも多いことを示している。
交互作用a3b3が負の領域にあることは、「高校生」×「心理学＞臨床心理学」の度数が、主効果で説明されるよりも少ないことを表している。
交互作用を検討することにより、主効果で表される傾向からのずれの大きいものがわかる。「心理学という分野を知った時期が小学生のとき」という回答者では「心理学は臨床心理学を含む分野である」ことを理解している回答者が多くなる傾向があるが、「心理学という分野を知った時期が高校生のとき」という回答者では、「心理学≒臨床心理学」と考えていた回答者が多くなる傾向にあると考えられる。
図３のWindowを閉じると、スクリプトの実行終了である。
なお、現在（2026.07.06）、matplotlib 3.11を用いているとエラーが出る。Matplotlib 3.10なら大丈夫であった。コマンド
conda install
matplotlib=3.10
を実行すると、matplotlib 3.10の最新版に置き換えられる。
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, 3rd. Ed. CRC Press.
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University
Kirk, R.E. (1995). Experimental design: Procedures for the behavioral sciences, third ed. Brooks/Cole Publishing Company.
Myers, J. L. & Well, A. D. (2003). Research design and statistical analysis, second ed. Lawrence Erlbaum Associates, Publishers.
岡本安晴（2019） いまさら聞けないPythonでデータ分析．丸善出版
Spiegelhalter, D. (2019). The art of statistics: Learning from data. Pelican
Winer, B. J., Brown, D.R., & Michels, K.M. (1991). Statistical principles in experimental design, third ed. McGraw-Hill, Inc.
リスト２ クロス表のベイズ分析Pythonスクリプト（xtablePoisson.py）
from cmdstanpy import CmdStanModel
import pandas as pd
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import arviz as az
import seaborn as sb
flnm = input('Data file(*.xlsx) = ')
pd_data = pd.read_excel(flnm)
F = pd_data.values[:,1:]
print(F)
NR, NC = np.shape(F)
print(NR,NC)
model =
CmdStanModel(stan_file="xtable_poisson.stan") # コンパイル
fit = model.sample(data={'nr':NR,
'nc':NC, 'F':F}) # サンプリング
print(fit.summary())
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
# トレースプロット
az.plot_trace_dist(inf_data)
plt.tight_layout()
plt.show()
df_sample = fit.draws_pd() # PandasのDataFrame型に変換
Lmbd = []
for ir in range(NR):
t_lmbd = []
for ic in range(NC):
t_lmbd.append(df_sample[f"lmbd[{ir+1},{ic+1}]"])
Lmbd.append(t_lmbd)
Lmbd = np.array(Lmbd)
g = []
a = [] # row effects
b = [] # column effects
ab = [] # interaction
for i in range(len(Lmbd[0,0,:])):
v = Lmbd[:,:,i]
g.append(np.mean(v))
va = np.mean(v, axis=1)
a.append(va - g[-1])
vb = np.mean(v, axis=0)
b.append(vb - g[-1])
t_v0 = []
for ir in range(NR):
t_v
= []
for
ic in range(NC):
t_v.append(v[ir,ic] - va[ir] - vb[ic] + g[-1])
t_v0.append(t_v)
ab.append(t_v0)
g = np.array(g)
a = np.array(a)
b = np.array(b)
vab = np.array(ab)
a = a.T
b = b.T
ab = np.empty((NR, NC, len(g)))
for ir in range(NR):
for ic in range(NC):
for
i in range(len(g)):
ab[ir][ic][i] = vab[i][ir][ic]
plt.figure(figsize=(12,5))
plt.subplot(121)
for ir in range(NR):
sb.kdeplot(a[ir],
label=f"a{ir+1}")
plt.xlabel('a', fontsize=14)
plt.legend()
plt.tight_layout()
plt.subplot(122)
for ic in range(NC):
sb.kdeplot(b[ic],
label=f"b{ic+1}")
plt.xlabel('b', fontsize=14)
plt.legend()
plt.tight_layout()
plt.show()
xdata = {}
for ir in range(NR):
for ic in range(NC):
xdata[f"a{ir+1}b{ic+1}"] = ab[ir][ic]
plt.figure(figsize=(12,5))
plt.plot([-1, len(xdata)], [0, 0],
ls='--')
sb.boxplot(xdata)
plt.title('Boxplots of Posterior
Distributions', fontsize=16)
plt.tight_layout()
plt.show()