交通事故と年齢のベイズ分析
クロス表の主効果の分析も行う
クロス表の分析は、カイ2乗検定が標準であるが、カイ２乗検定は交互作用を見るもので、主効果は無視されている。クロス表のデータを分散分析のデザインでモデル構成を行うことにより、交互作用を具体的に調べるとともに主効果も検討することができる（Kruschke, 2011）。
分析例として、年齢別の交通事故のデータを「警視庁の統計 令和5年（2023年）：https://www.keishicho.metro.tokyo.lg.jp/about_mpd/jokyo_tokei/tokei/k_tokei05.html」からダウンロードしたデータ（2025.05.25）を用いる。ダウンロードしたファイルkt16.xlsxから、分析に使用するデータを抜き出して図１のファイルを作成した。なお、本ウェブサイトのファイルは圧縮ファイルtrafficaccifiles.zipとした。
図１
分析結果が分かり易いものになるように、図１のデータを年齢によりグループ分けする。各年齢ごとに事故の比率を求め、内積表現（http://y-okamoto-psy1949.la.coocan.jp/Python/tips/CrossTableSVD/）を求め、図２の結果を得た。
図２
図２の布置から、70歳以上のグループと70歳未満のグループに分けることが考えられる。
内積表現のクラスター分析の結果（図3）からも、70歳以上と70歳未満にグループ分けすることが考えられる。
図３
図１のデータを年齢で２グループに分けてまとめたものが図４のファイルである。
図４
図４のデータは分散分析のデザインで分析できる（Kruschke, 2011）。
クロス表（図４）におけるセル
の値
がパラメータ
のポアッソン分布に従うとする。
パラメータの値は、分散分析の考え方で表１のように与える（岡本、2019、 pp. 202-203）。
|
表１ クロス表のモデル： |
|||
|
|
死者 |
重傷者 |
軽症者 |
|
70歳未満 |
|
|
|
|
70歳以上 |
|
|
|
「70歳未満、死者」を基準としてパラメータab0を与える。「70歳未満」に対する「70以上」の主効果をa1で表す。「死者」に対する「重傷者」、「軽症者」の主効果をb1、b2で表す。これらの主効果に対して、「70以上×重傷者」および「70歳以上×軽症者」の交互作用をc11およびc12で表す。
上のモデルに基づくベイズ分析のStanスクリプトをリスト１に用意した。リスト１のStanスクリプトを用いて図４のデータのベイズ分析を行うPythonスクリプトをリスト２のように作成した。ファイルは、圧縮ファイルtrafficaccifiles.zipにまとめた。
リスト１のスクリプトでは、対象となる年齢群の母集団人口に対するパラメータ
の調整を次のコードで行っている。
x[i][j] ~ poisson(n_pop[i]*lmbd[i][j]);
n_pop[1]は70歳未満の総人口、n_pop[2]は70歳以上の総人口であり、lmbd[i][j]は1人当たりの事故発生率を表す。
リスト２を実行すると、先ず、主効果aの事後分布のグラフが表示される（図５）。
図５
主効果aは、「70歳未満、死者」ab0を基準にしている。図５の結果は、a1が0より大きいので、「70歳以上、死者」の1人当たりの度数がより多いことを示している。
事故の内容による主効果bのグラフは、図６のようになる。
図６
軽症の方が重傷者より人口1人当たりの度数が多いことが分かる。また死者を基準にした効果であるので、重傷者と軽症者の分布が正の領域にあることから、これらの1人当たりの度数は死傷者より多いことが分かる。事後分布が分離されており、かつ、いずれの分布の領域も０を含まないので、これらの差は統計学的に明確な差であると考えられる。
図７は交互作用cの事後分布である。
図７
「70歳以上の軽症者」の事後分布が負の領域にあることから、70以上の軽症者の1人当たりの度数が少ないことが分かる。70歳以上の重傷者の事後分布も負の領域に寄っていることから、70歳以上の重傷者の1人当たりの度数も小さい傾向があると考えられる。これらは交互作用であるので、主効果と比べて1人当たりの度数が少ないことを示していることに注意。
図８にab0+b1（70歳未満、重傷者）の事後分布とab0+a1+b1+c11（70歳以上、重傷者）の事後分布を示してみた。
図８
「70歳以上、重傷者」の事後分布の方が右に位置している。「70歳未満、重傷者」と「70歳以上、重傷者」を主効果も含めて考えた場合、「70歳以上、重傷者」の方が1人当たりの事故頻度が高いということが明瞭に示されている。
Kruschke, J. K. (2011) Doing Bayesian data analysis. Elsevier.
岡本安晴（2019） いまさら聞けないPythonでデータ分析．丸善出版
リスト１ 表１のモデルのStanスクリプト（poissonmodel.stan）
data {
int n_col;
int n_row;
array[n_row,
n_col] int x;
array[n_row]
real n_pop;
}
parameters {
real ab0;
array[n_col-1] real b;
array[n_row-1] real a;
array[n_row-1, n_col-1] real
c;
}
transformed parameters {
array[n_row,
n_col] real lmbd;
for (i
in 1:n_row)
for
(j in 1:n_col) {
if (i == 1) {
if (j == 1) {
lmbd[i][j] =
exp(ab0);
}
else {
lmbd[i][j] = exp(ab0
+ b[j-1]);
}
}
else {
if (j == 1) {
lmbd[i][j] = exp(ab0
+ a[i-1]);
}
else {
lmbd[i][j] = exp(ab0
+ a[i-1] + b[j-1] + c[i-1][j-1]);
}
}
}
}
model {
ab0 ~ normal(0, 10);
a ~ normal(0, 10);
b ~ normal(0, 10);
for (i
in 1:n_row-1)
for
(j in 1:n_col-1)
c[i][j] ~ normal(0, 10);
for (i
in 1:n_row)
for
(j in 1:n_col) {
x[i][j] ~ poisson(n_pop[i]*lmbd[i][j]);
}
}
リスト２ 図４のデータの分析スクリプト（AnalTraffic.py）
from cmdstanpy
import CmdStanModel
import pandas as pd
import numpy
as np
import matplotlib.pyplot
as plt
import arviz
as az
import japanize_matplotlib
import seaborn as sb
data_frame = pd.read_excel('table_2groups.xlsx')
print(data_frame)
col_names = data_frame.keys()[1:]
row_names = data_frame['年齢'].values
print('row_names
=\n', row_names)
print('col_names
=\n', col_names)
x = data_frame.values[:,1:]
print(x)
n_pops = [10781259, 2479294]
Data = {'n_row':len(row_names), 'n_col':len(col_names), 'x':x, 'n_pop':n_pops}
sm = CmdStanModel(stan_file='poissonmodel.stan')
fit = sm.sample(data=Data)
print(fit.summary())
i_data = az.from_cmdstanpy(fit)
smmry = az.summary(i_data)
pd.set_option('display.max_rows',
None)
print(smmry)
fit_df = fit.draws_pd()
ab0s = fit_df['ab0']
a = {}
for i in
range(len(row_names)-1):
a[i+1] = fit_df[f'a[{i+1}]']
print(a.keys())
print(f'ab0-{row_names[0]}
= {np.median(ab0s):.3f}')
for i in
range(len(row_names)-1):
sb.kdeplot(a[i+1],
label=f'a{i+1}:{row_names[i+1]}')
print(f'a[{i+1}]-{row_names[i+1]} = {np.median(a[i+1]):.3f}')
plt.xlabel('main effect a', fontsize=14)
plt.ylabel('Density', fontsize=14)
plt.legend()
plt.show()
b = {}
for j in range(len(col_names)-1):
b[j+1] = fit_df[f'b[{j+1}]']
for j in range(len(col_names)-1):
sb.kdeplot(b[j+1],
label=f'b{j+1}:{col_names[j+1]}')
plt.xlabel('main effect b', fontsize=14)
plt.ylabel('Density', fontsize=14)
plt.legend()
plt.show()
c = {}
for i in
range(len(row_names)-1):
c[i+1] = {}
for j in range(len(col_names)-1):
c[i+1][j+1] = fit_df[f'c[{i+1},{j+1}]']
sb.kdeplot(c[i+1][j+1], label=f'c[{row_names[i+1]}][{col_names[j+1]}]')
plt.legend()
plt.xlabel('Interaction c', fontsize=14)
plt.ylabel('Density', fontsize=14)
plt.show()
sb.kdeplot(ab0s+b[1], label='70歳未満、重傷者')
sb.kdeplot(ab0s+a[1]+b[1]+c[1][1], label='70歳以上、重傷者')
plt.title('ab0+b[1] vs ab0+a[1]+b[1]+c[1][1]', fontsize=14)
plt.legend()
plt.show()