日本におけるMRI設置の現状と
過剰導入の実証的検証

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
from scipy.optimize import minimize
from scipy.special import ndtr   # Φ(x)

warnings.filterwarnings('ignore')

# ── パス設定 ──────────────────────────────────────────────────────────
DATA_DIR = 'data/2025_U4'
FIG_DIR  = 'html/figures'
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,
})

# ── データ読み込み ─────────────────────────────────────────────────────
df = pd.read_csv(os.path.join(DATA_DIR, '2025_U4_panel.csv'))
df = df.dropna(subset=['population', 'mri_total'])
df['ln_pop'] = np.log(df['population'].clip(lower=1))
df['ln_pop65'] = np.log((df['population'] * 0.28).clip(lower=1))  # 高齢化率約28%を推計に使用

print(f"データ: {len(df)}レコード, {df['year'].nunique()}年度")
print(f"MRI台数統計:\n{df['mri_total'].describe()}\n")

# ── Figure 1: MRI台数分布（年度別） ──────────────────────────────────
fig, axes = plt.subplots(1, 3, figsize=(14, 4), sharey=True)
YEARS = [2017, 2018, 2020]
colors = ['#2196F3', '#4CAF50', '#FF9800']
bins_label = ['0', '1-2', '3-5', '6-10', '11+']

for ax, year, col in zip(axes, YEARS, colors):
    d = df[df['year'] == year]
    counts = [
        (d['mri_total'] == 0).sum(),
        ((d['mri_total'] >= 1) & (d['mri_total'] <= 2)).sum(),
        ((d['mri_total'] >= 3) & (d['mri_total'] <= 5)).sum(),
        ((d['mri_total'] >= 6) & (d['mri_total'] <= 10)).sum(),
        (d['mri_total'] >= 11).sum(),
    ]
    bars = ax.bar(bins_label, counts, color=col, alpha=0.85, edgecolor='white', linewidth=0.8)
    for bar, cnt in zip(bars, counts):
        ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1,
                str(cnt), ha='center', va='bottom', fontsize=9)
    ax.set_title(f'{year}年度', fontsize=12, fontweight='bold')
    ax.set_xlabel('MRI台数（二次医療圏）', fontsize=10)
    ax.set_ylim(0, max(counts) * 1.18)

axes[0].set_ylabel('二次医療圏数', fontsize=10)
fig.suptitle('二次医療圏別MRI台数分布（2017〜2020年度）', fontsize=13, fontweight='bold', y=1.01)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2025_U4_dist.png'), bbox_inches='tight', dpi=150)
plt.close(fig)

print("Figure 1 保存: 2025_U4_dist.png")

# ── Figure 2: 人口 vs MRI台数 散布図（2020年） ───────────────────────
d20 = df[df['year'] == 2020].copy()

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 左: 散布図（生値）
ax = axes[0]
sc = ax.scatter(d20['population'] / 1e4, d20['mri_total'],
                alpha=0.5, s=30, c=d20['mri_total'], cmap='YlOrRd',
                vmin=0, vmax=40, edgecolors='none')
plt.colorbar(sc, ax=ax, label='MRI台数')
ax.set_xlabel('人口（万人）', fontsize=11)
ax.set_ylabel('MRI台数', fontsize=11)
ax.set_title('人口規模とMRI台数（2020年）', fontsize=12, fontweight='bold')
ax.set_xlim(0, d20['population'].max() / 1e4 * 1.05)

# 右: 対数スケール
ax = axes[1]
ax.scatter(d20['ln_pop'], d20['mri_total'],
           alpha=0.5, s=30, c='#1976D2', edgecolors='none')
z = np.polyfit(d20['ln_pop'].dropna(), d20.loc[d20['ln_pop'].notna(), 'mri_total'], 1)
p = np.poly1d(z)
xr = np.linspace(d20['ln_pop'].min(), d20['ln_pop'].max(), 100)
ax.plot(xr, p(xr), 'r-', linewidth=1.5, label=f'線形近似 (β={z[0]:.2f})')
ax.set_xlabel('ln(人口)', fontsize=11)
ax.set_ylabel('MRI台数', fontsize=11)
ax.set_title('対数人口とMRI台数（2020年）', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)

fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2025_U4_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig)

print("Figure 2 保存: 2025_U4_scatter.png")

# ── 順序プロビット推定（Ordered Probit）──────────────────────────────
print("\n" + "=" * 55)
print("■ 順序プロビット推定（Bresnahan-Reiss型）")
print("=" * 55)

# 2020年データで推定
d = d20.copy()
d = d.dropna(subset=['mri_bin', 'ln_pop'])
y = d['mri_bin'].values.astype(int)   # 0,1,2,3,4
X = d['ln_pop'].values

# K=5カテゴリのordered probit
# P(Y=k) = Φ(μ_k - β*X) - Φ(μ_{k-1} - β*X)
# パラメータ: β, μ_1, δ_2=μ_2-μ_1>0, δ_3>0, δ_4>0

def op_loglik(params, y, X):
    beta = params[0]
    mu1  = params[1]
    d2   = np.exp(params[2])
    d3   = np.exp(params[3])
    d4   = np.exp(params[4])
    mu2  = mu1 + d2
    mu3  = mu2 + d3
    mu4  = mu3 + d4
    thresholds = [-np.inf, mu1, mu2, mu3, mu4, np.inf]
    lp = beta * X
    ll = 0.0
    for k in range(5):
        idx = y == k
        if idx.sum() == 0:
            continue
        p = ndtr(thresholds[k+1] - lp[idx]) - ndtr(thresholds[k] - lp[idx])
        p = np.clip(p, 1e-15, 1.0)
        ll += np.log(p).sum()
    return -ll

# 初期値
p0 = [1.0, 5.0, np.log(1.5), np.log(1.5), np.log(1.5)]
result = minimize(op_loglik, p0, args=(y, X),
                  method='Nelder-Mead',
                  options={'maxiter': 20000, 'xatol': 1e-8, 'fatol': 1e-8})

beta_hat = result.x[0]
mu1_hat  = result.x[1]
mu2_hat  = mu1_hat + np.exp(result.x[2])
mu3_hat  = mu2_hat + np.exp(result.x[3])
mu4_hat  = mu3_hat + np.exp(result.x[4])

print(f"\n  係数 β (ln_pop): {beta_hat:.4f}")
print(f"  閾値 μ₁: {mu1_hat:.4f}")
print(f"  閾値 μ₂: {mu2_hat:.4f}")
print(f"  閾値 μ₃: {mu3_hat:.4f}")
print(f"  閾値 μ₄: {mu4_hat:.4f}")
print(f"  対数尤度: {-result.fun:.2f}")

# 必要人口閾値（βが正なら μ_k/β で求まる）
# P(Y>=k) >= 0.5 となる X = μ_k / β
if beta_hat > 0:
    thresholds_pop = {
        1: np.exp(mu1_hat / beta_hat),
        2: np.exp(mu2_hat / beta_hat),
        3: np.exp(mu3_hat / beta_hat),
        4: np.exp(mu4_hat / beta_hat),
    }
    print(f"\n  【必要人口閾値（P(MRI≥N)=0.5）】")
    for n, pop in thresholds_pop.items():
        print(f"  MRI {n}台以上: {pop/1e4:.1f}万人")

    # 競争効果: 1台目→2台目の必要人口比 (BR比)
    # 1台あたりの必要人口増加率が1を超えると競争で利潤が圧縮されている
    ratios = {}
    pops = list(thresholds_pop.values())
    for i in range(1, 4):
        ratios[i+1] = pops[i] / pops[i-1]
    print(f"\n  【競争効果（BR比）: N台→N+1台の必要人口比】")
    for n, r in ratios.items():
        print(f"  {n-1}台→{n}台: {r:.3f} (>1.0 → 競争で利潤圧縮)")

