"""
教育用再現コード: 2024年 統計データ分析コンペティション 優秀賞（高校生）
=================================================================
論文タイトル：福祉支援を通じた過疎化対策の提案
著者：黒木喬士郎、井上和幸、髙山大綺、玉田章人（大分工業高等専門学校）

【分析概要】
  データ：SSDSE-B-2026（政府統計の総合窓口 / 統計数理研究所）

  Step1. 変数構築・記述統計
  Step2. 相関分析（多重共線性確認）
  Step3. VIF確認
  Step4. 重回帰分析①（目的変数：人口増加率）
  Step5. 重回帰分析②（目的変数：保育所利用率）

  Key findings:
    - 人口増加率 R² は標準地価・雇用状況・保育所整備状況で説明
    - 保育所利用率 R² は保育所数・雇用・人口規模で説明
    - 保育施設整備が地域の福祉水準と人口動態に関連

【データサイエンス学習ポイント】
  1. 目的変数を2つ持つ分析（人口増加率 vs 保育所利用率）
  2. VIFによる多重共線性の確認と対処
  3. 標準化回帰係数による変数間の影響力比較
  4. 相関ヒートマップによる変数関係の可視化
  5. 過疎化という政策課題への統計的アプローチ

【データ】SSDSE-B-2026.csv（実データ）を使用。合成データ不使用。
=================================================================
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・SSDSE-B-2026.csv
#       → data/raw/SSDSE-B-2026.csv に配置
#
#   ダウンロード先:
#     https://www.nstac.go.jp/use/literacy/ssdse/
#     （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード）
#
#   フォルダ配置（プロジェクトルートからの相対パス）:
#     code/                ← このスクリプトの場所
#     data/raw/            ← CSV ファイルをここに配置
#     html/figures/        ← 図の出力先（自動生成）
#
#   実行方法（ファイルを一切編集せず実行可能）:
#     python3 code/2024_H2_yushu.py
# ============================================================


import os
import warnings
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from scipy import stats

warnings.filterwarnings('ignore')

# ──────────────────────────────────────────────────────────────
# 共通設定
# ──────────────────────────────────────────────────────────────
plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

FIG_DIR = 'html/figures'
DATA_PATH = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

# ================================================================
# ■ Step 0. 実データ読み込み（SSDSE-B-2026）
# ================================================================

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

df_b = pd.read_csv(DATA_PATH, encoding='cp932', header=1)

# 都道府県行のみ抽出（地域コード R + 5桁数字）
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()

# 列名を変数コードへマッピング（実列名を使用）
COL = {
    'A1101':   '総人口',
    'A1302':   '15～64歳人口',
    'F3102':   '月間有効求職者数（一般）',
    'F3105':   '就職件数（一般）',
    'E2101':   '小学校数',
    'I510120': '一般病院数',
    'L3221':   '消費支出（二人以上の世帯）',
    'C5401':   '標準価格（平均価格）（住宅地）',
    'J2503':   '保育所等数',
    'J2506':   '保育所等在所児数',
}

# 必要列を数値化
for code, colname in COL.items():
    df_b[colname] = pd.to_numeric(df_b[colname], errors='coerce')

# ── 2022 年断面 ──────────────────────────────────────────────
df22 = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)
df21 = df_b[df_b['年度'] == 2021][['都道府県', '総人口']].rename(
    columns={'総人口': '総人口_2021'}
)

# ── 目的変数①：人口増加率（2021→2022）──────────────────────
df22 = df22.merge(df21, on='都道府県', how='left')
df22['pop_growth'] = (
    (df22['総人口'] - df22['総人口_2021']) / df22['総人口_2021'] * 100
)

# ── 目的変数②：保育所利用率（在所児数 / 総人口 × 1000）───────
df22['welfare_rate'] = df22['保育所等在所児数'] / df22['総人口'] * 1000

# ── 説明変数 ─────────────────────────────────────────────────
# 正規雇用率 proxy: 就職件数 / 月間有効求職者数 × 100
df22['empl_rate'] = df22['就職件数（一般）'] / df22['月間有効求職者数（一般）'] * 100

# 失業率 proxy: 月間有効求職者数 / 15-64歳人口 × 100
df22['unemp_rate'] = df22['月間有効求職者数（一般）'] / df22['15～64歳人口'] * 100

# 学校数（人口比）: 小学校数 / 総人口 × 10000
df22['school_per10k'] = df22['小学校数'] / df22['総人口'] * 10000

# 病院数（人口比）: 一般病院数 / 総人口 × 10000
df22['hospital_per10k'] = df22['一般病院数'] / df22['総人口'] * 10000

# 一人当たり消費支出（直接使用）
df22['consumption'] = df22['消費支出（二人以上の世帯）']

# 標準地価（住宅地）
df22['land_price'] = df22['標準価格（平均価格）（住宅地）']

# 保育所数（人口比）: 保育所等数 / 総人口 × 10000
df22['daycare_per10k'] = df22['保育所等数'] / df22['総人口'] * 10000

# ── 分析データフレーム構築 ────────────────────────────────────
EXPLAIN_VARS = [
    'empl_rate', 'unemp_rate', 'school_per10k', 'hospital_per10k',
    'consumption', 'land_price', 'daycare_per10k',
]
EXPLAIN_LABELS = {
    'empl_rate':       '正規雇用率（就職/求職）',
    'unemp_rate':      '失業率 proxy（求職/15-64歳）',
    'school_per10k':   '小学校数（1万人比）',
    'hospital_per10k': '一般病院数（1万人比）',
    'consumption':     '一人当たり消費支出',
    'land_price':      '標準地価（住宅地）',
    'daycare_per10k':  '保育所数（1万人比）',
}

TARGET_VARS = ['pop_growth', 'welfare_rate']
all_vars = TARGET_VARS + EXPLAIN_VARS
df_ana = df22[['都道府県'] + all_vars].dropna().reset_index(drop=True)

print(f"分析対象: {len(df_ana)} 都道府県（欠損除去後）")
print(f"年度: 2022年（人口増加率は 2021→2022 年変化率）")
print()
print("記述統計（目的変数）:")
print(df_ana[TARGET_VARS].rename(
    columns={'pop_growth': '人口増加率（%）', 'welfare_rate': '保育所利用率（‰）'}
).describe().round(4))

# ================================================================
# ■ Step 1. 標準化・VIF計算
# ================================================================

print("\n" + "=" * 65)
print("■ VIF（分散拡大係数）の確認")
print("=" * 65)

# 説明変数の標準化（z-score）
df_std = df_ana[EXPLAIN_VARS].copy()
for v in EXPLAIN_VARS:
    mu, sig = df_std[v].mean(), df_std[v].std(ddof=1)
    df_std[v + '_z'] = (df_std[v] - mu) / sig
Z_VARS = [v + '_z' for v in EXPLAIN_VARS]

X_vif = df_std[Z_VARS].values
vif_vals = {}
for i, zv in enumerate(Z_VARS):
    vif = variance_inflation_factor(X_vif, i)
    vif_vals[zv] = vif
    label = EXPLAIN_LABELS.get(zv.replace('_z', ''), zv)
    flag = '  ★多重共線性の疑い' if vif > 5 else ''
    print(f"  {label:<28} VIF = {vif:5.2f}{flag}")

# ================================================================
# ■ Step 2. 重回帰分析①（目的変数：人口増加率）
# ================================================================

print("\n" + "=" * 65)
print("■ 重回帰分析①：人口増加率（2021→2022）")
print("=" * 65)

X_reg = sm.add_constant(df_std[Z_VARS])
y1 = df_ana['pop_growth']
reg1 = sm.OLS(y1, X_reg).fit(cov_type='HC1')

print(reg1.summary())
print()
print(f"  R²          = {reg1.rsquared:.4f}")
print(f"  調整済み R² = {reg1.rsquared_adj:.4f}")
print(f"  F 統計量    = {reg1.fvalue:.3f}  (p = {reg1.f_pvalue:.4f})")
print()
print("  標準化回帰係数:")
for zv in Z_VARS:
    b = reg1.params[zv]
    p = reg1.pvalues[zv]
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '  '
    label = EXPLAIN_LABELS.get(zv.replace('_z', ''), zv)
    print(f"    {label:<28} β = {b:+.4f}  p = {p:.4f} {sig}")

# ================================================================
# ■ Step 3. 重回帰分析②（目的変数：保育所利用率）
# ================================================================

print("\n" + "=" * 65)
print("■ 重回帰分析②：保育所利用率（在所児数 / 総人口 × 1000）")
print("=" * 65)

y2 = df_ana['welfare_rate']
reg2 = sm.OLS(y2, X_reg).fit(cov_type='HC1')

print(reg2.summary())
print()
print(f"  R²          = {reg2.rsquared:.4f}")
print(f"  調整済み R² = {reg2.rsquared_adj:.4f}")
print(f"  F 統計量    = {reg2.fvalue:.3f}  (p = {reg2.f_pvalue:.4f})")
print()
print("  標準化回帰係数:")
for zv in Z_VARS:
    b = reg2.params[zv]
    p = reg2.pvalues[zv]
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '  '
    label = EXPLAIN_LABELS.get(zv.replace('_z', ''), zv)
    print(f"    {label:<28} β = {b:+.4f}  p = {p:.4f} {sig}")

# ================================================================
# ■ 図の生成（4枚）
# ================================================================

print("\n" + "=" * 65)
print("■ 図の生成（4枚）")
print("=" * 65)

short_labels = [EXPLAIN_LABELS[v] for v in EXPLAIN_VARS]
vif_labels   = [EXPLAIN_LABELS.get(zv.replace('_z', ''), zv) for zv in Z_VARS]

# ─── 図1: 全変数の相関ヒートマップ ──────────────────────────────
print("図1: 相関ヒートマップを作成中...")

all_corr_vars = TARGET_VARS + EXPLAIN_VARS
all_corr_labels = {
    'pop_growth':  '人口増加率',
    'welfare_rate': '保育所利用率',
    **EXPLAIN_LABELS,
}
corr_labels_list = [all_corr_labels[v] for v in all_corr_vars]
corr_mat = df_ana[all_corr_vars].corr()

fig1, ax1 = plt.subplots(figsize=(12, 10))
im = ax1.imshow(corr_mat.values, cmap='RdBu_r', vmin=-1, vmax=1)
plt.colorbar(im, ax=ax1, shrink=0.8, label='相関係数')
ax1.set_xticks(range(len(all_corr_vars)))
ax1.set_yticks(range(len(all_corr_vars)))
ax1.set_xticklabels(corr_labels_list, fontsize=8, rotation=35, ha='right')
ax1.set_yticklabels(corr_labels_list, fontsize=8)
for i in range(len(all_corr_vars)):
    for j in range(len(all_corr_vars)):
        v = corr_mat.values[i, j]
        col = 'white' if abs(v) > 0.6 else 'black'
        ax1.text(j, i, f'{v:.2f}', ha='center', va='center',
                 fontsize=7.5, color=col, fontweight='bold')
# 高相関ペアに枠
high_corr_pairs = [
    (i, j) for i in range(len(all_corr_vars))
    for j in range(i + 1, len(all_corr_vars))
    if abs(corr_mat.values[i, j]) > 0.7
]
for (i, j) in high_corr_pairs:
    ax1.add_patch(plt.Rectangle((j - 0.5, i - 0.5), 1, 1,
                                fill=False, edgecolor='gold', lw=2.5))
    ax1.add_patch(plt.Rectangle((i - 0.5, j - 0.5), 1, 1,
                                fill=False, edgecolor='gold', lw=2.5))
ax1.set_title('全変数間の相関係数ヒートマップ\n（目的変数・説明変数、金枠：|r|>0.7）',
              fontsize=12, fontweight='bold', pad=15)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2024_H2_fig1_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("  -> 2024_H2_fig1_corr.png 保存完了")

# ─── 図2: 両回帰モデルの標準化係数比較 ──────────────────────────
print("図2: 標準化回帰係数比較を作成中...")

coefs1 = [reg1.params[zv] for zv in Z_VARS]
coefs2 = [reg2.params[zv] for zv in Z_VARS]
pvals1 = [reg1.pvalues[zv] for zv in Z_VARS]
pvals2 = [reg2.pvalues[zv] for zv in Z_VARS]

x_pos = np.arange(len(Z_VARS))
width = 0.38

fig2, ax2 = plt.subplots(figsize=(13, 6))
bars1 = ax2.bar(x_pos - width / 2, coefs1, width,
                label=f'人口増加率モデル（R²={reg1.rsquared:.3f}）',
                color='#1565C0', alpha=0.82, edgecolor='white')
bars2 = ax2.bar(x_pos + width / 2, coefs2, width,
                label=f'保育所利用率モデル（R²={reg2.rsquared:.3f}）',
                color='#43A047', alpha=0.82, edgecolor='white')
# 有意な係数に *印
for bar, p in zip(bars1, pvals1):
    if p < 0.05:
        ax2.text(bar.get_x() + bar.get_width() / 2,
                 bar.get_height() + (0.01 if bar.get_height() >= 0 else -0.03),
                 '*', ha='center', va='bottom', fontsize=11, color='#1565C0', fontweight='bold')
for bar, p in zip(bars2, pvals2):
    if p < 0.05:
        ax2.text(bar.get_x() + bar.get_width() / 2,
                 bar.get_height() + (0.01 if bar.get_height() >= 0 else -0.03),
                 '*', ha='center', va='bottom', fontsize=11, color='#43A047', fontweight='bold')
ax2.set_xticks(x_pos)
ax2.set_xticklabels(short_labels, fontsize=9, rotation=25, ha='right')
ax2.axhline(0, color='black', linewidth=1.0)
ax2.set_ylabel('標準化回帰係数（β）', fontsize=11)
ax2.set_title('両回帰モデルの標準化回帰係数比較\n（*: p<0.05）', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(axis='y', alpha=0.3)
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2024_H2_fig2_coef.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("  -> 2024_H2_fig2_coef.png 保存完了")

# ─── 図3: VIF棒グラフ ────────────────────────────────────────────
print("図3: VIF棒グラフを作成中...")

vif_values_list = [vif_vals[zv] for zv in Z_VARS]
bar_colors_vif = [
    '#E53935' if v > 5 else '#FF8F00' if v > 3 else '#43A047'
    for v in vif_values_list
]

fig3, ax3 = plt.subplots(figsize=(9, 5))
ax3.barh(range(len(Z_VARS)), vif_values_list,
         color=bar_colors_vif, alpha=0.85, edgecolor='white')
ax3.axvline(5, color='#E53935', linestyle='--', linewidth=2, label='VIF=5（警戒線）')
ax3.axvline(3, color='#FF8F00', linestyle='--', linewidth=1.5, label='VIF=3（注意線）')
ax3.set_yticks(range(len(Z_VARS)))
ax3.set_yticklabels(vif_labels, fontsize=9)
ax3.set_xlabel('VIF（分散拡大係数）', fontsize=11)
ax3.set_title('VIF（多重共線性確認）\n緑: VIF≤3 / 橙: VIF>3 / 赤: VIF>5',
              fontsize=11, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(axis='x', alpha=0.3)
ax3.invert_yaxis()
for i, v in enumerate(vif_values_list):
    ax3.text(v + 0.05, i, f'{v:.2f}', va='center', fontsize=8.5)
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2024_H2_fig3_vif.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("  -> 2024_H2_fig3_vif.png 保存完了")

# ─── 図4: 保育所数 vs 保育所利用率（散布図＋回帰線）────────────
print("図4: 散布図（保育所数 vs 保育所利用率）を作成中...")

x4 = df_ana['daycare_per10k']
y4 = df_ana['welfare_rate']
slope, intercept, r_val, p_val, se = stats.linregress(x4, y4)
x_line = np.linspace(x4.min(), x4.max(), 100)
y_line = intercept + slope * x_line

fig4, ax4 = plt.subplots(figsize=(8, 6))
ax4.scatter(x4, y4, color='#1565C0', alpha=0.75, s=70,
            edgecolors='white', linewidth=0.5, zorder=3)
ax4.plot(x_line, y_line, color='#E53935', linewidth=2,
         label=f'回帰直線 (r={r_val:.3f}, p={p_val:.4f})')
# 都道府県名ラベル（上位・下位）
for i, row in df_ana.iterrows():
    if row['daycare_per10k'] > x4.quantile(0.85) or row['daycare_per10k'] < x4.quantile(0.15):
        ax4.annotate(row['都道府県'], (row['daycare_per10k'], row['welfare_rate']),
                     fontsize=7.5, xytext=(4, 3), textcoords='offset points', alpha=0.8)
ax4.set_xlabel('保育所等数（1万人比）', fontsize=12)
ax4.set_ylabel('保育所等利用率（在所児数 / 総人口 × 1000）', fontsize=12)
ax4.set_title('保育所数と保育所利用率の関係\n（2022年、47都道府県）',
              fontsize=12, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2024_H2_fig4_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("  -> 2024_H2_fig4_scatter.png 保存完了")

print("\n" + "=" * 65)
print("全図の生成完了（4枚）")
print("=" * 65)
print(f"\n保存先: {os.path.abspath(FIG_DIR)}")
print("  2024_H2_fig1_corr.png    - 全変数の相関ヒートマップ")
print("  2024_H2_fig2_coef.png    - 両回帰モデルの標準化係数比較")
print("  2024_H2_fig3_vif.png     - VIF棒グラフ")
print("  2024_H2_fig4_scatter.png - 保育所数 vs 保育所利用率")