医療資源の地域格差
医療機関・医師数の分布と決定要因分析

import os
import warnings
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from scipy import stats

warnings.filterwarnings('ignore')

# ── パス設定 ──────────────────────────────────────────────────────────────────
FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

plt.rcParams.update({
    'font.family': 'Hiragino Sans',
    'axes.unicode_minus': False,
    'figure.dpi': 150,
    'axes.spines.top': False,
    'axes.spines.right': False,
})

# ── データ読み込み ────────────────────────────────────────────────────────────
raw = pd.read_csv(DATA_B, encoding='cp932', header=0)
# 行0はJapanese label row; 実データは行1以降
df_all = raw.iloc[1:].copy()
df_all.columns = raw.columns

YEAR_COL = df_all.columns[0]  # 'SSDSE-B-2026'

# 都道府県コードのみ（R01000〜R47000）
df_all = df_all[df_all['Code'].str.match(r'^R\d{5}$')].copy()

# 数値変換
NUM_COLS = ['I510120', 'I5102', 'I5103', 'A1101', 'A1303', 'L3221', 'L322106']
for col in NUM_COLS:
    df_all[col] = pd.to_numeric(df_all[col], errors='coerce')

df_all['year'] = df_all[YEAR_COL].astype(int)

# ── 派生指標の計算 ────────────────────────────────────────────────────────────
df_all['病院密度']          = df_all['I510120'] / df_all['A1101'] * 10000
df_all['診療所密度']        = df_all['I5102']   / df_all['A1101'] * 10000
df_all['歯科密度']          = df_all['I5103']   / df_all['A1101'] * 10000
df_all['医療機関密度_総合'] = (df_all['I510120'] + df_all['I5102']) / df_all['A1101'] * 10000
df_all['高齢化率']          = df_all['A1303']   / df_all['A1101'] * 100
df_all['消費支出_log']      = np.log(df_all['L3221'].clip(lower=1))
df_all['保健医療費_千円']   = df_all['L322106'] / 1000

print(f"読み込み完了: {len(df_all)}レコード, {df_all['year'].nunique()}年度, {df_all['Prefecture'].nunique()}都道府県")

# 2022年断面データ
d22 = df_all[df_all['year'] == 2022].copy()
d22 = d22.dropna(subset=['病院密度', '診療所密度', '歯科密度', '高齢化率']).sort_values('医療機関密度_総合', ascending=False)

def gini_coefficient(values):
    """ローレンツ曲線に基づくジニ係数を計算する。"""
    v = np.sort(np.array(values, dtype=float))
    n = len(v)
    idx = np.arange(1, n + 1)
    return (2 * np.dot(idx, v)) / (n * v.sum()) - (n + 1) / n

g_hosp  = gini_coefficient(d22['病院密度'].dropna())
g_clin  = gini_coefficient(d22['診療所密度'].dropna())
g_dent  = gini_coefficient(d22['歯科密度'].dropna())
g_total = gini_coefficient(d22['医療機関密度_総合'].dropna())
print(f"\nジニ係数 (2022年, 47都道府県):")
print(f"  病院密度:         {g_hosp:.4f}")
print(f"  診療所密度:       {g_clin:.4f}")
print(f"  歯科密度:         {g_dent:.4f}")
print(f"  医療機関密度_総合: {g_total:.4f}")

fig1, ax1 = plt.subplots(figsize=(14, 8))

prefs    = d22['Prefecture'].values
hosp_v   = d22['病院密度'].values
clin_v   = d22['診療所密度'].values
dent_v   = d22['歯科密度'].values
x        = np.arange(len(prefs))
bar_w    = 0.72

b1 = ax1.bar(x, hosp_v, bar_w, label='一般病院', color='#1565C0', alpha=0.9)
b2 = ax1.bar(x, clin_v, bar_w, bottom=hosp_v, label='一般診療所', color='#42A5F5', alpha=0.9)
b3 = ax1.bar(x, dent_v, bar_w, bottom=hosp_v + clin_v, label='歯科診療所', color='#E65100', alpha=0.85)

ax1.set_xticks(x)
ax1.set_xticklabels(prefs, rotation=90, fontsize=8.5)
ax1.set_ylabel('人口1万人当たり医療機関数（機関）', fontsize=11)
ax1.set_title('Fig 1｜都道府県別 医療機関密度（2022年）\n病院・診療所・歯科の積み上げ',
              fontsize=13, fontweight='bold', pad=14)
ax1.legend(loc='upper right', fontsize=10)
ax1.yaxis.set_major_formatter(mticker.FormatStrFormatter('%.1f'))

# 全国平均線
nat_avg = (d22['医療機関密度_総合'] + d22['歯科密度']).mean()
ax1.axhline(nat_avg, color='gray', linestyle='--', linewidth=1.2, alpha=0.8)
ax1.text(len(prefs) - 0.5, nat_avg + 0.08, f'全国平均\n{nat_avg:.2f}', fontsize=8.5,
         color='gray', ha='right')

