潜在患者数に対する医師偏在の可視化
入院患者数あたり医師数のジニ係数による都道府県比較

# ===== ステップ1: SSDSE-A(市区町村)から従来指標を実データ計算 =====
# 1行目=変数コード、2行目=年度、3行目=日本語名 → コード行を列名に使う
df = pd.read_csv(DATA_A, encoding='cp932', header=0, skiprows=[1, 2])
df = df.rename(columns={df.columns[0]: 'code',
                        'Prefecture': '都道府県',
                        'Municipality': '市区町村'})
df['総人口'] = pd.to_numeric(df['A1101'], errors='coerce')   # 2020年
df['医師数'] = pd.to_numeric(df['I6100'], errors='coerce')   # 2022年
print(f"【SSDSE-A-2025】{len(df)}市区町村(医師数=2022年、総人口=2020年)")
print('  ※ 原論文が使った SSDSE-2019A は医師数2016年・総人口2015年。年次が異なる')
print('  ※ 市区町村→二次医療圏の対応表(平成28年医療施設調査)と')
print('     二次医療圏別の推計入院患者数(平成29年患者調査)は SSDSE 未収録')
print('     → 原論文の中核「二次医療圏別ジニ係数」は再計算不可(報告値で提示)')

pref = df.groupby('都道府県', sort=False).agg(
    医師数=('医師数', 'sum'), 総人口=('総人口', 'sum'))
pref['人口10万対医師数'] = pref['医師数'] / pref['総人口'] * 100000
pref_sorted = pref.sort_values('人口10万対医師数', ascending=False)
print('\n=== 従来指標「人口10万人あたり医師数」(SSDSE-A-2025 実データ計算) ===')
print('上位5:')
print(pref_sorted.head(5)[['人口10万対医師数']].round(1).to_string())
print('下位5:')
print(pref_sorted.tail(5)[['人口10万対医師数']].round(1).to_string())
print(f"全国計: 医師数 {int(pref['医師数'].sum()):,}人 / "
      f"人口 {int(pref['総人口'].sum()):,}人 = "
      f"10万対 {pref['医師数'].sum()/pref['総人口'].sum()*1e5:.1f}人")
print('  → 上位と下位で約2倍の開き。原論文はこの従来指標に')
print('    「医療需要(ニーズ)が反映されていない」ことを問題にした')

# ===== ステップ2: 図1 医師数の推移(原論文 図1 の報告値) =====
# 出典: 医師・歯科医師・薬剤師調査(原論文 図1 に記載された値)
years = [2000, 2002, 2004, 2006, 2008, 2010, 2012, 2014, 2016]
docs  = [255792, 262687, 270371, 277927, 286699,
         295049, 303268, 311205, 319480]

fig, ax = plt.subplots(figsize=(8.6, 5.2))
bars = ax.bar([str(y) + '年' for y in years], docs,
              color='#5B9BD5', edgecolor='white', width=0.62)
for b, v in zip(bars, docs):
    ax.text(b.get_x() + b.get_width() / 2, v + 3500, f'{v:,}',
            ha='center', va='bottom', fontsize=8.2)
ax.set_ylim(0, 360000)
ax.yaxis.set_major_formatter(lambda x, _: f'{int(x):,}')
ax.set_ylabel('医師数(人)', fontsize=11)
ax.set_title('図1:医師数の推移(2000〜2016年)\n'
             '原論文 図1 の報告値の可視化(出典:医師・歯科医師・薬剤師調査。再計算ではない)',
             fontsize=11.5, fontweight='bold')
ax.grid(axis='y', alpha=0.3)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2019_U5_1_fig1.png'),
            dpi=150, bbox_inches='tight')
plt.close(fig)
print('\n[保存] 図1: 医師数の推移(原論文 図1 の報告値の可視化)')
print(f'  2000年 {docs[0]:,}人 → 2016年 {docs[-1]:,}人(+{docs[-1]-docs[0]:,}人)')
print('  → 医師総数は一貫して増加。それでも「偏在」は解消していない、が出発点')

