"""
教育用再現コード: 2024年 統計データ分析コンペティション 審査員奨励賞（大学生・一般）
=================================================================================
論文タイトル：小学生の運動能力と環境要因の分析
              ―スポーツ活動参加率の男女別規定要因―
受賞区分  ：審査員奨励賞 ［大学生・一般の部］
著者      ：中原智哉、山崎柊丞（早稲田大学基幹理工学部）

【分析概要】
  データ：SSDSE-D-2023（社会生活基本調査：生活時間・行動者率）+
          SSDSE-B-2026（都道府県別）+ SSDSE-E-2026（クロスセクション）
  目的   ：相関分析→VIF確認→男女別重回帰で
            スポーツ活動参加率の地域差を規定する要因を特定する

  Step1. 実データ読み込み（SSDSE-D 男女別 47都道府県）
  Step2. 相関分析（Pearson相関係数ヒートマップ）
  Step3. VIF（分散拡大係数）計算
  Step4. 重回帰分析（男女別）

【変数説明】
  目的変数: スポーツの総数（行動者率%）… 男別・女別
  説明変数:
    - 保育所千対    : 保育所等数 / 総人口 × 10000（SSDSE-B）
    - 年平均気温    : SSDSE-B
    - 県民所得      : 1人当たり県民所得（SSDSE-E）
    - 高齢化率      : 65歳以上人口 / 総人口（SSDSE-B）
    - ウォーキング率: ウォーキング・軽い体操 行動者率（SSDSE-D）
    - 睡眠時間      : 睡眠（分）（SSDSE-D）

【データ出典】
  独立行政法人統計センター「SSDSE（教育用標準データセット）」
  https://www.nstac.go.jp/use/literacy/ssdse/

【図の出力】
  html/figures/2024_U5_5_fig1_corr.png          ... 相関係数ヒートマップ
  html/figures/2024_U5_5_fig2_vif.png           ... VIF棒グラフ
  html/figures/2024_U5_5_fig3_coef.png          ... 男女別重回帰係数比較
  html/figures/2024_U5_5_fig4_scatter.png       ... ウォーキング率 vs スポーツ参加率散布図
=================================================================================
"""

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


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

from scipy import stats as scipy_stats
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor

# ── パス設定 ─────────────────────────────────────────────────────────────────
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,
})

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step1. 実データ読み込み
# ═══════════════════════════════════════════════════════════════════════════════
print("=" * 60)
print("■ 実データ読み込み（SSDSE-D-2023 / SSDSE-B-2026 / SSDSE-E-2026）")
print("=" * 60)

# SSDSE-D（社会生活基本調査）
df_d = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-D-2023.csv'),
                    encoding='cp932', header=1)
df_d = df_d[
    (df_d['男女の別'].isin(['0_総数', '1_男', '2_女'])) &
    (df_d['地域コード'] != 'R00000')
].copy()

# 男女別に分割
df_d_m = df_d[df_d['男女の別'] == '1_男'].set_index('都道府県').copy()
df_d_f = df_d[df_d['男女の別'] == '2_女'].set_index('都道府県').copy()
df_d_t = df_d[df_d['男女の別'] == '0_総数'].set_index('都道府県').copy()

# スポーツ・睡眠変数を数値化
for df_tmp in [df_d_m, df_d_f, df_d_t]:
    for c in ['スポーツの総数', 'ウォーキング・軽い体操', '睡眠']:
        df_tmp[c] = pd.to_numeric(df_tmp[c], errors='coerce')

print(f"SSDSE-D: 男={len(df_d_m)}件, 女={len(df_d_f)}件")

# 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')

df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口']
df_b['保育所千対'] = df_b['保育所等数'] / df_b['総人口'] * 10000

# 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('都道府県')
df_e['県民所得'] = pd.to_numeric(df_e['1人当たり県民所得（平成27年基準）'], errors='coerce')

# ─ データ統合 ─
common_prefs = (set(df_d_m.index) & set(df_d_f.index) &
                set(df_b.index) & set(df_e.index))
prefs = sorted(list(common_prefs))

sport_m  = df_d_m.loc[prefs, 'スポーツの総数'].values.astype(float)
sport_f  = df_d_f.loc[prefs, 'スポーツの総数'].values.astype(float)
walk_t   = df_d_t.loc[prefs, 'ウォーキング・軽い体操'].values.astype(float)
sleep_t  = df_d_t.loc[prefs, '睡眠'].values.astype(float)
nursery  = df_b.loc[prefs, '保育所千対'].values.astype(float)
temp     = df_b.loc[prefs, '年平均気温'].values.astype(float)
aging    = df_b.loc[prefs, '高齢化率'].values.astype(float)
income   = df_e.loc[prefs, '県民所得'].values.astype(float)

N = len(prefs)
VAR_NAMES = ['ウォーキング率', '睡眠時間', '保育所千対', '年平均気温', '県民所得', '高齢化率']
X_vars = np.column_stack([walk_t, sleep_t, nursery, temp, income, aging])

