XAIを用いた介護業界における
地域別の従業者数の就業要因に関する一考察

import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import warnings
warnings.filterwarnings('ignore')

from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import cross_val_score
from sklearn.inspection import partial_dependence
from scipy import stats as scipy_stats

try:
    import shap
    HAS_SHAP = True
except ImportError:
    HAS_SHAP = False
    print("SHAPが利用不可 → 代替実装（permutation importance近似）を使用")

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

print("=" * 60)
print("■ 実データ読み込み（SSDSE-E-2026 / SSDSE-B-2026）")
print("=" * 60)

# SSDSE-E
df_e_raw = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-E-2026.csv'),
                        encoding='cp932', header=0)
df_e = df_e_raw.iloc[1:].copy()
df_e.columns = df_e_raw.iloc[0].values
df_e = df_e.iloc[1:].copy()
df_e.columns = df_e_raw.iloc[1].values
df_e = df_e[df_e['都道府県'] != '全国'].set_index('都道府県').copy()

num_cols_e = ['総人口', '65歳以上人口', '総面積（北方地域及び竹島を除く）',
              '1人当たり県民所得（平成27年基準）', '医師数',
              '従業者数（民営）（医療、福祉）']
for c in num_cols_e:
    df_e[c] = pd.to_numeric(df_e[c], errors='coerce')

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

for c in ['総人口', '65歳以上人口', '年平均気温', '婚姻件数']:
    df_b[c] = pd.to_numeric(df_b[c], errors='coerce')

# ─ 変数作成 ─
common_prefs = sorted(set(df_e.index) & set(df_b.index))
PREFS = common_prefs

# 目的変数：医療・福祉従業者密度（人口1万対）
care_density = (df_e.loc[PREFS, '従業者数（民営）（医療、福祉）'] /
                df_e.loc[PREFS, '総人口'] * 10000).values.astype(float)

# 説明変数
aging_rate = (df_e.loc[PREFS, '65歳以上人口'] /
              df_e.loc[PREFS, '総人口']).values.astype(float)
income     = df_e.loc[PREFS, '1人当たり県民所得（平成27年基準）'].values.astype(float)
area       = df_e.loc[PREFS, '総面積（北方地域及び竹島を除く）'].values.astype(float)
pop        = df_e.loc[PREFS, '総人口'].values.astype(float)
pop_density = pop / (area / 100)  # /km²
doctor_rate = (df_e.loc[PREFS, '医師数'] / pop * 100000).values.astype(float)
marriage    = (df_b.loc[PREFS, '婚姻件数'] / df_b.loc[PREFS, '総人口'] * 1000).values.astype(float)
temp        = df_b.loc[PREFS, '年平均気温'].values.astype(float)

FEATURE_NAMES = ['高齢化率', '県民所得', '人口密度（対数）', '医師数10万対', '婚姻率', '年平均気温']
X = np.column_stack([aging_rate, income, np.log1p(pop_density), doctor_rate, marriage, temp])
n_features = len(FEATURE_NAMES)

# 欠損除去
valid = ~np.any(np.isnan(X), axis=1) & ~np.isnan(care_density)
X = X[valid]
care_density = care_density[valid]
PREFS_V = [PREFS[i] for i in range(len(PREFS)) if valid[i]]
N = len(PREFS_V)

print(f"分析対象: {N}都道府県")
print(f"医療・福祉従業者密度: mean={care_density.mean():.1f}, std={care_density.std():.1f} (/万人)")

print("\n図1: RF特徴量重要度を作成中...")

fig1, axes1 = plt.subplots(1, 2, figsize=(13, 5))
fig1.suptitle('Random Forest 特徴量重要度と累積重要度\n（医療・福祉従業者密度 / データ：SSDSE実データ）',
              fontsize=12, fontweight='bold')

imp_sorted = imp_df.sort_values('重要度')
top_feat = imp_df.iloc[0]['変数']
colors1 = ['#1565C0' if f == top_feat else '#2E7D32' if imp_df[imp_df['変数']==f]['重要度'].values[0] > 0.15
            else '#90CAF9' for f in imp_sorted['変数']]
axes1[0].barh(imp_sorted['変数'], imp_sorted['重要度'], color=colors1, edgecolor='white', alpha=0.88)
axes1[0].set_xlabel('特徴量重要度（不純度減少の平均）', fontsize=11)
axes1[0].set_title('Random Forest 特徴量重要度\n青：最重要変数', fontsize=11, fontweight='bold')
axes1[0].grid(axis='x', alpha=0.3)
for bar, val in zip(axes1[0].patches, imp_sorted['重要度']):
    axes1[0].text(val + 0.003, bar.get_y() + bar.get_height()/2,
                  f'{val:.3f}', va='center', fontsize=9)