# ===== ステップ3: 原論文 表2(都道府県別ジニ係数と差分)の転記 =====
# 値はすべて原論文 表2 の報告値(入院患者数あたり, 人口あたり, 差分)
GINI = {
    '北海道': (0.195, 0.173, 0.022), '青森県': (0.196, 0.175, 0.021),
    '岩手県': (0.218, 0.181, 0.037), '宮城県': (0.188, 0.117, 0.071),
    '秋田県': (0.195, 0.181, 0.014), '山形県': (0.107, 0.125, -0.018),
    '福島県': (0.177, 0.129, 0.049), '茨城県': (0.290, 0.257, 0.032),
    '栃木県': (0.293, 0.229, 0.064), '群馬県': (0.250, 0.218, 0.032),
    '埼玉県': (0.113, 0.109, 0.004), '千葉県': (0.171, 0.124, 0.047),
    '東京都': (0.359, 0.326, 0.033), '神奈川県': (0.106, 0.112, -0.006),
    '新潟県': (0.141, 0.138, 0.003), '富山県': (0.117, 0.104, 0.013),
    '石川県': (0.175, 0.146, 0.030), '福井県': (0.234, 0.217, 0.018),
    '山梨県': (0.175, 0.154, 0.021), '長野県': (0.179, 0.169, 0.010),
    '岐阜県': (0.151, 0.134, 0.017), '静岡県': (0.131, 0.113, 0.018),
    '愛知県': (0.176, 0.200, -0.025), '三重県': (0.123, 0.107, 0.016),
    '滋賀県': (0.215, 0.191, 0.024), '京都府': (0.167, 0.170, -0.003),
    '大阪府': (0.146, 0.126, 0.020), '兵庫県': (0.165, 0.123, 0.043),
    '奈良県': (0.109, 0.085, 0.024), '和歌山県': (0.205, 0.185, 0.020),
    '鳥取県': (0.197, 0.160, 0.037), '島根県': (0.292, 0.205, 0.088),
    '岡山県': (0.150, 0.105, 0.046), '広島県': (0.152, 0.083, 0.069),
    '山口県': (0.161, 0.138, 0.023), '徳島県': (0.154, 0.083, 0.071),
    '香川県': (0.211, 0.136, 0.074), '愛媛県': (0.244, 0.147, 0.097),
    '高知県': (0.119, 0.111, 0.008), '福岡県': (0.224, 0.160, 0.064),
    '佐賀県': (0.210, 0.171, 0.039), '長崎県': (0.150, 0.167, -0.017),
    '熊本県': (0.302, 0.220, 0.082), '大分県': (0.196, 0.104, 0.091),
    '宮崎県': (0.293, 0.192, 0.100), '鹿児島県': (0.276, 0.205, 0.071),
    '沖縄県': (0.159, 0.126, 0.033),
}
gini = pd.DataFrame(GINI, index=['G入院患者数', 'G人口', '差分']).T

print('\n=== 原論文 表2 の報告値(都道府県別ジニ係数)の検算 ===')
for col, label, rep in [('G人口', '人口あたり', (0.326, 0.083, 0.147)),
                        ('G入院患者数', '入院患者数あたり', (0.359, 0.106, 0.177))]:
    mx, mn, md = gini[col].max(), gini[col].min(), gini[col].median()
    ok = '一致' if (round(mx, 3), round(mn, 3), round(md, 3)) == rep else '不一致'
    print(f'  {label}: 最大 {mx:.3f} / 最小 {mn:.3f} / 中央値 {md:.3f} '
          f'→ 原論文本文の報告値 {rep} と{ok}')
print('  ※ 表2の値から取った要約であり、ジニ係数そのものの再計算ではない')

# ===== ステップ5: 図2・図3 ジニ係数の都道府県比較(報告値の再表現) =====
norm_g = Normalize(vmin=0.0, vmax=0.4)   # 原論文の凡例 0.000〜0.400 に合わせる
tile_map(gini['G人口'],
         '図2:人口あたりの医師数でみたジニ係数の都道府県比較\n'
         '原論文 表2・図2 の報告値の再表現(タイルマップ。再計算ではない)',
         'Blues', norm_g, 'ジニ係数(0=平等、0.4=不平等の目安)',
         '2019_U5_1_fig2.png')
print('\n[保存] 図2: 人口あたり医師数のジニ係数(原論文 表2 の報告値の再表現)')
print('  最大: 東京都 0.326 / 最小: 広島県・徳島県 0.083 / 中央値 0.147')

tile_map(gini['G入院患者数'],
         '図3:入院患者数あたり医師数でみたジニ係数の都道府県比較\n'
         '原論文 表2・図3 の報告値の再表現(タイルマップ。再計算ではない)',
         'Blues', norm_g, 'ジニ係数(0=平等、0.4=不平等の目安)',
         '2019_U5_1_fig3.png')
print('[保存] 図3: 入院患者数あたり医師数のジニ係数(同上)')
print('  最大: 東京都 0.359 / 最小: 神奈川県 0.106 / 中央値 0.177')
print('  → 最大・最小・中央値のすべてが人口あたり(図2)より高い')

# ===== ステップ6: 図4 2つのジニ係数の差分(正=青・負=赤) =====
norm_d = TwoSlopeNorm(vmin=-0.10, vcenter=0.0, vmax=0.10)
tile_map(gini['差分'],
         '図4:2つのジニ係数の差分の都道府県比較(正=青・負=赤)\n'
         '差分 = Gini(入院患者数) − Gini(人口)。'
         '原論文 表2・図4 の報告値の再表現(再計算ではない)',
         'RdBu', norm_d, '差分(青: 潜在的医療需要が高い / 赤: 低い)',
         '2019_U5_1_fig4.png', fmt='{:+.3f}')
print('[保存] 図4: 2つのジニ係数の差分(同上)')

