香川県の交通事故発生の要因を
交通違反件数を基に分析する

import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

FIG_DIR = 'html/figures'
DATA_B = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}', na=False)].copy()
df_b = df_b[df_b['地域コード'] != 'R00000'].copy()   # 全国行を除外
df_b['年度'] = df_b['年度'].astype(int)

# 原論文が比較した平成29年（2017年度）の総人口
pop17 = df_b[df_b['年度'] == 2017][['都道府県', '総人口']].copy()
pop17['総人口'] = pd.to_numeric(pop17['総人口'], errors='coerce')
pop17['人口順位'] = pop17['総人口'].rank(ascending=False).astype(int)
pop17 = pop17.sort_values('総人口', ascending=False).reset_index(drop=True)

print('=== SSDSE-B 2017年度（平成29年）の総人口 ===')
for ken in ['香川県', '和歌山県']:
    row = pop17[pop17['都道府県'] == ken].iloc[0]
    print(f"{ken}: {row['総人口']:,.0f} 人（全国 {row['人口順位']} 位）")
diff = abs(pop17.set_index('都道府県').loc['香川県', '総人口']
           - pop17.set_index('都道府県').loc['和歌山県', '総人口'])
print(f"人口差: {diff:,.0f} 人")

# 出典: e-Stat「統計でみる都道府県のすがた2019」（原論文 3.3節・図1の報告値）
# ※ SSDSE には交通事故件数が未収録のため、原論文の報告値をそのまま描く（再計算ではない）
jiko_10man = {'香川県': 633.5, '和歌山県': 274.2}   # 件 / 人口10万人
jiko_rank = {'香川県': '全国6位', '和歌山県': '全国36位'}

fig1, ax = plt.subplots(figsize=(6.5, 5))
colors = ['#C62828', '#1565C0']
bars = ax.bar(jiko_10man.keys(), jiko_10man.values(), color=colors, width=0.55)
for bar, (ken, v) in zip(bars, jiko_10man.items()):
    ax.text(bar.get_x() + bar.get_width()/2, v + 12,
            f'{v:.1f}件\n（{jiko_rank[ken]}）', ha='center', fontsize=12, fontweight='bold')
ax.set_ylabel('人口10万人あたり交通事故件数（件）')
ax.set_ylim(0, 780)
ax.set_title('人口10万人あたりの交通事故件数（原論文 図1の報告値）\n'
             '出典: e-Stat「統計でみる都道府県のすがた2019」', fontsize=12)
ax.text(0.5, -0.13, '※ 原論文の報告値の可視化（本ページでの再計算ではない）',
        transform=ax.transAxes, ha='center', fontsize=9, color='#777')
fig1.savefig(os.path.join(FIG_DIR, '2019_H5_2_fig1.png'), bbox_inches='tight')
plt.close(fig1)
print('図1 saved:', os.path.join(FIG_DIR, '2019_H5_2_fig1.png'))
print(f"香川県 {jiko_10man['香川県']}件 は 和歌山県 {jiko_10man['和歌山県']}件 の "
      f"約{jiko_10man['香川県']/jiko_10man['和歌山県']:.1f}倍")

# 原論文 図3「都道府県別人口（20位〜47位）」の再現
sub = pop17.iloc[19:].copy()          # 20位〜47位
sub = sub.sort_values('総人口')       # 少ない順に下から

bar_colors = ['#C62828' if k == '香川県' else '#1565C0' if k == '和歌山県'
              else '#B0BEC5' for k in sub['都道府県']]
fig2, ax = plt.subplots(figsize=(8, 9))
ax.barh(sub['都道府県'], sub['総人口'] / 1e4, color=bar_colors)
for y, (k, v, r) in enumerate(zip(sub['都道府県'], sub['総人口'], sub['人口順位'])):
    if k in ('香川県', '和歌山県'):
        ax.text(v / 1e4 + 2, y, f'{r}位 {v/1e4:.0f}万人',
                va='center', fontsize=10, fontweight='bold')
ax.set_xlabel('総人口（万人・2017年度）')
ax.set_title('都道府県別総人口 20位〜47位（原論文 図3の再現）\n'
             'SSDSE-B-2026 の実データから作成', fontsize=12)
fig2.savefig(os.path.join(FIG_DIR, '2019_H5_2_fig2.png'), bbox_inches='tight')
plt.close(fig2)
print('図2 saved:', os.path.join(FIG_DIR, '2019_H5_2_fig2.png'))

# 目的変数: 交通事故発生件数（平成14〜29年の年次データ）
# 説明変数: 各交通違反の取り締まり件数（警察庁データ）
# ※ 警察庁データは SSDSE 未収録のため係数は原論文の報告値（再計算ではない）
variables = ['速度違反', '飲酒運転', '整備不良', '一時停止', '駐停車']

