医師数偏在の要因と
医師数偏在の診療科偏在への影響

print("=" * 65)
print("■ データ読み込み・構築（東京都除く46都道府県）")
print("=" * 65)

PREFS_46 = [
    '北海道','青森県','岩手県','宮城県','秋田県','山形県','福島県',
    '茨城県','栃木県','群馬県','埼玉県','千葉県','神奈川県',
    '新潟県','富山県','石川県','福井県','山梨県','長野県','岐阜県',
    '静岡県','愛知県','三重県','滋賀県','京都府','大阪府','兵庫県',
    '奈良県','和歌山県','鳥取県','島根県','岡山県','広島県','山口県',
    '徳島県','香川県','愛媛県','高知県','福岡県','佐賀県','長崎県',
    '熊本県','大分県','宮崎県','鹿児島県','沖縄県'
]

# ── SSDSE-B-2026 読み込み（2022年度）──────────────────────────
df_b_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'),
    encoding='cp932', header=1
)
df_b_raw['年度'] = pd.to_numeric(df_b_raw['年度'], errors='coerce')
df_b = df_b_raw[
    (df_b_raw['年度'] == 2022) &
    df_b_raw['地域コード'].str.match(r'^R\d{5}$', na=False)
].copy()
df_b = df_b[df_b['都道府県'].isin(PREFS_46)].set_index('都道府県')

for col in ['総人口', '65歳以上人口', '15～64歳人口', '死亡数',
            '月間有効求人数（一般）', '月間有効求職者数（一般）',
            '消費支出（二人以上の世帯）', '大学数', '一般病院数', '一般診療所数',
            '年平均気温']:
    df_b[col] = pd.to_numeric(df_b[col], errors='coerce')

# ── SSDSE-E-2026 読み込み（横断面データ）──────────────────────
df_e_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-E-2026.csv'),
    encoding='cp932', header=1
)
df_e = df_e_raw.iloc[1:].copy()
df_e.columns = df_e_raw.iloc[0].values
df_e = df_e[df_e['都道府県'] != '全国'].reset_index(drop=True)
df_e = df_e[df_e['都道府県'].isin(PREFS_46)].set_index('都道府県')

for col in ['医師数', '総人口', '65歳以上人口', '一般病院数', '一般診療所数',
            '歯科診療所数', '1人当たり県民所得（平成27年基準）',
            '総面積（北方地域及び竹島を除く）', '消費支出（二人以上の世帯）']:
    df_e[col] = pd.to_numeric(df_e[col], errors='coerce')

# ── 主要指標の計算 ────────────────────────────────────────────
df = pd.DataFrame({'pref': PREFS_46})
df = df.set_index('pref')

# 医師数_10万対（人口10万人あたり医師数） — SSDSE-E 医師数 / SSDSE-E 総人口
df['医師数_10万対'] = (
    df_e['医師数'] / df_e['総人口'] * 100000
)

