Up

交通事故と年齢のベイズ分析

 

70歳以上の人の交通事故のパターンが他の年齢の人とは異なることが示された。また、クロス表の検定はカイ二乗検定が標準であるが、ベイズ分析を用いることにより、より豊かな情報を得ることができた。

年齢別の交通事故のデータを「警視庁の統計　令和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歳以上

 

上のモデルに基づくベイズ分析のスクリプトをリスト１およびリスト２に用意した。

リスト１はPython3.12/miniconda/Ubuntu24/WSL/Windows10で実行し、リスト２はPython3.13/Windows10で実行した。スクリプトファイルとデータファイルは同じフォルダに置き、リスト１の実行後、リスト２を実行した。使用したファイルは、圧縮ファイルtrafficaccifiles.zipにまとめた。

 

リスト１の実行により得られたMCMCサンプルの分析をリスト２のスクリプトで行っている。

リスト１のスクリプトでは、対象となる年齢群の総人口に対するパラメータの調整を次のコードで行っている。

 

x[i][j] ~ poisson(n_pop[i]*lmbd[i][j]);

 

リスト２を実行すると、先ず、主効果aの事後分布のグラフが表示される（図５）。

 

図５

 

主効果aは、「70歳未満、死者」ab0を基準にしているので、図５の結果は、「70歳以、死者」の度数がより多いことを示している。

各年齢群の交通事故死の確率を図４の表より求めると、図６のようになる。

 

図６

 

70歳以上の事故死の確率が高いことが分かる。

事故の内容による主効果bのグラフは、図７のようになる。

 

図７

 

軽症の方が重傷者より多いことが分かる。また死者を基準にした効果であるので、重傷者と軽症者の分布が正の領域にあることから、これらの確率は死傷者より高いことが分かる。事後分布が分離されており、かつ、いずれの分布の領域も０を含まないので、これらの差は統計学的に意味のある差であると考えられる。

図８は交互作用cの事後分布である。

 

図８

 

70歳以上の軽症者の事後分布が負の領域にあることから、70以上の軽症者の確率が小さいことが分かる。70歳以上の重傷者の事後分布も負の領域に寄っていることから、70歳以上の重傷者の比率も小さい傾向があると考えられる。

年齢群別に事故者の比率を求め、棒グラフを描くと図９のようになる。

 

図９

 

図８では、70歳以上、重傷者の事後分布が負の領域寄りであるが、図９では70歳以上の重傷者の比率は70歳未満の重傷者の比率より高い。70歳以上と70歳未満を比較するときは、図５に示されている年齢の主効果a1を踏まえて考える必要がある。図10にab0+b1（70歳未満、重傷者）の事後分布とa0+a1+b1+c11（70歳以上、重傷者）の事後分布を示す。70歳以上、重傷者の事後分布の方が右に位置している。

 

図10

 

 

 

 

 

引用文献

 

Kruschke, J. K. (2011) Doing Bayesian data analysis. Elsevier.

岡本安晴（2019）　いまさら聞けないPythonでデータ分析．丸善出版

 

 

リスト１　Ubuntu24/WSL/Windows10上で実行した。

 

import stan

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

import arviz as az

import pickle

 

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)

 

stan_code = """

    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]);

            }

    }

"""

 

n_pops = [10781259, 2479294]

 

Data = {'n_row':len(row_names), 'n_col':len(col_names), 'x':x, 'n_pop':n_pops}

sm = stan.build(stan_code, data=Data)

fit = sm.sample()

print(fit)

i_data = az.from_pystan(posterior=fit, posterior_model=sm)

smmry = az.summary(i_data)

 

pd.set_option('display.max_rows', None)

print(smmry)

 

fit_df = fit.to_frame()

print(fit)

 

with open('fit.pkl', 'wb') as f:

    pickle.dump(fit_df, f)

 

 

リスト２　リスト１の実行後、Python 3.13/Windows10で実行した

 

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

import japanize_matplotlib

import seaborn as sb

import pickle

 

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)

 

with open('fit.pkl', 'rb') as f:

    fit_df = pickle.load(f)

print(fit_df)

 

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()

 

 

 

 

Home