imp_desc = imp_df.sort_values('重要度', ascending=False)
cumsum = imp_desc['重要度'].cumsum().values
x_tick = range(1, n_features + 1)
axes1[1].bar(x_tick, imp_desc['重要度'].values, color='#1565C0', alpha=0.6, edgecolor='white', label='重要度')
axes1[1].plot(x_tick, cumsum, 'ro-', linewidth=2, markersize=8, label='累積重要度')
axes1[1].axhline(0.8, color='gray', linestyle='--', linewidth=1.0, label='累積80%')
axes1[1].set_xticks(list(x_tick))
axes1[1].set_xticklabels([f.split('（')[0] for f in imp_desc['変数']], fontsize=8.5, rotation=15, ha='right')
axes1[1].set_ylabel('重要度', fontsize=11)
axes1[1].set_title('累積特徴量重要度', fontsize=11, fontweight='bold')
axes1[1].legend(fontsize=9)
axes1[1].grid(axis='y', alpha=0.3)

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

print("図2: SHAP summary plotを作成中...")

if HAS_SHAP:
    fig2, ax2 = plt.subplots(figsize=(10, 7))
    shap.summary_plot(shap_values, X, feature_names=FEATURE_NAMES,
                      show=False, plot_type='dot', color_bar=True)
    plt.title('SHAP Summary Plot（Beeswarm）\n医療・福祉従業者密度への各変数の影響',
              fontsize=13, fontweight='bold')
    plt.tight_layout()
    fig2.savefig(os.path.join(FIG_DIR, '2024_U5_6_fig2_shap_summary.png'), bbox_inches='tight', dpi=150)
    plt.close(fig2)
else:
    fig2, ax2 = plt.subplots(figsize=(10, 7))
    mean_abs_shap = np.abs(shap_values).mean(axis=0)
    order = np.argsort(mean_abs_shap)

    cmap = cm.get_cmap('RdBu_r')
    for plot_idx, feat_idx in enumerate(order):
        shap_col = shap_values[:, feat_idx]
        feat_col = X[:, feat_idx]
        feat_norm = (feat_col - feat_col.min()) / (feat_col.max() - feat_col.min() + 1e-8)
        y_jitter = np.linspace(-0.15, 0.15, N) + plot_idx
        ax2.scatter(shap_col, y_jitter, c=feat_norm, cmap='RdBu_r',
                    s=35, alpha=0.75, vmin=0, vmax=1)

    ax2.set_yticks(range(n_features))
    ax2.set_yticklabels([FEATURE_NAMES[i] for i in order], fontsize=10)
    ax2.axvline(0, color='black', linewidth=0.8)
    ax2.set_xlabel('SHAP値（医療・福祉従業者密度への影響）', fontsize=12)
    ax2.set_title('SHAP Summary Plot（代替実装）\n右=従業者密度増加に寄与、左=減少に寄与',
                  fontsize=12, fontweight='bold')
    ax2.grid(axis='x', alpha=0.2)
    sm_plot = plt.cm.ScalarMappable(cmap='RdBu_r', norm=plt.Normalize(0, 1))
    cbar = plt.colorbar(sm_plot, ax=ax2)
    cbar.set_label('特徴量値（低→高）', fontsize=10)

    plt.tight_layout()
    fig2.savefig(os.path.join(FIG_DIR, '2024_U5_6_fig2_shap_summary.png'), bbox_inches='tight', dpi=150)
    plt.close(fig2)

print("  → 2024_U5_6_fig2_shap_summary.png 保存完了")

rf = RandomForestRegressor(n_estimators=200, max_depth=6, min_samples_leaf=2,
                            random_state=0, n_jobs=-1)
rf.fit(X, care_density)

cv_r2 = cross_val_score(rf, X, care_density, cv=5, scoring='r2')
print(f"\n【Random Forest】")
print(f"  訓練R² = {rf.score(X, care_density):.3f}")
print(f"  5-fold CV R² = {cv_r2.mean():.3f} (±{cv_r2.std():.3f})")

importance = rf.feature_importances_
imp_df = pd.DataFrame({'変数': FEATURE_NAMES, '重要度': importance}).sort_values('重要度', ascending=False)
print(f"\n【特徴量重要度】")
print(imp_df.round(4))

if HAS_SHAP:
    explainer = shap.TreeExplainer(rf)
    shap_values = explainer.shap_values(X)
    shap_available = True
    print("\n■ SHAP値を計算しました")
else:
    y_pred = rf.predict(X)
    y_mean = y_pred.mean()
    shap_values = np.zeros((N, n_features))
    for i in range(N):
        deviation = y_pred[i] - y_mean
        for j in range(n_features):
            shap_values[i, j] = deviation * importance[j]
    shap_available = False
    print("\n■ SHAP近似値（代替実装）を計算しました")

print("図3: 部分従属プロットを作成中...")

fig3, axes3 = plt.subplots(1, 2, figsize=(13, 5))
fig3.suptitle('部分従属プロット（PDP）：非線形効果の可視化', fontsize=13, fontweight='bold')

pdp_features = [0, 3]  # 高齢化率, 医師数10万対
pdp_labels = ['高齢化率（65歳以上割合）', '医師数10万対']
pdp_colors = ['#1565C0', '#E65100']