# 欠損チェック
valid_mask = ~np.any(np.isnan(X_vars), axis=1) & ~np.isnan(sport_m) & ~np.isnan(sport_f)
X_vars = X_vars[valid_mask]
sport_m = sport_m[valid_mask]
sport_f = sport_f[valid_mask]
prefs_v  = [prefs[i] for i in range(N) if valid_mask[i]]
N = len(prefs_v)

print(f"\n分析対象: {N}都道府県")
print(f"スポーツ参加率（男）: mean={sport_m.mean():.1f}%, std={sport_m.std():.1f}%")
print(f"スポーツ参加率（女）: mean={sport_f.mean():.1f}%, std={sport_f.std():.1f}%")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step2. 相関分析
# ═══════════════════════════════════════════════════════════════════════════════
df_corr = pd.DataFrame(X_vars, columns=VAR_NAMES)
df_corr['スポーツ参加率（男）'] = sport_m
df_corr['スポーツ参加率（女）'] = sport_f
corr_matrix = df_corr.corr()

print("\n【相関行列（スポーツ参加率との相関）】")
print(corr_matrix[['スポーツ参加率（男）', 'スポーツ参加率（女）']].round(3))

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step3. VIF計算
# ═══════════════════════════════════════════════════════════════════════════════
X_for_vif = sm.add_constant(X_vars)
vif_values = [variance_inflation_factor(X_for_vif, i + 1) for i in range(len(VAR_NAMES))]
vif_df = pd.DataFrame({'変数': VAR_NAMES, 'VIF': vif_values})
print("\n【VIF（分散拡大係数）】")
print(vif_df.round(3))

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step4. 重回帰分析（男女別）
# ═══════════════════════════════════════════════════════════════════════════════
X_reg = sm.add_constant(X_vars)
ols_m = sm.OLS(sport_m, X_reg).fit()
ols_f = sm.OLS(sport_f, X_reg).fit()

print("\n【重回帰結果 男性】")
print(f"  R² = {ols_m.rsquared:.3f}, AdjR² = {ols_m.rsquared_adj:.3f}")
print(ols_m.summary().tables[1])

print("\n【重回帰結果 女性】")
print(f"  R² = {ols_f.rsquared:.3f}, AdjR² = {ols_f.rsquared_adj:.3f}")
print(ols_f.summary().tables[1])

coef_m = ols_m.params[1:]
coef_f = ols_f.params[1:]
pval_m = ols_m.pvalues[1:]
pval_f = ols_f.pvalues[1:]

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

# ────────────────────────────────────────────────────────────────────────────
# 図1：相関係数ヒートマップ
# ────────────────────────────────────────────────────────────────────────────
print("\n図1: 相関係数ヒートマップを作成中...")

fig1, ax1 = plt.subplots(figsize=(9, 8))
corr_vals = corr_matrix.values
n_vars_all = len(corr_matrix.columns)
im = ax1.imshow(corr_vals, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax1, label='Pearson相関係数')
ax1.set_xticks(range(n_vars_all))
ax1.set_yticks(range(n_vars_all))
ax1.set_xticklabels(corr_matrix.columns, fontsize=8.5, rotation=25, ha='right')
ax1.set_yticklabels(corr_matrix.columns, fontsize=8.5)
ax1.set_title('Pearson相関係数行列\n（スポーツ参加率 × 地域環境変数）\nデータ：SSDSE実データ',
              fontsize=12, fontweight='bold')
for i in range(n_vars_all):
    for j in range(n_vars_all):
        val = corr_vals[i, j]
        sig = ''
        if i != j:
            try:
                _, p = scipy_stats.pearsonr(df_corr.iloc[:, i].dropna(),
                                             df_corr.iloc[:, j].dropna())
                sig = '*' if p < 0.05 else ''
            except Exception:
                pass
        text_color = 'white' if abs(val) > 0.6 else 'black'
        ax1.text(j, i, f'{val:.2f}{sig}', ha='center', va='center',
                 fontsize=8, fontweight='bold', color=text_color)

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

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