# ── Figure 3: 必要人口閾値と予測確率 ─────────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# 左: 予測確率曲線
ax = axes[0]
ln_range = np.linspace(8, 14, 200)
pop_range = np.exp(ln_range) / 1e4
mus = [mu1_hat, mu2_hat, mu3_hat, mu4_hat]
mu_ext = [-np.inf] + mus + [np.inf]
labels_p = ['0台', '1-2台', '3-5台', '6-10台', '11台+']
cm = plt.cm.get_cmap('tab10')
for k in range(5):
    p_k = ndtr(mu_ext[k+1] - beta_hat * ln_range) - ndtr(mu_ext[k] - beta_hat * ln_range)
    ax.plot(pop_range, p_k, label=labels_p[k], linewidth=2, color=cm(k))
ax.set_xlabel('人口（万人）', fontsize=11)
ax.set_ylabel('確率', fontsize=11)
ax.set_title('MRI台数カテゴリの予測確率\n（順序プロビット）', fontsize=11, fontweight='bold')
ax.legend(fontsize=9, loc='center right')
ax.set_xlim(0, 250)
ax.set_ylim(0, 1)

# 右: 必要人口閾値の棒グラフ（競争効果可視化）
ax = axes[1]
if beta_hat > 0:
    ns = list(thresholds_pop.keys())
    pops_wan = [v / 1e4 for v in thresholds_pop.values()]
    bars = ax.bar([f'{n}台以上' for n in ns], pops_wan,
                  color=['#42A5F5', '#26C6DA', '#66BB6A', '#FFA726'],
                  edgecolor='white', linewidth=0.8, alpha=0.9)
    for bar, v in zip(bars, pops_wan):
        ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.5,
                f'{v:.1f}万人', ha='center', va='bottom', fontsize=9)
    ax.set_ylabel('必要人口（万人）', fontsize=11)
    ax.set_title('MRI設置に必要な人口閾値\n（Bresnahan-Reiss推定値）',
                 fontsize=11, fontweight='bold')

fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2025_U4_br.png'), bbox_inches='tight', dpi=150)
plt.close(fig)

print("\nFigure 3 保存: 2025_U4_br.png")

# ── Figure 4: 地域別過剰設置マップ（散布図で代替）──────────────────
fig, ax = plt.subplots(figsize=(10, 6))

d20 = d20.copy()
d20['expected_bin'] = (beta_hat * d20['ln_pop'] - mu1_hat) / max(abs(beta_hat), 0.01)
d20['excess'] = d20['mri_total'] - d20['mri_bin']   # 実際 - ビン中央

# 人口 vs mri_per_100k（人口100万人あたりMRI台数）で過剰を可視化
sc = ax.scatter(d20['population'] / 1e4, d20['mri_per_100k'],
                c=d20['mri_per_100k'], cmap='RdYlGn_r',
                s=50, alpha=0.7, edgecolors='none',
                vmin=0, vmax=d20['mri_per_100k'].quantile(0.95))

plt.colorbar(sc, ax=ax, label='MRI台数/人口10万人')

# 日本平均線
avg_per100k = df[df['year'] == 2020]['mri_per_100k'].mean()
ax.axhline(avg_per100k, color='red', linewidth=1.5, linestyle='--',
           label=f'全国平均 {avg_per100k:.1f}台/10万人')

ax.set_xlabel('人口（万人）', fontsize=11)
ax.set_ylabel('MRI台数（人口10万人あたり）', fontsize=11)
ax.set_title('人口規模と人口当たりMRI密度（2020年）\n色が赤いほどMRI密度が高い',
             fontsize=11, fontweight='bold')
ax.legend(fontsize=10)
ax.set_xlim(0, d20['population'].max() / 1e4 * 1.05)

# 上位二次医療圏にラベル
top = d20.nlargest(5, 'mri_per_100k')
for _, row in top.iterrows():
    ax.annotate(row['sec_name'], (row['population']/1e4, row['mri_per_100k']),
                fontsize=7, alpha=0.8,
                xytext=(5, 2), textcoords='offset points')

fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2025_U4_excess.png'), bbox_inches='tight', dpi=150)
plt.close(fig)

print("Figure 4 保存: 2025_U4_excess.png")

# ── Figure 5: 年度変化（MRI台数の推移） ──────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 共通の二次医療圏のみ追跡（2017〜2020で重複する sec_code）
common_secs = set(df[df['year']==2017]['sec_code']) & \
              set(df[df['year']==2018]['sec_code']) & \
              set(df[df['year']==2020]['sec_code'])