kagawa = pd.DataFrame({
    '回帰係数': [-0.012, 2.134, 3.636, -0.517, -0.622],
    '標準誤差': [0.056, 0.964, 0.578, 0.072, 0.198],
    't値':     [-0.221, 2.214, 6.294, -7.205, -3.149],
    'p値':     [0.829, 0.051, 0.000, 0.000, 0.010],
}, index=variables)

wakayama = pd.DataFrame({
    '回帰係数': [0.112, 3.388, 3.447, -0.315, -0.845],
    '標準誤差': [0.067, 1.015, 0.662, 0.170, 0.320],
    't値':     [1.686, 3.337, 5.210, -1.851, -2.640],
    'p値':     [0.123, 0.008, 0.000, 0.094, 0.025],
}, index=variables)

print('=== 表1 香川県（原論文の報告値, R2=0.96, 補正R2=0.94）===')
print(kagawa)
print()
print('=== 表2 和歌山県（原論文の報告値, R2=0.91, 補正R2=0.865）===')
print(wakayama)
print()
print('p<0.05で有意なのに回帰係数が負の変数（香川）:',
      list(kagawa[(kagawa['p値'] < 0.05) & (kagawa['回帰係数'] < 0)].index))

def stars(p):
    return '**' if p < 0.01 else '*' if p < 0.05 else ''

y = np.arange(len(variables))
h = 0.38
fig3, ax = plt.subplots(figsize=(8.5, 5.5))
ax.barh(y + h/2, kagawa['回帰係数'], height=h, color='#C62828', label='香川県')
ax.barh(y - h/2, wakayama['回帰係数'], height=h, color='#1565C0', label='和歌山県')
for yi, v in enumerate(variables):
    ck, cw = kagawa.loc[v], wakayama.loc[v]
    ax.text(ck['回帰係数'] + (0.08 if ck['回帰係数'] >= 0 else -0.08), yi + h/2,
            f"{ck['回帰係数']:+.3f}{stars(ck['p値'])}", va='center',
            ha='left' if ck['回帰係数'] >= 0 else 'right', fontsize=10)
    ax.text(cw['回帰係数'] + (0.08 if cw['回帰係数'] >= 0 else -0.08), yi - h/2,
            f"{cw['回帰係数']:+.3f}{stars(cw['p値'])}", va='center',
            ha='left' if cw['回帰係数'] >= 0 else 'right', fontsize=10)
ax.axvline(0, color='#555', lw=1)
ax.set_yticks(y)
ax.set_yticklabels(variables)
ax.invert_yaxis()
ax.set_xlabel('回帰係数（目的変数: 交通事故発生件数）')
ax.set_xlim(-1.8, 4.8)
ax.legend(loc='lower right')
ax.set_title('重回帰の回帰係数の比較（原論文 表1・表2の報告値の可視化）\n'
             '* p<0.05, ** p<0.01 ／ 警察庁データはSSDSE未収録のため再計算ではない',
             fontsize=11)
fig3.savefig(os.path.join(FIG_DIR, '2019_H5_2_fig3.png'), bbox_inches='tight')
plt.close(fig3)
print('図3 saved:', os.path.join(FIG_DIR, '2019_H5_2_fig3.png'))

# 出典: グーネット「交通事故・交通死亡事故の発生原因は何？」（原論文が引用）
gense = {'安全不確認': 31, 'わき見運転': 16, '動静不注視': 11,
         '漠然運転': 9, '運転操作不適': 6, 'その他': 27}

labels = list(gense.keys())
vals = list(gense.values())
cols = ['#C62828' if k == '安全不確認' else '#90A4AE' if k == 'その他'
        else '#1565C0' for k in labels]

fig4, ax = plt.subplots(figsize=(8, 5))
ax.barh(labels, vals, color=cols)
for yi, v in enumerate(vals):
    ax.text(v + 0.5, yi, f'{v}%', va='center', fontsize=11, fontweight='bold')
ax.invert_yaxis()
ax.set_xlabel('交通事故件数に占める割合（%）')
ax.set_xlim(0, 36)
ax.set_title('法令違反別 交通事故件数の構成比（原論文 図16の報告値の可視化）\n'
             '出典: グーネット（原論文が引用）／「速度違反」は上位項目に現れない', fontsize=11)
fig4.savefig(os.path.join(FIG_DIR, '2019_H5_2_fig4.png'), bbox_inches='tight')
plt.close(fig4)
print('図4 saved:', os.path.join(FIG_DIR, '2019_H5_2_fig4.png'))
print('安全不確認（31%）がトップ。一時停止を怠る行動は「安全不確認」型の事故に直結する。')