欠測値を含むパネルデータを用いた
健康寿命の要因分析

print("=" * 65)
print("■ SSDSE-B-2026 データ読み込み")
print("=" * 65)

df = pd.read_csv(DATA_B, encoding='cp932', header=1)
df = df[df['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()

print(f"都道府県数: {df['都道府県'].nunique()}")
print(f"年度範囲  : {df['年度'].min()}〜{df['年度'].max()}")
print(f"総行数    : {len(df)} (47都道府県 × 12年)")

df = df.sort_values(['都道府県', '年度']).reset_index(drop=True)

df['保健医療費割合'] = df['保健医療費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）'] * 100

# 説明変数
df['高齢化率']   = df['65歳以上人口'] / df['総人口'] * 100
df['食料費割合'] = df['食料費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）'] * 100
df['病院密度']   = df['一般病院数'] / df['総人口'] * 100000

print("\n■ 主要変数 基本統計")
stats_cols = ['保健医療費割合', '高齢化率', '食料費割合', '病院密度']
print(df[stats_cols].describe().round(3).to_string())

print("\n" + "=" * 65)
print("■ 図1：欠測値可視化")
print("=" * 65)

# 分析に関連する変数群の欠測カウント
vis_cols = [
    '保健医療費（二人以上の世帯）',
    '消費支出（二人以上の世帯）',
    '食料費（二人以上の世帯）',
    '住居費（二人以上の世帯）',
    '光熱・水道費（二人以上の世帯）',
    '教育費（二人以上の世帯）',
    '教養娯楽費（二人以上の世帯）',
    '総人口',
    '65歳以上人口',
    '15歳未満人口',
    '一般病院数',
    '一般診療所数',
    '保育所等数',
    '合計特殊出生率',
    '標準価格（平均価格）（住宅地）',
    '標準価格（平均価格）（商業地）',
    '大学数',
    '保育所等利用待機児童数',
]

miss_counts = df[vis_cols].isnull().sum()
n_total = len(df)

# 変数名ラベル（短縮版）
short_labels = {
    '保健医療費（二人以上の世帯）'   : '保健医療費',
    '消費支出（二人以上の世帯）'     : '消費支出',
    '食料費（二人以上の世帯）'       : '食料費',
    '住居費（二人以上の世帯）'       : '住居費',
    '光熱・水道費（二人以上の世帯）' : '光熱・水道費',
    '教育費（二人以上の世帯）'       : '教育費',
    '教養娯楽費（二人以上の世帯）'   : '教養娯楽費',
    '総人口'                        : '総人口',
    '65歳以上人口'                  : '65歳以上人口',
    '15歳未満人口'                  : '15歳未満人口',
    '一般病院数'                    : '一般病院数',
    '一般診療所数'                  : '一般診療所数',
    '保育所等数'                    : '保育所等数',
    '合計特殊出生率'                : '合計特殊出生率',
    '標準価格（平均価格）（住宅地）' : '住宅地価格',
    '標準価格（平均価格）（商業地）' : '商業地価格',
    '大学数'                        : '大学数',
    '保育所等利用待機児童数'        : '待機児童数',
}

labels_short = [short_labels[c] for c in vis_cols]
miss_pcts = miss_counts.values / n_total * 100

# 色：欠測ゼロ＝青緑（完全）、欠測あり＝オレンジ
bar_colors_miss = ['#78909C' if v == 0 else '#E65100' for v in miss_counts.values]

fig1, ax1 = plt.subplots(figsize=(11, 5.5))
bars1 = ax1.bar(np.arange(len(vis_cols)), miss_pcts,
                color=bar_colors_miss, alpha=0.85, edgecolor='white', width=0.65)

ax1.set_xticks(np.arange(len(vis_cols)))
ax1.set_xticklabels(labels_short, rotation=45, ha='right', fontsize=9)
ax1.set_ylabel('欠測率（%）', fontsize=11)
ax1.set_title('SSDSE-B-2026（都道府県別パネルデータ）各変数の欠測状況\n'
              f'全{n_total}観測（47都道府県 × 12年）',
              fontsize=12, fontweight='bold')
ax1.set_ylim(0, max(miss_pcts.max() + 5, 15))
ax1.axhline(5, color='#C62828', linestyle='--', linewidth=1.0, alpha=0.7,
            label='欠測率 5% 目安')
ax1.grid(axis='y', alpha=0.3)

# 欠測数・割合ラベル
for i, (v, pct) in enumerate(zip(miss_counts.values, miss_pcts)):
    label = f'{v}件\n({pct:.1f}%)' if v > 0 else '0件\n(完全)'
    ax1.text(i, max(pct + 0.3, 0.5), label,
             ha='center', va='bottom', fontsize=7.5, color='#333')

# 凡例
patch_ok  = mpatches.Patch(color='#78909C', alpha=0.85, label='欠測なし（完全観測）')
patch_ng  = mpatches.Patch(color='#E65100', alpha=0.85, label='欠測あり')
ax1.legend(handles=[patch_ok, patch_ng], fontsize=9, loc='upper right')

# 教育的注釈
ax1.text(0.01, 0.97,
         '★ 欠測パターンの確認は前処理の第一歩（MCAR/MAR/MNARの判断に必要）',
         transform=ax1.transAxes, fontsize=8.5, color='#1565C0',
         va='top', style='italic')

plt.tight_layout()
fig1_path = os.path.join(FIG_DIR, '2022_U5_8_fig1_missing.png')
fig1.savefig(fig1_path, bbox_inches='tight', dpi=150)
plt.close(fig1)
print(f"図1保存: {fig1_path}")

print("\n■ 図2：保健医療費割合 時系列")

nat_avg_ts = df.groupby('年度')['保健医療費割合'].mean()

# 2022年の値で上位5・下位5県を選定
df_2022 = df[df['年度'] == 2022].copy()
top5_prefs  = df_2022.nlargest(5,  '保健医療費割合')['都道府県'].tolist()
bot5_prefs  = df_2022.nsmallest(5, '保健医療費割合')['都道府県'].tolist()

# 色マップ
cmap_top = plt.cm.Reds(np.linspace(0.5, 0.9, 5))
cmap_bot = plt.cm.Blues(np.linspace(0.5, 0.9, 5))
color_map = {}
for i, p in enumerate(top5_prefs):
    color_map[p] = cmap_top[i]
for i, p in enumerate(bot5_prefs):
    color_map[p] = cmap_bot[i]

highlight_prefs = top5_prefs + bot5_prefs

fig2, ax2 = plt.subplots(figsize=(11, 5.5))

# 全47都道府県（グレー背景）
for pref in df['都道府県'].unique():
    if pref not in highlight_prefs:
        tmp = df[df['都道府県'] == pref].sort_values('年度')
        ax2.plot(tmp['年度'], tmp['保健医療費割合'],
                 color='#BDBDBD', linewidth=0.6, alpha=0.40, zorder=1)

# ハイライト：高位5県（赤系）
for i, pref in enumerate(top5_prefs):
    tmp = df[df['都道府県'] == pref].sort_values('年度')
    val = df_2022[df_2022['都道府県'] == pref]['保健医療費割合'].values[0]
    ax2.plot(tmp['年度'], tmp['保健医療費割合'],
             color=color_map[pref], linewidth=1.8, marker='o', markersize=3.5,
             label=f'{pref}（{val:.1f}%）', zorder=4)

# ハイライト：低位5県（青系）
for i, pref in enumerate(bot5_prefs):
    tmp = df[df['都道府県'] == pref].sort_values('年度')
    val = df_2022[df_2022['都道府県'] == pref]['保健医療費割合'].values[0]
    ax2.plot(tmp['年度'], tmp['保健医療費割合'],
             color=color_map[pref], linewidth=1.8, marker='s', markersize=3.5,
             label=f'{pref}（{val:.1f}%）', zorder=4)

# 全国平均（太線・黒破線）
ax2.plot(nat_avg_ts.index, nat_avg_ts.values,
         color='#212121', linewidth=2.5, linestyle='--',
         label=f'全国平均（2022: {nat_avg_ts[2022]:.1f}%）', zorder=5)

ax2.set_xlabel('年度', fontsize=12)
ax2.set_ylabel('保健医療費割合（%）', fontsize=12)
ax2.set_title('保健医療費割合の時系列推移\n全47都道府県（2012〜2023年度）— 上位5県（赤）・下位5県（青）・全国平均（黒破線）',
              fontsize=12, fontweight='bold')
ax2.set_xticks(sorted(df['年度'].unique()))
ax2.tick_params(axis='x', rotation=45)
ax2.grid(True, alpha=0.3)
ax2.legend(fontsize=8, loc='upper left', ncol=2, framealpha=0.9)

# COVID注釈
ymin, ymax = ax2.get_ylim()
ax2.axvspan(2019.5, 2021.5, color='#FFECB3', alpha=0.35, zorder=0)
ax2.text(2020.5, ymin + (ymax - ymin) * 0.04,
         'COVID-19', ha='center', fontsize=8, color='#795548', style='italic')

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

print("\n■ 図3：FE係数プロット")

coefs  = [fe_params[v]  for v in XVARS]
ses    = [fe_stderr[v]  for v in XVARS]
pvals  = [fe_pvalues[v] for v in XVARS]
labels = [var_labels[v] for v in XVARS]

bar_colors_fe = []
for p in pvals:
    if p < 0.01:
        bar_colors_fe.append('#C62828')
    elif p < 0.05:
        bar_colors_fe.append('#FF8F00')
    else:
        bar_colors_fe.append('#9E9E9E')

fig3, ax3 = plt.subplots(figsize=(9, 4.5))
ypos = np.arange(len(XVARS))

ax3.barh(ypos, coefs, color=bar_colors_fe, alpha=0.82,
         edgecolor='white', height=0.5)
ax3.errorbar(coefs, ypos, xerr=1.96 * np.array(ses),
             fmt='none', color='#333', capsize=5, linewidth=1.8)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(ypos)
ax3.set_yticklabels(labels, fontsize=12)
ax3.set_xlabel('固定効果モデル 推定係数（±1.96 SE）', fontsize=11)
ax3.set_title(f'保健医療費割合の決定要因（パネルOLS 双方向固定効果モデル）\n'
              f'Within R²={fe_rsq:.3f}、N=47都道府県×12年（クラスタリング標準誤差）',
              fontsize=12, fontweight='bold')
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)

# p値注釈
x_range = max(coefs) - min(coefs) if len(coefs) > 1 else abs(coefs[0]) * 2
margin = max(x_range * 0.04, 0.001)
for i, (c, p, se) in enumerate(zip(coefs, pvals, ses)):
    sig = ('***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.')
    x_ann = c + 1.96 * se + margin if c >= 0 else c - 1.96 * se - margin
    ha = 'left' if c >= 0 else 'right'
    ax3.text(x_ann, i, sig, va='center', ha=ha,
             fontsize=10, color=bar_colors_fe[i], fontweight='bold')

legend_patches = [
    mpatches.Patch(color='#C62828', alpha=0.82, label='p < 0.01'),
    mpatches.Patch(color='#FF8F00', alpha=0.82, label='p < 0.05'),
    mpatches.Patch(color='#9E9E9E', alpha=0.82, label='n.s. (p ≥ 0.05)'),
]
ax3.legend(handles=legend_patches, fontsize=9, loc='lower right')

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

import matplotlib
matplotlib.use('Agg')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib import rcParams
from linearmodels import PanelOLS
import warnings
warnings.filterwarnings('ignore')
import os

rcParams['font.family'] = 'Hiragino Sans'
rcParams['axes.unicode_minus'] = False
rcParams['figure.dpi'] = 150

_dir = os.path.dirname(os.path.abspath(__file__))
DATA_B  = os.path.join(_dir, '..', 'data', 'raw', 'SSDSE-B-2026.csv')
FIG_DIR = os.path.join(_dir, '..', 'html', 'figures')
os.makedirs(FIG_DIR, exist_ok=True)

print("\n" + "=" * 65)
print("■ パネルOLS 固定効果モデル推定（双方向固定効果）")
print("=" * 65)

YVAR  = '保健医療費割合'
XVARS = ['高齢化率', '食料費割合', '病院密度']

df_panel = df.copy().set_index(['都道府県', '年度'])
panel_data = df_panel[[YVAR] + XVARS].dropna()

print(f"推定サンプル: {len(panel_data)} obs（{panel_data.index.get_level_values(0).nunique()}都道府県）")

mod = PanelOLS(
    dependent=panel_data[YVAR],
    exog=panel_data[XVARS],
    entity_effects=True,
    time_effects=True,
)
res = mod.fit(cov_type='clustered', cluster_entity=True)
print(res.summary)

fe_params  = res.params
fe_stderr  = res.std_errors
fe_pvalues = res.pvalues
fe_rsq     = res.rsquared

print("\n■ 推定結果まとめ")
print(f"  Within R²: {fe_rsq:.4f}")
var_labels = {
    '高齢化率'   : '高齢化率（%）',
    '食料費割合' : '食料費割合（%）',
    '病院密度'   : '病院密度（人口10万対）',
}
for var in XVARS:
    sig = ('***' if fe_pvalues[var] < 0.001 else
           '**'  if fe_pvalues[var] < 0.01  else
           '*'   if fe_pvalues[var] < 0.05  else 'n.s.')
    print(f"  {var:<12} β={fe_params[var]:+.4f}  SE={fe_stderr[var]:.4f}  "
          f"p={fe_pvalues[var]:.4f}  {sig}")

print("\n■ 図4：高齢化率 vs 保健医療費割合 散布図（2022年）")

df_2022 = df[df['年度'] == 2022].copy()

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

fig4, ax4 = plt.subplots(figsize=(10, 7))

# 地方別散布
for region, grp in df_2022.groupby('地方'):
    ax4.scatter(grp['高齢化率'], grp['保健医療費割合'],
                color=region_colors_map.get(region, '#999'),
                s=70, alpha=0.75, zorder=3, label=region,
                edgecolors='white', linewidths=0.8)

# 都道府県ラベル（全47）
for _, row in df_2022.iterrows():
    ax4.annotate(
        row['都道府県'].replace('県', '').replace('都', '').replace('道', '').replace('府', ''),
        xy=(row['高齢化率'], row['保健医療費割合']),
        xytext=(3, 3), textcoords='offset points',
        fontsize=6.5, color='#444', zorder=5,
    )

# 回帰直線
x_v = df_2022['高齢化率'].values
y_v = df_2022['保健医療費割合'].values
z   = np.polyfit(x_v, y_v, 1)
xs  = np.linspace(x_v.min(), x_v.max(), 200)
ax4.plot(xs, np.poly1d(z)(xs), color='#333', linewidth=1.8,
         linestyle='--', alpha=0.8, zorder=4, label='回帰直線')

from scipy import stats as sp_stats
r_val, p_val = sp_stats.pearsonr(x_v, y_v)

ax4.set_xlabel('高齢化率（%）', fontsize=12)
ax4.set_ylabel('保健医療費割合（%）', fontsize=12)
ax4.set_title(f'高齢化率 vs 保健医療費割合（2022年度・47都道府県）\n'
              f'ピアソン相関: r = {r_val:.3f}（p = {p_val:.4f}）',
              fontsize=12, fontweight='bold')
ax4.legend(fontsize=9, loc='upper right', framealpha=0.9)
ax4.grid(True, alpha=0.3)

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

print("\n" + "=" * 65)
print("■ 全図の生成完了（4枚）")
print("=" * 65)
print(f"  2022_U5_8_fig1_missing.png : 欠測値棒グラフ（変数別欠測率）")
print(f"  2022_U5_8_fig2_ts.png      : 保健医療費割合 時系列（全国 + 高低5県）")
print(f"  2022_U5_8_fig3_fe_coef.png : FE係数プロット（保健医療費割合の決定要因）")
print(f"  2022_U5_8_fig4_scatter.png : 高齢化率 vs 保健医療費割合 散布図（2022年）")
print(f"\n  パネルOLS結果（双方向固定効果、クラスタリング標準誤差）：")
for var in XVARS:
    sig = ('***' if fe_pvalues[var] < 0.001 else
           '**'  if fe_pvalues[var] < 0.01  else
           '*'   if fe_pvalues[var] < 0.05  else 'n.s.')
    print(f"    {var_labels[var]:<22} β={fe_params[var]:+.4f}  "
          f"p={fe_pvalues[var]:.4f}  {sig}")
print(f"\n  Within R² = {fe_rsq:.4f}")