df_common = df[df['sec_code'].isin(common_secs)].copy()
pivot = df_common.pivot_table(index='sec_code', columns='year', values='mri_total')

ax = axes[0]
yr_totals = df.groupby('year')['mri_total'].sum()
ax.bar([str(y) for y in yr_totals.index], yr_totals.values,
       color=['#2196F3', '#4CAF50', '#FF9800'], alpha=0.85, edgecolor='white')
for i, (yr, v) in enumerate(yr_totals.items()):
    ax.text(i, v + 20, f'{v:,}台', ha='center', va='bottom', fontsize=11)
ax.set_ylabel('MRI台数合計', fontsize=11)
ax.set_title('全国MRI台数の推移（二次医療圏集計）', fontsize=11, fontweight='bold')

ax = axes[1]
change_17_20 = pivot[2020] - pivot[2017]
ax.hist(change_17_20.dropna(), bins=20, color='#7B1FA2', alpha=0.75, edgecolor='white')
ax.axvline(0, color='red', linewidth=1.5, linestyle='--')
ax.set_xlabel('MRI台数変化（2017→2020）', fontsize=11)
ax.set_ylabel('二次医療圏数', fontsize=11)
ax.set_title('MRI台数変化の分布（2017→2020）', fontsize=11, fontweight='bold')
mean_change = change_17_20.mean()
ax.axvline(mean_change, color='blue', linewidth=1.5, linestyle=':',
           label=f'平均変化 {mean_change:.1f}台')
ax.legend(fontsize=10)

fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2025_U4_trend.png'), bbox_inches='tight', dpi=150)
plt.close(fig)

print("Figure 5 保存: 2025_U4_trend.png")

# ── 記述統計サマリー ──────────────────────────────────────────────────
print("\n" + "=" * 55)
print("■ 記述統計（2020年）")
print("=" * 55)
desc_cols = ['mri_total', 'mri_3t', 'mri_15t', 'population', 'mri_per_100k']
print(d20[desc_cols].describe().round(2).to_string())

print(f"\n  [MRI密度] 全国平均: {avg_per100k:.2f}台/10万人")
print(f"  [最大]  {d20.loc[d20['mri_per_100k'].idxmax(), 'sec_name']}: "
      f"{d20['mri_per_100k'].max():.1f}台/10万人")
print(f"  [台数最大] {d20.loc[d20['mri_total'].idxmax(), 'sec_name']}: "
      f"{d20['mri_total'].max()}台")

print("\n■ 完了。全figureを html/figures/ に保存しました。")

a = pd.read_csv('data/raw/SSDSE-A-2025.csv', encoding='cp932', header=2)
a = a[a['地域コード'].str.match(r'^R\d+', na=False)].copy()
a = a[pd.to_numeric(a['総人口'], errors='coerce') > 0]
num = lambda c: pd.to_numeric(a[c], errors='coerce')

d2 = pd.DataFrame({'市区町村': (a['都道府県'] + a['市区町村']).values})
d2['病院密度'] = (num('一般病院数') / num('総人口') * 100000).values
d2['診療所密度'] = (num('一般診療所数') / num('総人口') * 100000).values
d2['高齢化率'] = (num('65歳以上人口') / num('総人口') * 100).values
d2['人口'] = num('総人口').values
d2 = d2.dropna()
big = d2[d2['人口'] >= 200000]; small = d2[d2['人口'] < 50000]
print(f'市区町村数: {len(d2)}')
print(f'人口20万以上 ({len(big)}): 病院密度 平均 {big["病院密度"].mean():.1f} / 10万人')
print(f'人口5万未満 ({len(small)}): 病院密度 平均 {small["病院密度"].mean():.1f} / 10万人')
t, p = stats.ttest_ind(big['病院密度'], small['病院密度'])
print(f'差の t 検定: p = {p:.4f}')
print('→ 小規模自治体ほど人口当たり病院数は多い（病院の最小規模性）。MRI の「過剰導入」を人口比だけで判定できない構造')