fig2, ax2 = plt.subplots(figsize=(8, 5))
colors2 = ['#C62828' if v > 5 else '#FB8C00' if v > 2.5 else '#1565C0' for v in vif_values]
bars2 = ax2.bar(VAR_NAMES, vif_values, color=colors2, edgecolor='white', alpha=0.88)
ax2.axhline(5, color='#C62828', linestyle='--', linewidth=1.5, label='VIF=5（問題の目安）')
ax2.axhline(2.5, color='#FB8C00', linestyle=':', linewidth=1.5, label='VIF=2.5（注意）')
ax2.set_xticklabels(VAR_NAMES, fontsize=10, rotation=15, ha='right')
ax2.set_ylabel('VIF（分散拡大係数）', fontsize=12)
ax2.set_title('VIF（多重共線性の確認）\nVIF<5 → 多重共線性の問題なし', fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(axis='y', alpha=0.3)
for bar, val in zip(bars2, vif_values):
    ax2.text(bar.get_x() + bar.get_width()/2, val + 0.05,
             f'{val:.2f}', ha='center', va='bottom', fontsize=10, fontweight='bold')

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

# ────────────────────────────────────────────────────────────────────────────
# 図3：男女別重回帰係数比較
# ────────────────────────────────────────────────────────────────────────────
print("図3: 男女別係数比較を作成中...")

x = np.arange(len(VAR_NAMES))
width = 0.35

fig3, ax3 = plt.subplots(figsize=(11, 6))
bars_m = ax3.bar(x - width/2, coef_m, width, label='男性', color='#1565C0', alpha=0.82, edgecolor='white')
bars_f = ax3.bar(x + width/2, coef_f, width, label='女性', color='#E91E63', alpha=0.82, edgecolor='white')
ax3.axhline(0, color='black', linewidth=0.8)
ax3.set_xticks(x)
ax3.set_xticklabels(VAR_NAMES, fontsize=10, rotation=10, ha='right')
ax3.set_ylabel('回帰係数（非標準化）', fontsize=12)
ax3.set_title(f'男女別 重回帰係数の比較\n（*: p<0.05）　男性R²={ols_m.rsquared:.3f}  女性R²={ols_f.rsquared:.3f}',
              fontsize=12, fontweight='bold')
ax3.legend(fontsize=11)
ax3.grid(axis='y', alpha=0.3)

for i, (bar_m, bar_f, pm, pf) in enumerate(zip(bars_m, bars_f, pval_m, pval_f)):
    hm = bar_m.get_height()
    hf = bar_f.get_height()
    if pm < 0.05:
        ax3.text(bar_m.get_x() + bar_m.get_width()/2,
                 hm + (0.02 if hm >= 0 else -0.05), '*',
                 ha='center', fontsize=14, color='#1565C0', fontweight='bold')
    if pf < 0.05:
        ax3.text(bar_f.get_x() + bar_f.get_width()/2,
                 hf + (0.02 if hf >= 0 else -0.05), '*',
                 ha='center', fontsize=14, color='#E91E63', fontweight='bold')

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

# ────────────────────────────────────────────────────────────────────────────
# 図4：ウォーキング率 vs スポーツ参加率 散布図（男女別）
# ────────────────────────────────────────────────────────────────────────────
print("図4: ウォーキング率 vs スポーツ参加率散布図を作成中...")

walk_vals = X_vars[:, 0]  # ウォーキング率
r_m, p_m = scipy_stats.pearsonr(walk_vals, sport_m)
r_f, p_f = scipy_stats.pearsonr(walk_vals, sport_f)

fig4, axes4 = plt.subplots(1, 2, figsize=(13, 5))
fig4.suptitle('ウォーキング率とスポーツ参加率の関係（都道府県別実データ）', fontsize=13, fontweight='bold')

for ax, score, gender, clr, r_v, p_v in zip(
    axes4,
    [sport_m, sport_f],
    ['男性', '女性'],
    ['#1565C0', '#E91E63'],
    [r_m, r_f],
    [p_m, p_f],
):
    ax.scatter(walk_vals, score, c=clr, alpha=0.75, s=55, edgecolors='white', linewidth=0.5)
    coef_fit = np.polyfit(walk_vals, score, 1)
    x_fit = np.linspace(walk_vals.min(), walk_vals.max(), 100)
    ax.plot(x_fit, np.polyval(coef_fit, x_fit), color='#333', linewidth=2, linestyle='--')
    # 代表的な都道府県ラベル
    for i, pref in enumerate(prefs_v):
        if pref in ['東京都', '秋田県', '北海道', '沖縄県', '長野県']:
            ax.annotate(pref.replace('県','').replace('府','').replace('都','').replace('道',''),
                        (walk_vals[i], score[i]),
                        textcoords='offset points', xytext=(5, 3), fontsize=8)
    ax.set_xlabel('ウォーキング・軽い体操 行動者率（%）', fontsize=11)
    ax.set_ylabel('スポーツ参加率（%）', fontsize=11)
    ax.set_title(f'{gender}（r = {r_v:.3f}, {"*" if p_v < 0.05 else "ns"}）', fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.2)
    ax.text(0.05, 0.95, f'r = {r_v:.3f}\n回帰係数 = {coef_fit[0]:.3f}', transform=ax.transAxes,
            fontsize=10, va='top', bbox=dict(boxstyle='round', facecolor='#E3F2FD', alpha=0.8))

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

print("\n" + "=" * 60)
print("✓ 全図の生成完了（4枚）")
print("=" * 60)
print("\n【主要知見】")
print(f"  男性モデル R² = {ols_m.rsquared:.3f}")
print(f"  女性モデル R² = {ols_f.rsquared:.3f}")
print(f"  ウォーキング率との相関: 男性 r={r_m:.3f}, 女性 r={r_f:.3f}")
print(f"  VIF最大値: {max(vif_values):.2f}（多重共線性の問題{'あり' if max(vif_values)>5 else 'なし'}）")
print(f"  男性に有意な変数: {[VAR_NAMES[i] for i, p in enumerate(pval_m) if p < 0.05]}")
print(f"  女性に有意な変数: {[VAR_NAMES[i] for i, p in enumerate(pval_f) if p < 0.05]}")
print(f"  使用データ: SSDSE-D-2023, SSDSE-B-2026, SSDSE-E-2026 ({N}都道府県)")