fig1.tight_layout()
fig1_path = os.path.join(FIG_DIR, '2021_U4_fig1_ranking.png')
fig1.savefig(fig1_path, dpi=150, bbox_inches='tight')
plt.close(fig1)
print(f"\nFig1 保存: {fig1_path}")

d22_sc = d22.dropna(subset=['高齢化率', '医療機関密度_総合'])

fig2, ax2 = plt.subplots(figsize=(10, 7))

# 地域区分（総務省8地域）
region_map = {
    '北海道': '北海道',
    '青森県': '東北', '岩手県': '東北', '宮城県': '東北', '秋田県': '東北',
    '山形県': '東北', '福島県': '東北',
    '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東',
    '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国', '島根県': '中国', '岡山県': '中国', '広島県': '中国', '山口県': '中国',
    '徳島県': '四国', '香川県': '四国', '愛媛県': '四国', '高知県': '四国',
    '福岡県': '九州', '佐賀県': '九州', '長崎県': '九州', '熊本県': '九州',
    '大分県': '九州', '宮崎県': '九州', '鹿児島県': '九州', '沖縄県': '九州',
}
region_colors = {
    '北海道': '#1565C0', '東北': '#2E7D32', '関東': '#F57F17',
    '中部': '#6A1B9A', '近畿': '#C62828', '中国': '#00838F',
    '四国': '#E65100', '九州': '#37474F',
}

d22_sc = d22_sc.copy()
d22_sc['region'] = d22_sc['Prefecture'].map(region_map).fillna('その他')

plotted_regions = set()
for _, row in d22_sc.iterrows():
    reg = row['region']
    col = region_colors.get(reg, '#888888')
    label = reg if reg not in plotted_regions else '_nolegend_'
    ax2.scatter(row['高齢化率'], row['医療機関密度_総合'],
                color=col, s=60, alpha=0.85, label=label, zorder=3)
    plotted_regions.add(reg)
    ax2.annotate(row['Prefecture'].replace('県', '').replace('都', '').replace('府', '').replace('道', ''),
                 (row['高齢化率'], row['医療機関密度_総合']),
                 textcoords='offset points', xytext=(5, 3), fontsize=7.5, color='#333')

# 回帰直線
x_sc = d22_sc['高齢化率'].values
y_sc = d22_sc['医療機関密度_総合'].values
slope, intercept, r, p, _ = stats.linregress(x_sc, y_sc)
x_line = np.linspace(x_sc.min(), x_sc.max(), 100)
ax2.plot(x_line, slope * x_line + intercept, 'k--', linewidth=1.5, alpha=0.6,
         label=f'回帰直線 (r={r:.3f}, p={p:.3f})')

ax2.set_xlabel('高齢化率（65歳以上人口比率, %）', fontsize=11)
ax2.set_ylabel('医療機関密度_総合（人口1万人当たり）', fontsize=11)
ax2.set_title('Fig 2｜高齢化率 × 医療機関密度（2022年・47都道府県）',
              fontsize=13, fontweight='bold', pad=14)
ax2.legend(fontsize=9, loc='upper left', ncol=2)

fig2.tight_layout()
fig2_path = os.path.join(FIG_DIR, '2021_U4_fig2_scatter.png')
fig2.savefig(fig2_path, dpi=150, bbox_inches='tight')
plt.close(fig2)
print(f"Fig2 保存: {fig2_path}")

d_ols = d22.dropna(subset=['医療機関密度_総合', '高齢化率', '消費支出_log', '保健医療費_千円']).copy()

def standardize(s):
    return (s - s.mean()) / s.std()

y_std = standardize(d_ols['医療機関密度_総合'])
X_dict = {
    '高齢化率': standardize(d_ols['高齢化率']),
    '消費支出(log)': standardize(d_ols['消費支出_log']),
    '保健医療費(千円)': standardize(d_ols['保健医療費_千円']),
}

coef_results = {}
for name, xvar in X_dict.items():
    X_mat = np.column_stack([np.ones(len(xvar)), xvar.values])
    beta, res, rank, sv = np.linalg.lstsq(X_mat, y_std.values, rcond=None)
    # 標準誤差
    n_, k_ = X_mat.shape
    y_hat = X_mat @ beta
    sigma2 = np.sum((y_std.values - y_hat)**2) / (n_ - k_)
    cov_mat = sigma2 * np.linalg.inv(X_mat.T @ X_mat)
    se = np.sqrt(np.diag(cov_mat))
    t_stat = beta[1] / se[1]
    p_val = 2 * stats.t.sf(abs(t_stat), df=n_ - k_)
    ci95 = 1.96 * se[1]
    coef_results[name] = {'coef': beta[1], 'ci': ci95, 'p': p_val}

# 複合OLS（全変数投入）
X_all = np.column_stack([
    np.ones(len(d_ols)),
    X_dict['高齢化率'].values,
    X_dict['消費支出(log)'].values,
    X_dict['保健医療費(千円)'].values,
])
beta_all, _, _, _ = np.linalg.lstsq(X_all, y_std.values, rcond=None)
n_a, k_a = X_all.shape
y_hat_a = X_all @ beta_all
sigma2_a = np.sum((y_std.values - y_hat_a)**2) / (n_a - k_a)
cov_all = sigma2_a * np.linalg.inv(X_all.T @ X_all)
se_all = np.sqrt(np.diag(cov_all))