# 高齢化率 = 65歳以上人口 / 総人口 × 100
df['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100

# 粗死亡率 = 死亡数 / 総人口 × 1000（人口千対）
df['粗死亡率'] = df_b['死亡数'] / df_b['総人口'] * 1000

# 有効求人倍率（代理） = 月間有効求人数 / 月間有効求職者数
df['有効求人倍率'] = (
    df_b['月間有効求人数（一般）'] / df_b['月間有効求職者数（一般）']
)

# 1人当たり県民所得（SSDSE-E）
df['1人当たり県民所得'] = df_e['1人当たり県民所得（平成27年基準）']

# 年平均気温（SSDSE-B）
df['年平均気温'] = df_b['年平均気温']

# 面積当たり病院診療所数（SSDSE-E: 面積はha単位 → 100km² 換算）
# 総面積(ha) / 100 = km² × 0.01 → 万km² = 総面積(ha) / 1,000,000
# 病院診療所数 per 100km² = (一般病院数 + 一般診療所数) / (総面積ha / 10000)
df['面積100km2当たり病院診療所数'] = (
    (df_e['一般病院数'] + df_e['一般診療所数']) /
    (df_e['総面積（北方地域及び竹島を除く）'] / 10000)
)

# 大学数（SSDSE-B 2022）
df['大学数'] = df_b['大学数']

# 消費支出（SSDSE-E）
df['消費支出'] = df_e['消費支出（二人以上の世帯）']

df = df.dropna(subset=['医師数_10万対', '高齢化率'])
df = df.reset_index()

N = len(df)
print(f"分析対象: {N}都道府県")
print("\n記述統計（主要変数）:")
print(df[['医師数_10万対', '高齢化率', '粗死亡率', '年平均気温',
          '面積100km2当たり病院診療所数']].describe().round(3))

print("\n" + "=" * 65)
print("■ Step3. Shapiro-Wilk 検定（正規性） + Levene 検定（等分散性）")
print("=" * 65)

# 医師数_10万対の3群分類 → 粗死亡率を比較
q33 = df['医師数_10万対'].quantile(1/3)
q67 = df['医師数_10万対'].quantile(2/3)
df['医師数群'] = pd.cut(df['医師数_10万対'],
                        bins=[-np.inf, q33, q67, np.inf],
                        labels=['低医師群（医師少）', '中医師群', '高医師群（医師多）'])

groups = {g: df[df['医師数群'] == g]['粗死亡率'].dropna().values
          for g in ['低医師群（医師少）', '中医師群', '高医師群（医師多）']}

print("\nShapiro-Wilk 検定（各群の粗死亡率の正規性）:")
for gname, gvals in groups.items():
    sw_stat, sw_p = stats.shapiro(gvals)
    normal = '正規分布に従う' if sw_p > 0.05 else '正規分布に従わない'
    print(f"  {gname}: W={sw_stat:.4f}, p={sw_p:.4f} → {normal}")

lev_stat, lev_p = stats.levene(*groups.values())
print(f"\nLevene 検定（等分散性）: F={lev_stat:.4f}, p={lev_p:.4f}")
print(f"  → {'等分散が仮定できる' if lev_p > 0.05 else '等分散は仮定できない → Kruskal-Wallis 検定を使用'}")

print("図1: 相関分析を作成中...")

fig1, axes1 = plt.subplots(1, 2, figsize=(14, 6))
fig1.suptitle('医師数_10万対の相関分析', fontsize=13, fontweight='bold')

# (a) 相関係数棒グラフ
ax1a = axes1[0]
corr_vars_plot  = [CORR_LABELS.get(v, v) for v in CORR_VARS]
corr_vals_plot  = [corr_results[v][0] for v in CORR_VARS]
corr_pvals_plot = [corr_results[v][1] for v in CORR_VARS]
bar_colors_corr = ['#E53935' if r < 0 else '#43A047' for r in corr_vals_plot]
sorted_idx_corr = np.argsort(corr_vals_plot)
ax1a.barh(np.arange(len(CORR_VARS)),
          [corr_vals_plot[i] for i in sorted_idx_corr],
          color=[bar_colors_corr[i] for i in sorted_idx_corr],
          alpha=0.8, edgecolor='white')
ax1a.set_yticks(np.arange(len(CORR_VARS)))
ax1a.set_yticklabels([corr_vars_plot[i] for i in sorted_idx_corr], fontsize=9)
ax1a.axvline(0, color='black', linewidth=1.0)
ax1a.set_xlabel('Pearson 相関係数（医師数_10万対との相関）', fontsize=10)
ax1a.set_title('医師数_10万対との相関係数\n（赤=負の相関、緑=正の相関）', fontsize=10, fontweight='bold')
ax1a.grid(axis='x', alpha=0.3)
for i, (idx, r, p) in enumerate(zip(sorted_idx_corr,
                                    [corr_vals_plot[j] for j in sorted_idx_corr],
                                    [corr_pvals_plot[j] for j in sorted_idx_corr])):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    if sig:
        offset = 0.02 if r > 0 else -0.02
        ha = 'left' if r > 0 else 'right'
        ax1a.text(r + offset, i, sig, va='center', ha=ha, fontsize=9, fontweight='bold')

# (b) 医師数_10万対 vs 粗死亡率 散布図
ax1b = axes1[1]
sc = ax1b.scatter(df['医師数_10万対'], df['粗死亡率'],
                  c=df['高齢化率'], cmap='RdYlGn_r',
                  s=80, alpha=0.85, edgecolors='white', linewidth=0.5, zorder=3)
x_s = df['医師数_10万対'].dropna().values
y_s = df['粗死亡率'].loc[df['医師数_10万対'].notna()].values
slope, intercept, r_val, p_val, _ = stats.linregress(x_s, y_s)
x_line = np.linspace(x_s.min(), x_s.max(), 100)
ax1b.plot(x_line, intercept + slope * x_line, '-', color='#333', linewidth=2.0,
          label=f'回帰直線 r={r_val:.3f} (p={p_val:.3f})')

# 注目都道府県
highlight_prefs = ['京都府', '秋田県', '沖縄県']
for pref in highlight_prefs:
    row = df[df['pref'] == pref]
    if len(row) > 0:
        ax1b.annotate(pref.replace('府', '').replace('県', ''),
                     (row['医師数_10万対'].values[0], row['粗死亡率'].values[0]),
                     fontsize=8.5, fontweight='bold', color='#333',
                     xytext=(6, 4), textcoords='offset points',
                     arrowprops=dict(arrowstyle='->', color='#555', lw=0.8))

plt.colorbar(sc, ax=ax1b, label='高齢化率（%）')
ax1b.set_xlabel('医師数_10万対（人口10万対）', fontsize=11)
ax1b.set_ylabel('粗死亡率（人口千対）', fontsize=11)
ax1b.set_title('医師数_10万対 vs 粗死亡率\n色：高齢化率（緑=低、赤=高）',
               fontsize=10, fontweight='bold')
ax1b.legend(fontsize=9)
ax1b.grid(True, alpha=0.3)

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2025_H1_fig1_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("  → 2025_H1_fig1_corr.png 保存完了")

print("図2: 重回帰係数プロットを作成中...")

fig2, axes2 = plt.subplots(1, 2, figsize=(14, 6))
fig2.suptitle('重回帰分析結果：医師数_10万対の規定要因（東京都除く46都道府県）',
              fontsize=13, fontweight='bold')

# (a) 回帰係数（森林プロット）
ax2a = axes2[0]
coefs_reg  = [reg_res.params.get(zv, 0) for zv in Z_REG]
ses_reg    = [reg_res.bse.get(zv, 0) for zv in Z_REG]
pvals_reg  = [reg_res.pvalues.get(zv, 1) for zv in Z_REG]
labels_reg = [CORR_LABELS.get(v.replace('_z', ''), v) for v in Z_REG]

bar_colors2 = ['#E53935' if c < 0 and p < 0.05
               else '#43A047' if c > 0 and p < 0.05
               else '#9E9E9E'
               for c, p in zip(coefs_reg, pvals_reg)]

ax2a.barh(np.arange(len(Z_REG)), coefs_reg,
          xerr=[1.96 * s for s in ses_reg],
          color=bar_colors2, alpha=0.8, edgecolor='white',
          capsize=4, error_kw={'elinewidth': 1.5, 'ecolor': '#555'})
ax2a.axvline(0, color='black', linewidth=1.0)
ax2a.set_yticks(np.arange(len(Z_REG)))
ax2a.set_yticklabels(labels_reg, fontsize=9)
ax2a.set_xlabel('標準化回帰係数（±95%CI）', fontsize=10)
ax2a.set_title(f'重回帰係数（HC1ロバスト標準誤差）\nR²={reg_res.rsquared:.3f}',
               fontsize=10, fontweight='bold')
ax2a.grid(axis='x', alpha=0.3)
ax2a.invert_yaxis()

for i, (c, s, p) in enumerate(zip(coefs_reg, ses_reg, pvals_reg)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    if sig:
        ax2a.text(c + np.sign(c) * 1.96 * s + np.sign(c) * 0.01, i,
                  sig, va='center', ha='left' if c > 0 else 'right',
                  fontsize=10, fontweight='bold')

# (b) VIF 棒グラフ
ax2b = axes2[1]
vif_labels = [CORR_LABELS.get(v.replace('_z', ''), v) for v in Z_REG]
vif_values_plot = [vif_vals[v] for v in Z_REG]
bar_colors_vif = ['#E53935' if v > 5 else '#FF8F00' if v > 3 else '#43A047'
                  for v in vif_values_plot]
ax2b.barh(np.arange(len(Z_REG)), vif_values_plot,
          color=bar_colors_vif, alpha=0.8, edgecolor='white')
ax2b.axvline(5, color='#E53935', linestyle='--', linewidth=1.5, label='VIF=5（警戒線）')
ax2b.axvline(3, color='#FF8F00', linestyle='--', linewidth=1.0, label='VIF=3（注意線）')
ax2b.set_yticks(np.arange(len(Z_REG)))
ax2b.set_yticklabels(vif_labels, fontsize=9)
ax2b.set_xlabel('VIF（分散拡大係数）', fontsize=10)
ax2b.set_title('多重共線性確認（VIF）\nVIF<5: 問題なし', fontsize=10, fontweight='bold')
ax2b.legend(fontsize=9)
ax2b.grid(axis='x', alpha=0.3)
ax2b.invert_yaxis()
for i, v in enumerate(vif_values_plot):
    ax2b.text(v + 0.05, i, f'{v:.2f}', va='center', fontsize=8.5)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2025_H1_fig2_regression.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("  → 2025_H1_fig2_regression.png 保存完了")

print("図3: Kruskal-Wallis検定結果を作成中...")

fig3, axes3 = plt.subplots(1, 2, figsize=(13, 6))
fig3.suptitle('Kruskal-Wallis 検定 + Dunn 検定（Bonferroni補正）\n粗死亡率と医師数水準の関係',
              fontsize=12, fontweight='bold')

# (a) ボックスプロット + ストリッププロット
ax3a = axes3[0]
group_data = [groups[g] for g in ['低医師群（医師少）', '中医師群', '高医師群（医師多）']]
bp = ax3a.boxplot(group_data, patch_artist=True, notch=False, widths=0.55)
bp_colors = ['#E53935', '#FF8F00', '#43A047']
for patch, color in zip(bp['boxes'], bp_colors):
    patch.set_facecolor(color)
    patch.set_alpha(0.7)
for median in bp['medians']:
    median.set(color='white', linewidth=2.5)

# ジッタープロット（決定論的：x位置を均等分散）
for i, (gdata, color) in enumerate(zip(group_data, bp_colors), 1):
    n = len(gdata)
    jitter = np.linspace(-0.15, 0.15, n)
    ax3a.scatter(i + jitter, gdata, color=color, alpha=0.6, s=35, zorder=3)

# 有意差の記号
y_max = max(max(g) for g in group_data) * 1.05
for (gname1, gname2, u, p_raw, p_bonf) in dunn_results:
    if p_bonf < 0.05:
        i1 = ['低医師群（医師少）', '中医師群', '高医師群（医師多）'].index(gname1) + 1
        i2 = ['低医師群（医師少）', '中医師群', '高医師群（医師多）'].index(gname2) + 1
        ax3a.plot([i1, i2], [y_max, y_max], 'k-', linewidth=1.5)
        sig_text = '***' if p_bonf < 0.001 else '**' if p_bonf < 0.01 else '*'
        ax3a.text((i1 + i2) / 2, y_max * 1.005, sig_text,
                  ha='center', fontsize=12, fontweight='bold')
        y_max *= 1.08

ax3a.set_xticklabels(['低医師群\n（医師少）', '中医師群', '高医師群\n（医師多）'], fontsize=10)
ax3a.set_ylabel('粗死亡率（人口千対）', fontsize=11)
ax3a.set_title(f'医師数3群の粗死亡率比較\nKW: H={kw_stat:.2f}, p={kw_p:.4f}',
               fontsize=10, fontweight='bold')
ax3a.grid(axis='y', alpha=0.3)

# (b) 前提確認検定サマリー
ax3b = axes3[1]
ax3b.axis('off')
table_data = [['前提確認検定', '検定統計量', 'p値', '判定']]
for gname, gvals in groups.items():
    sw_stat_, sw_p_ = stats.shapiro(gvals)
    judge = '正規性OK' if sw_p_ > 0.05 else '非正規性'
    short = gname.replace('医師群（医師少）', '(少)').replace('医師群（医師多）', '(多)').replace('医師群', '')
    table_data.append([f'Shapiro-Wilk\n{short}', f'W={sw_stat_:.4f}', f'p={sw_p_:.3f}', judge])
table_data.append([f'Levene 検定\n（等分散性）', f'F={lev_stat:.4f}', f'p={lev_p:.3f}',
                   '等分散OK' if lev_p > 0.05 else '不等分散'])
table_data.append([f'Kruskal-Wallis\n（3群比較）', f'H={kw_stat:.4f}', f'p={kw_p:.4f}',
                   '有意**' if kw_p < 0.01 else '有意*' if kw_p < 0.05 else '非有意'])

tbl = ax3b.table(cellText=table_data[1:], colLabels=table_data[0],
                 loc='center', cellLoc='center', bbox=[0.02, 0.0, 0.98, 1.0])
tbl.auto_set_font_size(False)
tbl.set_fontsize(9.5)
for (row, col), cell in tbl.get_celld().items():
    if row == 0:
        cell.set_facecolor('#1565C0')
        cell.set_text_props(color='white', fontweight='bold')
    elif '有意' in str(cell.get_text().get_text()):
        cell.set_facecolor('#FFCDD2')
    elif 'OK' in str(cell.get_text().get_text()):
        cell.set_facecolor('#C8E6C9')
ax3b.set_title('前提確認 + Kruskal-Wallis 検定 サマリー',
               fontsize=11, fontweight='bold', pad=10)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2025_H1_fig3_kruskal.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("  → 2025_H1_fig3_kruskal.png 保存完了")

print("図4: 医師数分布と地域特性を作成中...")

fig4, axes4 = plt.subplots(1, 2, figsize=(14, 6))
fig4.suptitle('医師数_10万対の地域分布（東京都除く46都道府県）',
              fontsize=13, fontweight='bold')

# (a) 都道府県ランキング棒グラフ
ax4a = axes4[0]
nat_mean = df['医師数_10万対'].mean()
df_sorted = df.sort_values('医師数_10万対')
bar_col4 = ['#43A047' if v >= nat_mean else '#E53935' for v in df_sorted['医師数_10万対']]
ax4a.barh(np.arange(len(df_sorted)), df_sorted['医師数_10万対'] - nat_mean,
          color=bar_col4, alpha=0.8, edgecolor='white', left=nat_mean)
ax4a.axvline(nat_mean, color='black', linewidth=1.5,
             label=f'46都道府県平均（={nat_mean:.1f}）')
ax4a.set_yticks(np.arange(len(df_sorted)))
ax4a.set_yticklabels(
    [p.replace('県', '').replace('府', '').replace('道', '').replace('都', '')
     for p in df_sorted['pref']],
    fontsize=7.5
)
ax4a.set_xlabel('医師数_10万対（人/10万人）', fontsize=10)
ax4a.set_title('都道府県別 医師数_10万対\n（緑=平均以上、赤=平均以下）',
               fontsize=10, fontweight='bold')
ax4a.legend(fontsize=9)
ax4a.grid(axis='x', alpha=0.3)

# (b) 医師数_10万対 vs 面積100km²当たり病院診療所数 散布図
ax4b = axes4[1]
sc2 = ax4b.scatter(df['面積100km2当たり病院診療所数'], df['医師数_10万対'],
                   c=df['1人当たり県民所得'], cmap='plasma',
                   s=80, alpha=0.85, edgecolors='white', linewidth=0.5, zorder=3)
plt.colorbar(sc2, ax=ax4b, label='1人当たり県民所得（万円）')

x4 = df['面積100km2当たり病院診療所数'].dropna().values
y4 = df['医師数_10万対'].loc[df['面積100km2当たり病院診療所数'].notna()].values
slope4, intercept4, r4, p4, _ = stats.linregress(x4, y4)
x4_line = np.linspace(x4.min(), x4.max(), 100)
ax4b.plot(x4_line, intercept4 + slope4 * x4_line, '--', color='#333', linewidth=1.8,
          label=f'r={r4:.3f} (p={p4:.3f})')

ax4b.axhline(nat_mean, color='gray', linestyle=':', linewidth=1.0, alpha=0.7)
ax4b.set_xlabel('面積100km²当たり病院診療所数', fontsize=11)
ax4b.set_ylabel('医師数_10万対', fontsize=11)
ax4b.set_title('施設密度 vs 医師数_10万対\n色：1人当たり県民所得',
               fontsize=10, fontweight='bold')
ax4b.legend(fontsize=9)
ax4b.grid(True, alpha=0.3)

# 代表都道府県にラベル
for _, row4 in df.iterrows():
    if row4['pref'] in ['北海道', '沖縄県', '京都府', '秋田県', '鹿児島県', '東京都']:
        ax4b.annotate(row4['pref'].replace('県', '').replace('道', '').replace('府', ''),
                     (row4['面積100km2当たり病院診療所数'], row4['医師数_10万対']),
                     fontsize=8.5, fontweight='bold', color='#333',
                     xytext=(5, 4), textcoords='offset points')

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2025_H1_fig4_specialty.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("  → 2025_H1_fig4_specialty.png 保存完了")

print("\n" + "=" * 65)
print("全図の生成完了（4枚）")
print("=" * 65)
print(f"\n保存先: {FIG_DIR}")
print("  2025_H1_fig1_corr.png      - 相関分析（棒グラフ＋散布図）")
print("  2025_H1_fig2_regression.png - 重回帰係数とVIF")
print("  2025_H1_fig3_kruskal.png   - Kruskal-Wallis検定と前提確認")
print("  2025_H1_fig4_specialty.png  - 地域別医師数分布と施設密度")

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from scipy import stats
from itertools import combinations
import warnings
warnings.filterwarnings('ignore')

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

import os
DATA_DIR = 'data/raw'
FIG_DIR  = 'html/figures'
os.makedirs(FIG_DIR, exist_ok=True)

print("\n" + "=" * 65)
print("■ Step1. 相関分析（医師数_10万対 × 各説明変数）")
print("=" * 65)

CORR_VARS = ['高齢化率', '有効求人倍率', '1人当たり県民所得',
             '年平均気温', '面積100km2当たり病院診療所数', '大学数', '消費支出']
CORR_LABELS = {
    '高齢化率':             '高齢化率（%）',
    '有効求人倍率':         '有効求人倍率（代理）',
    '1人当たり県民所得':    '1人当たり県民所得（万円）',
    '年平均気温':           '年平均気温（℃）',
    '面積100km2当たり病院診療所数': '面積100km²当たり病院診療所数',
    '大学数':               '大学数',
    '消費支出':             '消費支出（二人以上世帯）',
}

print(f"\n  {'変数':<32} {'相関係数':>8} {'p値':>10} {'有意'}")
print("  " + "-" * 58)
corr_results = {}
for var in CORR_VARS:
    mask = df['医師数_10万対'].notna() & df[var].notna()
    x = df.loc[mask, '医師数_10万対'].values
    y = df.loc[mask, var].values
    r, p = stats.pearsonr(x, y)
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    label = CORR_LABELS.get(var, var)
    print(f"  {label:<32} {r:>8.4f} {p:>10.4f} {sig}")
    corr_results[var] = (r, p)

print("\n" + "=" * 65)
print("■ Step2. 重回帰分析（医師数_10万対 ~ 社会・経済・気候・設備要因）")
print("=" * 65)

REG_VARS = ['高齢化率', '有効求人倍率', '1人当たり県民所得',
            '年平均気温', '面積100km2当たり病院診療所数', '大学数']

df_reg = df[['医師数_10万対'] + REG_VARS].dropna().copy()
for v in REG_VARS:
    mu, sg = df_reg[v].mean(), df_reg[v].std()
    df_reg[v + '_z'] = (df_reg[v] - mu) / sg if sg > 0 else 0.0
Z_REG = [v + '_z' for v in REG_VARS]

X_reg = sm.add_constant(df_reg[Z_REG])
y_reg = df_reg['医師数_10万対']
reg_res = sm.OLS(y_reg, X_reg).fit(cov_type='HC1')

print(f"\n重回帰分析結果:")
print(f"  R² = {reg_res.rsquared:.4f}, 調整済みR² = {reg_res.rsquared_adj:.4f}")
print(f"  F統計量 = {reg_res.fvalue:.3f}, p = {reg_res.f_pvalue:.4f}")
print(f"\n  {'変数':<32} {'係数':>8} {'SE':>8} {'p値':>10} {'有意'}")
print("  " + "-" * 62)
for zv in Z_REG:
    b  = reg_res.params.get(zv, np.nan)
    se = reg_res.bse.get(zv, np.nan)
    p  = reg_res.pvalues.get(zv, 1.0)
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '†' if p < 0.1 else ''
    label = CORR_LABELS.get(zv.replace('_z', ''), zv)
    print(f"  {label:<32} {b:>8.4f} {se:>8.4f} {p:>10.4f} {sig}")

# VIF
print(f"\nVIF（多重共線性確認）:")
X_vif = df_reg[Z_REG].copy()
vif_vals = {}
for i, col in enumerate(Z_REG):
    vif = variance_inflation_factor(X_vif.values, i)
    vif_vals[col] = vif
    flag = ' ★多重共線性の疑い' if vif > 5 else ''
    label = CORR_LABELS.get(col.replace('_z', ''), col)
    print(f"  {label:<32} VIF={vif:.2f}{flag}")

print("\n" + "=" * 65)
print("■ Step4. Kruskal-Wallis検定 + Dunn検定（Bonferroni補正）")
print("=" * 65)

kw_stat, kw_p = stats.kruskal(*groups.values())
print(f"\nKruskal-Wallis 検定: H={kw_stat:.4f}, p={kw_p:.6f}")
if kw_p < 0.05:
    print("  → 有意差あり → Dunn 検定（事後検定）を実施")
else:
    print("  → 有意差なし")

# Dunn 検定（Bonferroni補正: Mann-Whitney U で代替）
group_names = list(groups.keys())
pairs = list(combinations(range(3), 2))
n_pairs = len(pairs)
print(f"\nDunn 検定（Bonferroni補正, n_pairs={n_pairs}）:")
print(f"  {'比較':<44} {'U統計量':>10} {'p値（補正前）':>14} {'p値（補正後）':>14} {'有意'}")
print("  " + "-" * 86)
dunn_results = []
for (i, j) in pairs:
    g1, g2 = groups[group_names[i]], groups[group_names[j]]
    u_stat, p_raw = stats.mannwhitneyu(g1, g2, alternative='two-sided')
    p_bonf = min(p_raw * n_pairs, 1.0)
    sig = '***' if p_bonf < 0.001 else '**' if p_bonf < 0.01 else '*' if p_bonf < 0.05 else 'n.s.'
    comp = f"{group_names[i]} vs {group_names[j]}"
    print(f"  {comp:<44} {u_stat:>10.1f} {p_raw:>14.4f} {p_bonf:>14.4f} {sig}")
    dunn_results.append((group_names[i], group_names[j], u_stat, p_raw, p_bonf))

print(f"\n群別 粗死亡率の記述統計:")
for gname, gvals in groups.items():
    print(f"  {gname}: n={len(gvals)}, 中央値={np.median(gvals):.3f}, 平均={gvals.mean():.3f}")

print("\n" + "=" * 65)
print("■ 図の生成（4枚）")
print("=" * 65)