比率のクロス表のベイズ分析
ロジスティック回帰分析
2025.12：CmdStanPy対応に改訂
ベイズ分析ソフトStanのPythonインターフェースCmdStanPyの簡単な解説を、このウェブサイトに用意した。
比率のクロス表にロジスティック回帰モデルを適用してベイズ分析を試みる。
京都大学のホームページに掲載されている平成19年度と平成29年度の一般入試諸統計から合格者が200名以上の学部の合格者数内数（「認」と「他」を除く）を男女別にまとめると図１のようになった。
図１ Excelデータ（data1929.xlsx）
女子の比率をロジスティック回帰モデルで分析してみる。ロジスティック回帰モデルについては、＜岡本安晴「心理学データ分析と測定」勁草書房＞第４章で説明している。
表１のように各条件（行i、列j）における注目している事象の確率（図１の場合、女子の確率）をpijで表す。
|
表１ 行条件iと列条件jにおける確率pij |
||||
|
|
列１ |
・・・ |
列ｊ |
・・・ |
|
行１ |
p11 |
・・・ |
|
・・・ |
|
|
|
|
|
|
|
行i |
|
・・・ |
|
・・・ |
|
|
|
|
|
|
行条件iと列条件jの組み合わせにおける下段の度数（図１では女子の数）が総数（図１では女子の数と男子の数の和）の２項分布における度数として確率pijの推定を行う。確率pijは、比率「女子の数／（女子の数＋男子の数）」に対応するが、ベイズ分析を行うために、総数「女子の数＋男子の数」に対する該当事象の度数が「女子の数」である２項分布を適用してpijの分析を行う。
線形回帰モデルを用いるために、確率pijの値域［０、１］をロジット変換した値gijについて考える（Kruschke,2015）。
である。
ロジット変換値gijを行主効果ai、列主効果bj、交互作用(ab)jkの和に分解する（Myers & Well, 2003）。
ここで、制約条件を
と設定する。条件（３）は、モデル（２）を次式のモデルと区別するために必要である。
すなわち、
を主効果
と
、および交互作用
に分解する
このとき、
である。KRとKCは行条件の数と列条件の数を表す。
ロジスティック回帰モデルでは、式（１）より
となる。
各セルにおける度数は、確率pijの２項分布に従うとして、本ウェブサイト後半に示すリスト１のようにStanスクリプトを用意した。
リスト１では、
はデータの収集された共通状況（上の場合、京都大学）に対するパラメータ値であるので、共通の分布からのサンプルであるとして、次のように設定している。
mu ~ normal(0.0, 1000);
// Hyper prior distribution
sgm
~ uniform(0.0001, 1000); // Hyper prior distribution
for (i
in 1:KR)
for
(j in 1:KC) {
g[i][j] ~ normal(mu, sgm);
// Prior distribution
すなわち、は、ハイパーパラメータ平均mu、標準偏差sgmの正規分布をハイパー事前分布とすると設定する。
である。
リスト１のStanスクリプトを利用して表１のベイズ分析を行うPythonスクリプトは、本ウェブサイト後半に示すリスト２として用意した。
MCMCサンプリングで得られたサンプル値から式（２）の主効果および交互作用は、リスト２の次のスクリプトにより求めている。
for ig in
range(n_samples):
gdd
= fit_g[:,:,ig].mean()
gdc
= fit_g[:,:,ig].mean(axis =
0)
grd
= fit_g[:,:,ig].mean(axis =
1)
r_e[i] = grd - gdd
# 行主効果
c_e[i] = gdc - gdd
# 列主効果
rc_e[i] = fit_g[:,:,ig] - grd.reshape((KR, 1)) - gdc + gdd # 交互作用
i
+= 1
リスト１のStanスクリプトファイル（LogisModel_h.stan）、リスト２のPythonスクリプトファイル（MainLogistic.py）、および入力データファイル（例えば、図１のファイルであればdata1929.xlsx）を同じフォルダに置く。そのフォルダにカレントディレクトリ（フォルダ）を移し、CmdStanPyのインストールされた環境において次のコマンドを実行する。サンプルファイルは、filesLogisticReg.zipにまとめている。CmdStanPyの簡単な説明を、このウェブサイトに用意した。
(stan) *****/filesLogisticReg$ python MainLogistic.py
上のコマンドを実行すると、入力データファイル名の設定を求めるダイアログウィンドウが表示される（図２）。
図２
入力データファイルは、図１の形式で用意する。第1行目に、列条件のカテゴリ（図１の場合は、学部名）を書き、第１列目に行条件のカテゴリ（図１では年度）、第２列目に度数のカテゴリを書く。度数データは２行目以降、３列目から該当するセルに書く。各行条件・列条件に対応するセルは、上に非該当カテゴリ（図１では「男子」、２項分布の説明では「０」で表されることが多い）の度数、下の段に該当カテゴリ（図１では「女子」、２項分布の説明では「１」で表されることが多い）の度数を書く。列条件あるいは行条件は図１のものより多くてもよい。列および行条件の個数は２以上であればよい。
図２のダイアログウィンドウにファイル名を設定して「Open」ボタンをクリックすると、画像保存用ファイル名の設定を求めるダイアログウィンドウが表示される（図３）。
図３
グラフを保存するファイル名であるが、ファイル拡張子が「.png」である任意のファイル名を設定する。設定したファイル名の「.png」の部分がグラフの種類に合わせて変更された名前のファイルに画像が保存される。
図３のダイアログウィンドウにファイル名を設定して「Save」ボタンをクリックすると、計算結果を保存するためのテキストファイル名の設定を求めるダイアログウィンドウが表示される（図４）。
図４
ファイル拡張子が「.txt」の任意のテキストファイル名を設定する。
テキストファイル名を設定して「Save」ボタンをクリックすると、データの読み込みが行われ、Stanスクリプトのコンパイルが始まる。コンパイルには、多少時間が掛かる。っコンパイル後、MCMCサンプリングが行われ、サンプリング後、トレース図が表示される（図５）。
図５
左側に事後分布のグラフ、右側にサンプリングのトレース図が表示されている。
図５のWindowを閉じると、図６のグラフが表示される。
図６
式（９）のハイパー事前分布のグラフである。MCMCサンプルから100対のサンプル値（mu, sgm）をランダムに選び、それぞれのサンプル値（mu, sgm）に対して式（９）の正規分布のグラフを描いたものである。それら100本のグラフの平均値が赤の太い曲線で描かれている。曲線が負の領域に偏っていることから、パラメータgijは負の値をとる傾向があり、その値の式（８）による変換で得られるpijは0.5より小さい値となる傾向があることが分かる。
図６のWindowを閉じると、図７のグラフが表示される。
図７
主効果と交互作用の事後分布のグラフである。
左上のグラフから、平成19年度に比べて平成29年度は女子学生の比率が小さくなっていることが分かる。
右上のグラフから、学部別の女子学生の比率は、文学部と農学部が高く、理学部と工学部が低いことが分かる。法学部はその中間であるが、文学部と農学部の方に寄っている。
下段のグラフは交互作用であるが、次のグラフ（図８）を見ると分かりやすい。
図７のWindowを閉じると、図８のグラフが表示される。
図８
女子学生の確率pijの事後分布の中央値（小円盤）と95％CI（縦棒）のグラフである。Gelmanら（2021）は、事後分布の代表値として、平均値より安定した中央値を勧めている。縦棒は、各学部において左側が平成19年度、右側が平成29年度の確率の95％CIを表している。女子学生の確率は、平成19年度に比べて平成29年度は、文学部と農学部では減少しているが、法学部は増加している。しかし、この傾向は、95％CIの重なりを考えると断定はできない。
図８のWindowを閉じると、スクリプトの実行終了である。
実行終了後、図４で設定した出力ファイルを開くと以下のような内容である。
入力データファイル =
/mnt/c/WorkSpace/www/LaCoocan/Python/IntroPyStat/BysLogistic/filesLogisticReg/data1929.xlsx
行変数：['平成19年', '平成29年']
列変数：['文学部' '法学部' '理学部' '工学部' '農学部']
P = 女子 / (男子 + 女子)
Mean
MCSE StdDev ... ESS_tail ESS_bulk/s R_hat
lp__ -1830.460000 0.062440 2.557980 ... 2485.74 19736.8 1.001590
mu
-1.281640 0.004991 0.359806 ... 2642.32 66082.7 1.000770
sgm
1.056660 0.004967 0.311100 ... 2867.65 63336.6 1.000300
g[1,1] -0.068234 0.001502 0.134713 ... 3131.72 96478.3 1.001090
g[1,2] -1.273730 0.001439 0.134287 ... 2888.44 103848.0 1.000040
g[1,3] -1.985940 0.001855 0.174678 ... 2303.21 106691.0 0.999950
g[1,4] -2.398590 0.001299 0.115508 ... 3003.73 93471.1 1.002400
g[1,5] -0.479677 0.001189 0.113372 ... 3134.21 107971.0 1.001530
g[2,1] -0.383138 0.001475 0.138089 ... 3153.41 105571.0 1.002770
g[2,2] -0.958273 0.001414 0.126421 ... 2711.67 97203.8 1.000310
g[2,3] -2.142360 0.002030 0.180464 ... 2923.28 97600.0 1.000270
g[2,4] -2.385060 0.001301 0.116504 ... 3003.85 95477.3 1.000380
g[2,5] -0.700610 0.001320 0.122751 ... 2530.64 105540.0 1.000470
p[1,1] 0.483024 0.000373 0.033480 ... 3131.72 96478.4 1.001080
p[1,2] 0.219482 0.000245 0.022911 ... 2888.44 103848.0 1.000080
p[1,3] 0.121911 0.000197 0.018623 ... 2303.21 106691.0 0.999802
p[1,4] 0.083704 0.000099 0.008824 ... 3003.73 93470.9 1.002430
p[1,5] 0.382684 0.000281 0.026697 ... 3134.21 107971.0 1.001540
p[2,1] 0.405802 0.000354 0.033151 ... 3153.41 105571.0 1.002780
p[2,2] 0.277935 0.000280 0.025256 ... 2711.67 97203.8 1.000320
p[2,3] 0.106252 0.000195 0.017065 ... 2923.28 97600.0 1.000440
p[2,4] 0.084753 0.000101 0.009005 ... 3003.85 95477.3 1.000390
p[2,5] 0.332237 0.000291 0.027127 ... 2530.64 105540.0 1.000440
[23 rows x 11 columns]
行主効果
平成19年
平均値 = 0.036
中央値 = 0.037
Q1 = 0.007, Q3 = 0.066
95%確信区間 = [-0.048, 0.118]
平成29年
平均値 = -0.036
中央値 = -0.037
Q1 =
-0.066, Q3 = -0.007
95%確信区間 = [-0.118, 0.048]
列主効果
文学部
平均値 = 1.052
中央値 = 1.053
Q1 = 0.994, Q3 = 1.112
95%確信区間 = [0.876, 1.225]
法学部
平均値 = 0.162
中央値 = 0.161
Q1 = 0.103, Q3 = 0.220
95%確信区間 = [-0.007, 0.334]
理学部
平均値 = -0.787
中央値 = -0.785
Q1 =
-0.856, Q3 = -0.716
95%確信区間 = [-1.000, -0.580]
工学部
平均値 = -1.114
中央値 = -1.113
Q1 =
-1.167, Q3 = -1.062
95%確信区間 = [-1.265, -0.970]
農学部
平均値 = 0.687
中央値 = 0.687
Q1 = 0.636, Q3 = 0.739
95%確信区間 = [0.536, 0.839]
交互作用
平成19年×文学部
平均値 = 0.121
中央値 = 0.121
Q1 = 0.064, Q3 = 0.180
95%確信区間 = [-0.045, 0.285]
平成19年×法学部
平均値 = -0.194
中央値 = -0.193
Q1 =
-0.250, Q3 = -0.139
95%確信区間 = [-0.361, -0.027]
平成19年×理学部
平均値 = 0.042
中央値 = 0.042
Q1 =
-0.027, Q3 = 0.114
95%確信区間 = [-0.162, 0.247]
平成19年×工学部
平均値 = -0.043
中央値 = -0.044
Q1 =
-0.096, Q3 = 0.010
95%確信区間 = [-0.192, 0.105]
平成19年×農学部
平均値 = 0.074
中央値 = 0.074
Q1 = 0.024, Q3 = 0.126
95%確信区間 = [-0.075, 0.221]
平成29年×文学部
平均値 = -0.121
中央値 = -0.121
Q1 =
-0.180, Q3 = -0.064
95%確信区間 = [-0.285, 0.045]
平成29年×法学部
平均値 = 0.194
中央値 = 0.193
Q1 = 0.139, Q3 = 0.250
95%確信区間 = [0.027, 0.361]
平成29年×理学部
平均値 = -0.042
中央値 = -0.042
Q1 =
-0.114, Q3 = 0.027
95%確信区間 = [-0.247, 0.162]
平成29年×工学部
平均値 = 0.043
中央値 = 0.044
Q1 =
-0.010, Q3 = 0.096
95%確信区間 = [-0.105, 0.192]
平成29年×農学部
平均値 = -0.074
中央値 = -0.074
Q1 =
-0.126, Q3 = -0.024
95%確信区間 = [-0.221, 0.075]
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.
Myers, J.L. and Well, A.D. (2003) Research design and statistical analysis, 2nd ed. Lawrence Erlbaum Associates, Publishers.
リスト１ ロジスティック回帰分析用Stanスクリプト（LogisModel_h.stan）
data {
int KR;
int KC;
array[KR, KC, 2] int X;
}
parameters {
real mu;
// Hyper parameter
real<lower=0.0001> sgm;
// Hyper parameter
array[KR, KC] real g;
}
transformed parameters {
array[KR, KC] real<lower
= 0.0, upper = 1.0> p;
for (i
in 1:KR)
for
(j in 1:KC)
p[i][j] = inv_logit(g[i][j]);
}
model {
mu ~ normal(0.0, 1000);
// Hyper prior distribution
sgm
~ uniform(0.0001, 1000); // Hyper prior distribution
for (i
in 1:KR)
for
(j in 1:KC) {
g[i][j] ~ normal(mu, sgm);
// Prior distribution
}
for (i
in 1:KR)
for
(j in 1:KC) {
X[i][j][2] ~ binomial(X[i][j][1]
+ X[i][j][2], p[i][j]);
}
}
リスト２ リスト１のStanスクリプトを用いるPythonスクリプト（MainLogistic.py）
#
# Yasuharu Okamoto, 2020.08, 2025.12
#
import tkinter
as tk
import tkinter.filedialog
as fd
import pandas as pd
import numpy
as np
import scipy.stats
as ss
from cmdstanpy
import CmdStanModel
import matplotlib.pyplot
as plt
import japanize_matplotlib
import seaborn as sb
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
print(data)
KC = len(data[0])
- 2
print(KC)
KR = (len(data)
- 1) // 2
print(KR)
NameC = data[0][2:]
NameR = []
for i in
range(KR):
NameR.append(data[2*i + 1][0])
print(NameC)
print(NameR)
print('\nP =
{0} / ({1} + {0})\n'.format(data[2][1], data[1][1]))
fout.write('\n行変数：{}\n'.format(NameR))
fout.write('列変数：{}\n'.format(NameC))
fout.write('\n\nP =
{0} / ({1} + {0})\n'.format(data[2][1], data[1][1]))
X = np.empty((KR,
KC, 2), dtype = int)
for i in
range(KR):
for j in range(KC):
X[i][j][0] = data[i*2 + 1][j + 2]
X[i][j][1] = data[i*2 + 2][j + 2]
for i in
range(KR):
for j in range(KC):
print('{0:>10d}'.format(X[i][j][0]), end =
'')
print()
for j in range(KC):
print('{0:>10d}'.format(X[i][j][1]), end =
'')
print()
StanData = {'KR': KR, 'KC': KC, 'X': X}
model = CmdStanModel(stan_file='LogisModel_h.stan')
fit = model.sample(data=StanData)
print(fit.summary())
fout.write(f'{fit.summary()}')
infdata = az.from_cmdstanpy(fit) # Arviz InfereceData型へ変換
az.plot_trace(infdata)
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()
idx_slct = rng.choice(len(fit['mu']), size=100, replace=False)
xcoord = np.linspace(-8,
5, 1000)
curve_list = []
for idx in idx_slct:
ycoord
= ss.norm.pdf(xcoord, loc=fit['mu'][idx], scale=fit['sgm'][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.legend()
plt.title('Hyper Normal Distributions', fontsize=16)
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_h_normal.png')
SavedFiles.append(foutnm[:-4]
+ '_h_normal.png')
plt.show()
#
# 主効果と交互作用
#
fit_g = [[]]*KR
for ir in
range(KR):
fit_g[ir] = [[]]*KC
for ic
in range(KC):
print(np.shape(fit[f'g[{ir+1},{ic+1}]']))
fit_g[ir][ic]
= fit[f'g[{ir+1},{ic+1}]']
fit_g = np.array(fit_g)
print(KR, KC)
print(np.shape(fit_g))
n_samples = len(fit_g[0][0])
print('n_samples
=', n_samples)
r_e = np.empty((n_samples, KR))
c_e = np.empty((n_samples, KC))
rc_e = np.empty((n_samples, KR, KC))
i = 0
for ig in
range(n_samples):
gdd
= fit_g[:,:,ig].mean()
gdc
= fit_g[:,:,ig].mean(axis =
0)
grd
= fit_g[:,:,ig].mean(axis =
1)
r_e[i] = grd - gdd
# 行主効果
c_e[i] = gdc - gdd
# 列主効果
rc_e[i] = fit_g[:,:,ig] - grd.reshape((KR, 1)) - gdc + gdd # 交互作用
i
+= 1
plt.figure(figsize=(12,8))
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.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}'.format(
NameR[i], NameC[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]
+ '_effect.png')
SavedFiles.append(foutnm[:-4]
+ '_effect.png')
plt.show()
p_med = np.empty((KR,
KC))
p_025 = np.empty((KR,
KC))
p_975 = np.empty((KR,
KC))
for i in
range(KR):
for j in range(KC):
p_025[i][j], p_med[i][j], p_975[i][j] = \
np.percentile(fit[f'p[{i+1},{j+1}]'],
[2.5, 50, 97.5])
plt.figure(figsize=(7,5))
plt.title(r'$p_{ij}$の中央値, 95%CI' +
'\n' + f'{NameR}', fontsize
= 16)
for j in range(KC):
for i
in range(KR):
plt.plot([j-0.2 + 0.1*i, j-0.2 +
0.1*i],
[p_025[i][j], p_975[i][j]],
lw=5)
plt.plot(j-0.2 + 0.1*i, p_med[i][j], 'o', markersize=10)
plt.ylabel(r'$\hat{p}_{ij}$')
plt.xticks(range(KC), NameC)
plt.legend()
plt.tight_layout()
plt.savefig(foutnm[:-4]
+ '_pij.png')
SavedFiles.append(foutnm[:-4]
+ '_pij.png')
plt.show()
fout.close()
print('\n保存ファイル：')
for nm in SavedFiles:
print(' ', nm)
print()