multi_coefs = {
    '高齢化率': (beta_all[1], 1.96 * se_all[1],
                 2 * stats.t.sf(abs(beta_all[1] / se_all[1]), df=n_a - k_a)),
    '消費支出(log)': (beta_all[2], 1.96 * se_all[2],
                     2 * stats.t.sf(abs(beta_all[2] / se_all[2]), df=n_a - k_a)),
    '保健医療費(千円)': (beta_all[3], 1.96 * se_all[3],
                       2 * stats.t.sf(abs(beta_all[3] / se_all[3]), df=n_a - k_a)),
}

print("\n--- OLS結果 (標準化係数, N=47都道府県, 2022年) ---")
for vname, (c, ci, p) in multi_coefs.items():
    sig = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else 'n.s.'))
    print(f"  {vname:18s}: β={c:+.4f}  95%CI[{c-ci:+.4f}, {c+ci:+.4f}]  p={p:.4f} {sig}")

fig3, ax3 = plt.subplots(figsize=(8, 5))

var_labels = list(multi_coefs.keys())
coefs  = [multi_coefs[v][0] for v in var_labels]
cis    = [multi_coefs[v][1] for v in var_labels]
pvals  = [multi_coefs[v][2] for v in var_labels]
colors_ols = ['#1565C0' if c > 0 else '#C62828' for c in coefs]

y_pos = np.arange(len(var_labels))
ax3.barh(y_pos, coefs, xerr=cis, color=colors_ols, alpha=0.80,
         height=0.5, error_kw={'elinewidth': 2, 'capsize': 5, 'ecolor': '#333'})

for i, (c, p) in enumerate(zip(coefs, pvals)):
    sig_str = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else 'n.s.'))
    ax3.text(c + (cis[i] + 0.02) * np.sign(c), i, sig_str,
             va='center', fontsize=11, fontweight='bold',
             color='#333' if sig_str != 'n.s.' else '#999')

ax3.axvline(0, color='black', linewidth=1.0, linestyle='-')
ax3.set_yticks(y_pos)
ax3.set_yticklabels(var_labels, fontsize=11)
ax3.set_xlabel('標準化回帰係数（β）± 95%信頼区間', fontsize=11)
ax3.set_title('Fig 3｜医療機関密度の決定要因\nOLS標準化係数（2022年断面, N=47）', fontsize=13, fontweight='bold', pad=14)
ax3.text(0.99, 0.03, '***p<0.001  **p<0.01  *p<0.05  n.s.=有意でない',
         transform=ax3.transAxes, ha='right', va='bottom', fontsize=9, color='#555')

fig3.tight_layout()
fig3_path = os.path.join(FIG_DIR, '2021_U4_fig3_ols.png')
fig3.savefig(fig3_path, dpi=150, bbox_inches='tight')
plt.close(fig3)
print(f"Fig3 保存: {fig3_path}")

df_ts = df_all[df_all['year'].between(2012, 2023)].copy()
df_ts['region'] = df_ts['Prefecture'].map(region_map).fillna('その他')
df_ts = df_ts.dropna(subset=['医療機関密度_総合'])

region_ts = (
    df_ts.groupby(['year', 'region'])['医療機関密度_総合']
    .mean()
    .reset_index()
)

region_order = ['北海道', '東北', '関東', '中部', '近畿', '中国', '四国', '九州']
region_ls = ['-', '--', '-.', ':', '-', '--', '-.', ':']

fig4, ax4 = plt.subplots(figsize=(11, 6))

for reg, ls_ in zip(region_order, region_ls):
    sub = region_ts[region_ts['region'] == reg].sort_values('year')
    col = region_colors.get(reg, '#888')
    ax4.plot(sub['year'], sub['医療機関密度_総合'],
             marker='o', markersize=5, linewidth=2.2, linestyle=ls_,
             color=col, label=reg, alpha=0.9)

ax4.set_xlabel('年度', fontsize=11)
ax4.set_ylabel('医療機関密度_総合（人口1万人当たり）', fontsize=11)
ax4.set_title('Fig 4｜地域別 医療機関密度の時系列推移（2012〜2023年）\n8地域平均（一般病院＋一般診療所）',
              fontsize=13, fontweight='bold', pad=14)
ax4.set_xticks(range(2012, 2024))
ax4.xaxis.set_tick_params(rotation=45)
ax4.legend(loc='upper right', fontsize=10, ncol=2)
ax4.yaxis.set_major_formatter(mticker.FormatStrFormatter('%.1f'))

fig4.tight_layout()
fig4_path = os.path.join(FIG_DIR, '2021_U4_fig4_timeseries.png')
fig4.savefig(fig4_path, dpi=150, bbox_inches='tight')
plt.close(fig4)
print(f"Fig4 保存: {fig4_path}")

print("\n=== 全図の出力が完了しました ===")
print(f"出力先: {FIG_DIR}")