pos = gini[gini['差分'] > 0]
neg = gini[gini['差分'] < 0]
west = ['滋賀県', '京都府', '大阪府', '兵庫県', '奈良県', '和歌山県', '三重県',
        '鳥取県', '島根県', '岡山県', '広島県', '山口県', '徳島県', '香川県',
        '愛媛県', '高知県', '福岡県', '佐賀県', '長崎県', '熊本県', '大分県',
        '宮崎県', '鹿児島県', '沖縄県']
top5 = gini['差分'].sort_values(ascending=False).head(5)
print('\n=== 差分(原論文 表2 の報告値)の内訳 ===')
print(f'  差分が正(青): {len(pos)}都道府県 / 負(赤): {len(neg)}県'
      f'({"、".join(neg.index)})')
print('  差分の大きい上位5県: '
      + '、'.join(f'{p} {v:+.3f}' for p, v in top5.items()))
print(f'  上位5県のうち西日本(近畿以西): {sum(p in west for p in top5.index)}県')
print('  → 「潜在的医療需要に対する医師偏在はおもに西日本で発生」という')
print('    原論文の指摘(図4)と整合する')

print('\n====== まとめ ======')
print('・従来指標(人口10万対医師数)は SSDSE-A-2025 実データで再計算し、')
print('  都道府県間で約2倍の開きを確認(実データ計算)')
print('・二次医療圏別ジニ係数(図2〜4)は原論文 表2 の報告値の可視化であり、')
print('  再計算ではない(対応表・患者調査データが SSDSE 未収録のため)')
print('・入院患者数あたりのジニ係数(最大0.359・中央値0.177)は人口あたり')
print('  (最大0.326・中央値0.147)より一様に高く、差分の大きい県は西日本に集中')
print('\nDONE: 2019_U5_1_shorei')

# ===== ステップ4: タイルマップ(コロプレス図の再表現)関数 =====
# 都道府県を格子状に並べた「タイルマップ」。原論文の地図(コロプレス図)の
# 再表現(白地図データを使わずに塗り分けを表現する方法)
TILE = {  # (列, 行) 行は上から下へ
    '北海道': (13, 0), '青森県': (13, 1), '秋田県': (12, 2), '岩手県': (13, 2),
    '山形県': (12, 3), '宮城県': (13, 3), '新潟県': (10, 4), '福島県': (12, 4),
    '富山県': (9, 4), '石川県': (8, 4), '福井県': (8, 5), '長野県': (10, 5),
    '群馬県': (11, 5), '栃木県': (12, 5), '茨城県': (13, 5), '岐阜県': (9, 5),
    '山梨県': (11, 6), '埼玉県': (12, 6), '千葉県': (13, 6), '愛知県': (9, 6),
    '静岡県': (10, 6), '東京都': (12, 7), '神奈川県': (11, 7), '滋賀県': (8, 6),
    '京都府': (7, 5), '兵庫県': (6, 5), '鳥取県': (5, 5), '島根県': (4, 5),
    '大阪府': (7, 6), '岡山県': (5, 6), '広島県': (4, 6), '山口県': (3, 6),
    '奈良県': (8, 7), '和歌山県': (7, 7), '三重県': (9, 7),
    '香川県': (5, 7), '徳島県': (6, 7), '愛媛県': (4, 7), '高知県': (5, 8),
    '福岡県': (2, 6), '佐賀県': (1, 6), '長崎県': (0, 7), '大分県': (2, 7),
    '熊本県': (1, 7), '宮崎県': (2, 8), '鹿児島県': (1, 8), '沖縄県': (0, 10),
}


def tile_map(values, title, cmap, norm, cbar_label, fname, fmt='{:.3f}'):
    fig, ax = plt.subplots(figsize=(10.4, 8.2))
    sm = ScalarMappable(norm=norm, cmap=cmap)
    for p, (cx, cy) in TILE.items():
        v = values[p]
        fc = sm.to_rgba(v)
        lum = 0.299 * fc[0] + 0.587 * fc[1] + 0.114 * fc[2]
        tc = 'white' if lum < 0.55 else '#222222'
        ax.add_patch(Rectangle((cx, -cy), 0.96, 0.96, facecolor=fc,
                               edgecolor='white', linewidth=1.2))
        name = p if p == '北海道' else p[:-1]  # 「県」「府」「都」を省く
        ax.text(cx + 0.48, -cy + 0.60, name, ha='center', va='center',
                fontsize=7.8, color=tc, fontweight='bold')
        ax.text(cx + 0.48, -cy + 0.28, fmt.format(v), ha='center',
                va='center', fontsize=7.2, color=tc)
    ax.set_xlim(-0.4, 15.4)
    ax.set_ylim(-10.6, 1.6)
    ax.set_aspect('equal')
    ax.axis('off')
    cb = fig.colorbar(sm, ax=ax, shrink=0.62, pad=0.01)
    cb.set_label(cbar_label, fontsize=10)
    ax.set_title(title, fontsize=11.5, fontweight='bold')
    fig.tight_layout()
    fig.savefig(os.path.join(FIG_DIR, fname), dpi=150, bbox_inches='tight')
    plt.close(fig)