for ax, feat_idx, label, clr in zip(axes3, pdp_features, pdp_labels, pdp_colors):
    pdp_result = partial_dependence(rf, X, features=[feat_idx], grid_resolution=50)
    grid_vals = pdp_result['grid_values'][0]
    pdp_vals = pdp_result['average'][0]

    ax.plot(grid_vals, pdp_vals, color=clr, linewidth=2.5)
    ax.fill_between(grid_vals, pdp_vals - pdp_vals.std() * 0.5,
                    pdp_vals + pdp_vals.std() * 0.5, alpha=0.15, color=clr)
    ax.scatter(X[:, feat_idx], rf.predict(X),
               c='gray', alpha=0.3, s=15, zorder=2, label='実データ（予測値）')
    ax.set_xlabel(label, fontsize=11)
    ax.set_ylabel('医療・福祉従業者密度（予測値）', fontsize=11)
    ax.set_title(f'部分従属プロット：{label}', fontsize=11, fontweight='bold')
    ax.grid(True, alpha=0.2)
    ax.legend(fontsize=9)

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

print("図4: 高齢化率 vs 医療・福祉従業者密度散布図を作成中...")

aging_vals = X[:, 0] * 100  # %表示
r_val, p_val = scipy_stats.pearsonr(aging_vals, care_density)

fig4, axes4 = plt.subplots(1, 2, figsize=(13, 5))
fig4.suptitle(f'高齢化率と医療・福祉従業者密度（XAI最重要変数：{top_feat}）',
              fontsize=13, fontweight='bold')

# 左: 散布図（県民所得で色分け）
ax4a = axes4[0]
income_norm = (X[:, 1] - X[:, 1].min()) / (X[:, 1].max() - X[:, 1].min())
sc4 = ax4a.scatter(aging_vals, care_density, c=income_norm, cmap='RdYlBu',
                   s=60, alpha=0.85, edgecolors='white', linewidth=0.5)
plt.colorbar(sc4, ax=ax4a, label='県民所得（低→高）')
coef4 = np.polyfit(aging_vals, care_density, 1)
x_fit = np.linspace(aging_vals.min(), aging_vals.max(), 100)
ax4a.plot(x_fit, np.polyval(coef4, x_fit), 'k--', linewidth=2)
# 代表的な都道府県ラベル
care_rank = np.argsort(care_density)[::-1]
for idx in list(care_rank[:5]) + list(care_rank[-5:]):
    short = PREFS_V[idx].replace('県','').replace('府','').replace('都','').replace('道','')
    ax4a.annotate(short, (aging_vals[idx], care_density[idx]),
                  textcoords='offset points', xytext=(5, 3), fontsize=7.5, color='#333')
ax4a.set_xlabel('高齢化率（65歳以上割合 %）', fontsize=11)
ax4a.set_ylabel('医療・福祉従業者密度（/万人）', fontsize=11)
ax4a.set_title(f'高齢化率 → 医療・福祉従業者密度\nr = {r_val:.3f}', fontsize=11, fontweight='bold')
ax4a.grid(True, alpha=0.2)
ax4a.text(0.05, 0.95, f'r = {r_val:.3f}\n（SHAP最重要変数）',
          transform=ax4a.transAxes, fontsize=10, va='top',
          bbox=dict(boxstyle='round', facecolor='#E3F2FD', alpha=0.8))

# 右: SHAP値（高齢化率）vs 高齢化率
ax4b = axes4[1]
shap_aging = shap_values[:, 0]
ax4b.scatter(aging_vals, shap_aging, c='#1565C0', alpha=0.75, s=55)
coef4b = np.polyfit(aging_vals, shap_aging, 1)
ax4b.plot(x_fit, np.polyval(coef4b, x_fit), 'r--', linewidth=2, label='回帰直線')
ax4b.axhline(0, color='black', linewidth=0.8)
ax4b.set_xlabel('高齢化率（65歳以上割合 %）', fontsize=11)
ax4b.set_ylabel('SHAP値（医療・福祉従業者密度への影響）', fontsize=11)
ax4b.set_title('SHAP値 vs 高齢化率\n非線形効果の確認', fontsize=11, fontweight='bold')
ax4b.legend(fontsize=9)
ax4b.grid(True, alpha=0.2)

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

print("\n" + "=" * 60)
print("✓ 全図の生成完了（4枚）")
print("=" * 60)
print("\n【主要知見】")
print(f"  RF訓練R² = {rf.score(X, care_density):.3f}")
print(f"  5-fold CV R² = {cv_r2.mean():.3f}")
print(f"  最重要変数: {imp_df.iloc[0]['変数']} (重要度={imp_df.iloc[0]['重要度']:.3f})")
print(f"  高齢化率と医療・福祉従業者密度の相関: r = {r_val:.3f}")
print(f"  SHAPライブラリ利用: {HAS_SHAP}")
print(f"  使用データ: SSDSE-E-2026, SSDSE-B-2026 ({YEAR}年, {N}都道府県)")