度数クロス表データのベイズ分析
階層的ポアッソン回帰モデル
2025.12：CmdStanPy対応に改訂
ベイズ分析のソフトStanのPythonインターフェースであるCmdStanPy対応のスクリプトを書いてみた。CmdStanPyの簡単な説明は、このウェブサイトに用意した。
クロス表の形で与えられた度数データを、階層的ポアッソン回帰モデルにより分析してみた。
図１のデータは、＜芝ら（編）「統計用語辞典」新曜社＞の分割表の例として掲載されているもの（p. 230）を入力したものである。
図１ Excelデータ（datap230.xlsx）
図１の形式の度数のクロス表データに対してポアッソン回帰モデルによるベイズ分析を行う。度数データをポアッソン分布で表すとき、パラメータ
の対数に対して線形モデルが適用される（Reich & Ghosh, 2019）。ポアッソン分布のパラメータの対数を、クロス表の主効果と交互作用の和によって表す（Kruschke, 2015）。
クロス表の分析モデルを一般的に説明するために、表１のように各条件（行i、列j）における効果を
で表す。
|
表１ 主効果と交互作用 |
||||
|
|
列１ |
・・・ |
列ｊ |
・・・ |
|
行１ |
|
・・・ |
|
・・・ |
|
|
|
|
|
|
|
行i |
|
・・・ |
|
・・・ |
|
|
|
|
|
|
以下の議論では、各セルの度数は確率的に独立であるとする。
表１における効果を行主効果rai、列主効果caj、交互作用rcaijの和として
と表す。ここで、制約条件
をおく。この条件は、モデル（１）を次のモデル
と区別するために必要である。
制約条件（２）より、以下の式が導かれる。
上の効果に対してポアッソン分布のパラメータを対数関数の逆関数である指数関数によって次のようにおく。
各セルにおける度数は、パラメータが
であるポアッソン分布に従うとして、本ウェブサイト後半に示すリスト１のようにStanスクリプトを用意した。
リスト１では、データX[i][j]に対してlmbd[i][j]を求め、
X[i][j]
~ poisson(lmbd[i][j])
このlmbd[i][j]に対してz[i][j]を求めている。
z[i][j] = log(lmbd[i][j]);
これは、式（７）の示す計算の方向とは逆であるが、MCMCサンプリングにおける発散を避けるための処置である。
リスト１では、事前分布のパラメータは、ハイパー事前分布に従うとしている。これは、表１のデータが共通の枠組みからのサンプルであると考え、共通の枠組みをハイパー事前分布として設定したものである。
vr
~ exponential(0.01);
// Hyper prior distribution
for (i
in 1:KR)
for
(j in 1:KC) {
lmbd[i][j] ~
exponential(r); // Prior distribution
}
なお、指数分布のパラメータが０になるのを避けるため、次式の変換を行っている。
r = 1.0e-9 + vr;
リスト１のStanスクリプトによりzijの事後分布が得られる。このzijの事後分布から、式（４）～（６）により主効果および効果量の事後分布を本ウェブサイト後半に示すリスト２における次のスクリプトによって求めている。
for i in
range(n_samples):
zdd
= fit_z[:,:,i].mean()
# g..
zdc
= fit_z[:,:,i].mean(axis =
0) # g.j
zrd
= fit_z[:,:,i].mean(axis =
1) # gi.
r_e[i] = zrd - zdd
# 行主効果
c_e[i] = zdc - zdd
# 列主効果
rc_e[i] = fit_z[:,:,i] - zrd.reshape((KR, 1)) - zdc + zdd # 交互作用
スクリプトファイルなどのファイルは、圧縮ファイルfilesPoisson.zipとしてまとめたので、ダウンロード展開すれば自由に利用できる。なお、Stan（PyStan2）のPythonにおける利用方法は、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞の第10章～第12章で説明している。
リスト１のStanスクリプトファイル（PoissonModel_h.stan）、リスト２のPythonスクリプトファイル（MainPoisson.py）、および入力データアイル（例えば図１のファイルの場合はdatap230.xlsx）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移し、CmdStanyのインストールされた環境で、次のコマンドを実行する。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
(stan) *****/filesPoisson$ python MainPoisson.py
上のコマンドを実行すると、入力データファイル名の設定を求めるダイアログウィンドウが表示される（図２）。
図２
入力データファイルは、図１の形式で用意する。第1行目に、列条件のカテゴリ（図１の場合は、「自宅生」、「寮・下宿生」）を書き、第１列目に行条件のカテゴリ（図１では「男」、「女」）を書く。度数データは２行目以降、２列目から該当するセルに書く。列条件あるいは行条件は図１のものより多くてもよい。列および行条件の個数は２以上であればよい。
図２のダイアログウィンドウにファイル名を設定して「Open」ボタンをクリックすると、画像保存用ファイル名の設定を求めるダイアログウィンドウが表示される（図３）。
図３
事後分布のグラフを保存するファイル名であるが、ファイル拡張子が「.png」である任意のファイル名を設定する。設定したファイル名の「.png」の部分が「_trace.png」などに変更された名前のファイルに画像が保存される。
図３のダイアログウィンドウにファイル名を設定して「Save」ボタンをクリックすると、計算結果を保存するためのテキストファイル名の設定を求めるダイアログウィンドウが表示される（図４）。
図４
ファイル拡張子が「.txt」の任意のテキストファイル名を設定する。
テキストファイル名を設定して「Save」ボタンをクリックすると、データの読み込みが行われ、Stanスクリプトのコンパイルが始まる。コンパイルには多少時間が掛かる。コンパイル後、MCMCサンプリングが始まり、MCMCサンプリングが終了すると、トレース図が表示される（図５
）。
図５
左側のグラフに事後分布のグラフ、右側のグラフにMCMCサンプリングのトレース図が表示されている。
図５のWindowを閉じると、図６のグラフが表示される。
図６
式（７）の
は、与えられたデータ全体、例えば図１のデータ、の共通の枠組みからのサンプルであると考え、共通の枠組みを指数分布で与える。
lmbd[i][j] ~
exponential(r);
図６の左側のグラフは、MCMCサンプルからランダムにｒを100個選び、各ｒに対する指数分布のグラフを描いたものである。100本のグラフの平均値を赤の太い曲線で表している。スクリプトを以下に示す。
idx_slct = rng.choice(len(fit['r']), size=100, replace=False)
xcoord = np.linspace(0,
10, 1000)
curve_list = []
for idx in idx_slct:
ycoord
= ss.expon.pdf(xcoord, scale=1.0/fit['r'][idx])
plt.plot(xcoord,ycoord, alpha=0.3)
curve_list.append(ycoord)
mn_curve = np.mean(curve_list, axis=0)
plt.plot(xcoord, mn_curve, c='r', label='Mean values')
plt.title('Hyper Exponential Distributions', fontsize=16)
plt.xlabel('Lmbd', fontsize=14)
plt.legend()
右側のグラフは、ｒの事後分布のグラフである。
図６のWindowを閉じると、図７のグラフが表示される。
図７
行主効果、列主効果、交互作用の事後分布のグラフである。
左上のグラフは、行主効果の事後分布のグラフである。グラフの重なりが少ないので、男の方が多いことがかなり確信をもっていえる。
右上のグラフは，列主効果の事後分布のグラフである。グラフの重なりがかなりあるので、自宅生の方が量・下宿性より多い傾向があると言える。
下のグラフは、交互作用である。４本の交互作用のグラフであるが、２本ずつが重なっているので、２本に見える。グラフの注釈に示されている値は、それぞれの分布の中央値である。Gelman et al. (2021)は、事後分布の代表値として、平均値より安定した中央値を勧めている。男×自宅生（－0.29）と女×寮・下宿性（－0.29）が少なく、男×寮・下宿性（0.29）と女×自宅生（0.29）が多く、この２グループの事後分布は、それぞれのグループにはっきりと分離されている。
図７のWindowを閉じると、図８のグラフが表示される。
図８
式（１）の
の事後分布の中央値のグラフである。それぞれのグラフにおける２本の直線が平行でないことから、図７の下段のグラフとともに交互作用が明瞭に示されている。
図８のWindowを閉じると、図９のグラフが表示される。
図９
次式による事後予測度数のヒストグラムである。
X[i][j] ~ poisson(lmbd[i][j])
分布のパラメータlmbd[i][j]の値のMCMCサンプルに対して、与えられたポアッソン分布から度数X[i][j]をランダムに選んでいる。スクリプトを以下に示す。
for ir in
range(KR):
for ic
in range(KC):
smpls = []
for vlmbd in fit[f'lmbd[{ir+1},{ic+1}]']:
v = ss.poisson.rvs(vlmbd)
smpls.append(v)
ax[ir][ic].hist(smpls)
ax[ir][ic].set_title(NameR[ir] + ', ' + NameC[ic], fontsize=16)
グラフは、横軸および縦軸ともそれぞれ４つのグラフ間で共通である。
「男、寮・下宿生」が多く、「女、寮・下宿生」が少ないことがよく表されている。
図９のWindowを閉じると、スクリプトの実行終了である。
実行終了後、図４で設定した出力テキストファイルを開くと、以下のような内容である。
入力データファイル = /mnt/c/WorkSpace/www/LaCoocan/Python/IntroPyStat/BysPoisson/filesPoisson/datap230.xlsx
行変数：['男',
'女']
列変数：['自宅生' '寮・下宿生']
Mean
MCSE StdDev MAD ...
ESS_bulk
ESS_tail
ESS_bulk/s R_hat
lp__
96.905100 0.037562 1.634580 1.438750 ... 1881.83 2664.10 35506.2 1.00352
vr
0.086271 0.000590 0.041543 0.037806 ... 4967.17 2833.99 93720.1 1.00012
lmbd[1,1] 14.745100 0.053322 3.652570 3.500920 ... 4355.54 2857.02 82179.9 1.00180
lmbd[1,2] 20.236300 0.059104 4.466870 4.425940 ... 5663.97 3194.37 106867.0 1.00053
lmbd[2,1] 16.675400 0.053153 3.989930 3.833710 ... 5652.67 2729.76 106654.0 1.00073
lmbd[2,2] 7.349250 0.038303 2.630500 2.584860 ... 4578.93 2724.85 86394.8 1.00065
r
0.086271 0.000590 0.041543 0.037806 ... 4967.17 2833.99 93720.1 1.00012
z[1,1] 2.659790 0.003972 0.252526 0.242828 ... 4355.54 2857.02 82180.1 1.00212
z[1,2] 2.982970 0.003008 0.223101 0.221331 ... 5663.96 3194.37 106867.0 1.00053
z[2,1] 2.785230 0.003267 0.241620 0.236107 ... 5652.67 2729.76 106654.0 1.00103
z[2,2] 1.929000 0.005624 0.370066 0.371464 ... 4578.93 2724.85 86394.8 1.00112
[11 rows x 11 columns]
行主効果
男
平均値 = 0.232
中央値 = 0.230
Q1 = 0.134, Q3 = 0.325
95%確信区間 = [-0.032, 0.518]
女
平均値 = -0.232
中央値 = -0.230
Q1 =
-0.325, Q3 = -0.134
95%確信区間 = [-0.518, 0.032]
列主効果
自宅生
平均値 = 0.133
中央値 = 0.131
Q1 = 0.048, Q3 = 0.220
95%確信区間 = [-0.131, 0.408]
寮・下宿生
平均値 = -0.133
中央値 = -0.131
Q1 =
-0.220, Q3 = -0.048
95%確信区間 = [-0.408, 0.131]
交互作用
男×自宅生
平均値 = -0.295
中央値 = -0.290
Q1 =
-0.386, Q3 = -0.198
95%確信区間 = [-0.570, -0.030]
男×寮・下宿生
平均値 = 0.295
中央値 = 0.290
Q1 = 0.198, Q3 = 0.386
95%確信区間 = [0.030, 0.570]
女×自宅生
平均値 = 0.295
中央値 = 0.290
Q1 = 0.198, Q3 = 0.386
95%確信区間 = [0.030, 0.570]
女×寮・下宿生
平均値 = -0.295
中央値 = -0.290
Q1 =
-0.386, Q3 = -0.198
95%確信区間 = [-0.570, -0.030]
zij=log(λij)の中央値:
z[男, 自宅生] = 2.667
z[男, 寮・下宿生] = 2.994
z[女, 自宅生] = 2.791
z[女, 寮・下宿生] = 1.952
Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University
Kruschke, J. K. (2015). Doing Bayesian data analysis, 2nd ed. Academic Press.
Reich, B. J. & Ghosh, S. K. (2019). Beysian statistical methods. CRC Press.
芝 祐順、渡部 洋、石塚智一（編）（1984）．統計用語辞典．新曜社．
リスト１ ポアッソン回帰分析用Stanスクリプト（PoissonModel_h.stan）
data {
int KR;
int KC;
array[KR, KC] int X;
}
parameters {
real<lower=0> vr;
array[KR, KC] real<lower
= 0.0> lmbd;
}
transformed parameters {
real r;
r = 1.0e-9 + vr;
}
model {
vr
~ exponential(0.01);
// Hyper prior distribution
for (i
in 1:KR)
for
(j in 1:KC) {
lmbd[i][j] ~
exponential(r); // Prior distribution
}
for (i
in 1:KR)
for
(j in 1:KC) {
X[i][j] ~ poisson(lmbd[i][j]);
}
}
generated quantities {
array[KR, KC] real z;
for (i
in 1:KR)
for (j
in 1:KC) {
z[i][j] = log(lmbd[i][j]);
}
}
リスト２ リスト１のStanスクリプトを用いるPythonスクリプト（MainPoisson.py）
#
#
Yasuharu Okamoto, 2020.08
#
import tkinter
as tk
import tkinter.filedialog
as fd
import pandas as pd
import numpy
as np
from cmdstanpy
import CmdStanModel
import matplotlib.pyplot
as plt
import japanize_matplotlib
import seaborn as sb
import scipy.stats
as ss
import arviz
as az
def askinfilename(msg
= None, types = [('', '*.*')]):
"""
入力用ファイル名の設定
"""
rt = tk.Tk()
rt.withdraw()
flnm
= fd.askopenfilename(title = msg, filetypes = types)
rt.destroy()
return flnm
def askoutfilename(msg
= None, types = [('', '*.*')], deftype = '.txt'):
"""
出力用ファイル名の設定
"""
rt = tk.Tk()
rt.withdraw()
flnm
= fd.asksaveasfilename(title = msg, filetypes =
types,
defaultextension = deftype)
rt.destroy()
return flnm
finnm = askinfilename('Input
Data File', types = [('Excel(*.xlsx)', '*.xlsx')])
foutnm = askoutfilename(msg
= 'Output Image File',
types = [('Image(*.png)', '*.png')], deftype = '.png')
ftxtoutnm = askoutfilename(msg
= 'Output Text File',
types = [('Text (*.txt)', '*.txt')],
deftype = '.txt')
fout = open(ftxtoutnm,
'w')
SavedFiles = [ftxtoutnm]
fout.write('入力データファイル = {}\n'.format(finnm))
xlsx_fl = pd.ExcelFile(finnm)
data_xlsx = pd.read_excel(xlsx_fl, header
= None)
data = data_xlsx.values
KC = len(data[0])
- 1
print('KC = ', KC)
KR = len(data)
- 1
print('KR = ', KR)
NameC = data[0][1:]
NameR = []
for i in
range(KR):
NameR.append(data[i + 1][0])
print('NameC:
', NameC)
print('NameR:
', NameR)
fout.write('\n行変数：{}\n'.format(NameR))
fout.write('列変数：{}\n'.format(NameC))
X = np.empty((KR,
KC), dtype = int)
for i in
range(KR):
for j in range(KC):
X[i][j] = data[i + 1][j + 1]
for i in
range(KR):
for j in range(KC):
print('{0:>10d}'.format(X[i][j]), end = '')
print()
StanData = {'KR': KR, 'KC': KC, 'X': X}
model = CmdStanModel(stan_file='PoissonModel_h.stan')
fit = model.sample(data=StanData)
print(fit.diagnose())
print(fit.summary())
fout.write(f'{fit.summary()}')
infdata = az.from_cmdstanpy(fit) # Arviz InfereceData型へ変換
az.plot_trace(infdata, var_names=['r', 'lmbd'])
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_trace.png')
SavedFiles.append(foutnm[:-4]
+ '_trace.png')
plt.show()
fit = fit.draws_pd()
# Pandas DataFrame型への変換
print(fit.keys())
rng = np.random.default_rng()
plt.figure(figsize=(14,5))
plt.subplot(121)
idx_slct = rng.choice(len(fit['r']), size=100, replace=False)
xcoord = np.linspace(0,
10, 1000)
curve_list = []
for idx in idx_slct:
ycoord
= ss.expon.pdf(xcoord, scale=1.0/fit['r'][idx])
plt.plot(xcoord,ycoord, alpha=0.3)
curve_list.append(ycoord)
mn_curve = np.mean(curve_list, axis=0)
plt.plot(xcoord, mn_curve, c='r', label='Mean values')
plt.title('Hyper Exponential Distributions', fontsize=16)
plt.xlabel('Lmbd', fontsize=14)
plt.legend()
plt.subplot(122)
sb.kdeplot(fit['r'])
plt.title('Posterior distribution of r', fontsize=16)
plt.xlabel('r', fontsize=16)
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_exp_r.png')
plt.show()
#
# 主効果と交互作用
#
plt.figure(figsize=(12,
6))
fit_z = [[]]*KR
for ir in
range(KR):
fit_z[ir] = [[]]*KC
for ic
in range(KC):
print(np.shape(fit[f'z[{ir+1},{ic+1}]']))
fit_z[ir][ic]
= fit[f'z[{ir+1},{ic+1}]']
fit_z = np.array(fit_z)
print(KR, KC)
print(np.shape(fit_z))
n_samples = len(fit_z[0][0])
r_e = np.empty((n_samples, KR))
c_e = np.empty((n_samples, KC))
rc_e = np.empty((n_samples, KR, KC))
for i in
range(n_samples):
zdd
= fit_z[:,:,i].mean()
# g..
zdc
= fit_z[:,:,i].mean(axis =
0) # g.j
zrd
= fit_z[:,:,i].mean(axis =
1) # gi.
r_e[i] = zrd - zdd
# 行主効果
c_e[i] = zdc - zdd
# 列主効果
rc_e[i] = fit_z[:,:,i] - zrd.reshape((KR, 1)) - zdc + zdd # 交互作用
plt.subplot(2,4,(1,2))
plt.title('行主効果', fontsize = 20)
fout.write('\n行主効果\n')
for i in
range(KR):
sb.kdeplot(r_e[:, i], label = NameR[i])
ra_mean
= r_e[:, i].mean()
ra_025p, ra_25p, ra_50p,
ra_75p, ra_975p = \
np.percentile(r_e[:, i], [2.5, 25, 50, 75, 97.5])
fout.write('\n {}\n'.format(NameR[i]))
fout.write(' 平均値 = {0:>.3f}\n'.format(ra_mean))
fout.write(' 中央値 = {0:>.3f}\n'.format(ra_50p))
fout.write(' Q1 =
{0:>.3f}, Q3 =
{1:>.3f}\n'.format(ra_25p, ra_75p))
fout.write(' 95%確信区間 = [{0:>.3f},
{1:>.3f}]\n'.format(ra_025p, ra_975p))
plt.yticks([])
plt.legend()
plt.subplot(2,4,(3,4))
plt.title('列主効果', fontsize = 20)
fout.write('\n列主効果\n')
for j in range(KC ):
sb.kdeplot(c_e[:, j], label = NameC[j])
ca_mean
= c_e[:, j].mean()
ca_025p, ca_25p, ca_50p,
ca_75p, ca_975p = \
np.percentile(c_e[:, j], [2.5,
25, 50, 75, 97.5])
fout.write('\n {}\n'.format(NameC[j]))
fout.write(' 平均値 = {0:>.3f}\n'.format(ca_mean))
fout.write(' 中央値 = {0:>.3f}\n'.format(ca_50p))
fout.write(' Q1 =
{0:>.3f}, Q3 =
{1:>.3f}\n'.format(ca_25p, ca_75p))
fout.write(' 95%確信区間 = [{0:>.3f},
{1:>.3f}]\n'.format(ca_025p, ca_975p))
plt.yticks([])
plt.legend()
plt.tight_layout()
plt.subplot(2,4,(6,7))
plt.title('交互作用', fontsize = 20)
fout.write('\n交互作用\n')
for i in
range(KR):
for j in range(KC):
sb.kdeplot(rc_e[:, i, j], label = r'{0}$\times${1}: {2:.2f}'.
format(
NameR[i], NameC[j], np.median(rc_e[:, i, j])))
ga_mean = rc_e[:, i, j].mean()
ga_025p, ga_25p, ga_50p, ga_75p, ga_975p = \
np.percentile(rc_e[:,
i, j], [2.5, 25, 50, 75, 97.5])
fout.write('\n {0}×{1}\n'.format(
NameR[i], NameC[j]))
fout.write(' 平均値 = {0:>.3f}\n'.format(ga_mean))
fout.write(' 中央値 = {0:>.3f}\n'.format(ga_50p))
fout.write(' Q1 =
{0:>.3f}, Q3 =
{1:>.3f}\n'.format(ga_25p, ga_75p))
fout.write(' 95%確信区間 = [{0:>.3f},
{1:>.3f}]\n'.format(ga_025p, ga_975p))
plt.yticks([])
plt.legend()
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_rca.png')
SavedFiles.append(foutnm[:-4]
+ '_rca.png')
plt.show()
plt.figure(figsize=(12,
6))
plt.subplot(121)
fout.write('\n\nzij=log(λij)の中央値:\n')
z_med = np.empty((KR,
KC))
for i in
range(KR):
for j in range(KC):
z_med[i][j] = np.percentile(fit[f'z[{i+1},{j+1}]'], 50)
fout.write('\nz[{0}, {1}] =
{2:.3f}\n'.
format(NameR[i], NameC[j], z_med[i][j]))
plt.title(r'$\hat{z}_{ij}$の中央値', fontsize
= 20)
for i in
range(KR):
plt.plot(range(KC),
z_med[i], label = NameR[i])
plt.ylabel(r'$\hat{z}_{ij}$',
fontsize=14)
plt.xticks(range(KC), NameC,
fontsize=16)
plt.legend()
plt.subplot(122)
plt.title(r'$\hat{z}_{ij}$の中央値', fontsize
= 20)
for j in range(KC):
plt.plot(range(KR),
z_med[:, j], label = NameC[j])
plt.ylabel(r'$\hat{z}_{ij}$',
fontsize=14)
plt.xticks(range(KR), NameR,
fontsize=16)
plt.legend()
plt.savefig(foutnm[:-4]
+ '_zrc.png')
SavedFiles.append(foutnm[:-4]
+ '_zrc.png')
plt.show()
fig, ax = plt.subplots(KR,
KC, sharex=True, sharey=True,figsize=(12,8))
for ir in
range(KR):
for ic
in range(KC):
smpls = []
for vlmbd in fit[f'lmbd[{ir+1},{ic+1}]']:
v = ss.poisson.rvs(vlmbd)
smpls.append(v)
ax[ir][ic].hist(smpls)
ax[ir][ic].set_title(NameR[ir] + ', ' + NameC[ic], fontsize=16)
s_size = len(fit['lmbd[1,1]'])
fig.suptitle(f'事後予測度数、 サンプルサイズ={s_size}',
fontsize=18)
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_pred.png')
SavedFiles.append(foutnm[:-4]
+ '_pred.png')
plt.show()
fout.close()
SavedFiles.append(ftxtoutnm)
print('\n保存ファイル：')
for nm in SavedFiles:
print(' ', nm)
print()