日本における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/ に保存しました。")