"""
教育用再現コード: 2023年度 統計データ分析コンペティション 審査員奨励賞（大学生）
=================================================================================
論文タイトル：英語教育実施状況調査をエビデンス生成に繋げるための探究
受賞区分  ：審査員奨励賞 [大学生の部]

【分析概要】
  データ：SSDSE-B-2026.csv（47都道府県 × 2022年度）
  目的   ：英語教育の質・成果を直接測定する統計データは公的統計に存在しないため、
            都道府県レベルの教育・経済・人口指標を代理変数として用い、
            「大学進学率」を英語能力（学習達成度）の代理指標として設定し、
            その地域差を規定する要因を探索的に分析する。

  Step1. SSDSE-B-2026.csv から実データを読み込む（2022年度・47都道府県）
  Step2. 各変数のエンジニアリングと記述統計
  Step3. Pearson相関分析：大学進学率と各変数の相関係数
  Step4. 重回帰分析：大学進学率 ~ 教育費割合 + 住宅地価格 + 教員比
  Step5. 地理的比較：47都道府県の大学進学率ランキング
  Step6. 都市農村格差：住宅地価格グループ別の比較

【変数定義】
  大学進学率（英語能力 proxy）: E4602 / E4601 × 100
  教育費割合（教育投資）      : L322108 / L3221 × 100
  住宅地価格（経済格差）      : C5401（千円/m²）
  小学教員/児童比（学習環境） : E2401 / E2501
  交通通信費割合（ICT環境）   : L322107 / L3221
  大学生/15-64歳人口比        : E6302 / A1302 × 100
  高齢化率                    : A1303 / A1101 × 100
  国際旅券発行数per capita    : G5105 / A1101 × 10000

【図の出力】
  html/figures/2023_U5_3_fig1_corr.png    ... 大学進学率との相関係数棒グラフ
  html/figures/2023_U5_3_fig2_scatter.png ... 教育費割合 vs 大学進学率（回帰線付き）
  html/figures/2023_U5_3_fig3_ranking.png ... 47都道府県の大学進学率水平棒グラフ
  html/figures/2023_U5_3_fig4_box.png     ... 住宅地価格グループ別 大学進学率箱ひげ図
=================================================================================
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・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/2023_U5_3_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
import statsmodels.api as sm

# ── パス設定 ─────────────────────────────────────────────────────────────────
BASE_DIR  = os.path.join(_script_dir, '..')
FIG_DIR   = os.path.join(BASE_DIR, 'html', 'figures')
DATA_PATH = os.path.join(BASE_DIR, 'data', 'raw', 'SSDSE-B-2026.csv')
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. 実データ読み込み（SSDSE-B-2026, 2022年度・47都道府県）
# ═══════════════════════════════════════════════════════════════════════════════
print("=" * 60)
print("■ SSDSE-B-2026.csv 読み込み（2022年度）")

df_all = pd.read_csv(DATA_PATH, encoding='cp932', header=0, skiprows=[1])
df = df_all[df_all['SSDSE-B-2026'] == 2022].copy().reset_index(drop=True)

assert len(df) == 47, f"47都道府県分のデータが必要ですが {len(df)} 行しかありません"
print(f"  読み込み完了: {len(df)} 都道府県")

PREFS = list(df['Prefecture'])

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step2. 変数エンジニアリング
# ═══════════════════════════════════════════════════════════════════════════════
print("\n■ 変数エンジニアリング")

# 目的変数（英語能力 proxy）
y_univ = df['E4602'] / df['E4601'] * 100        # 大学進学率 (%)

# 説明変数群
edu_cost   = df['L322108'] / df['L3221'] * 100  # 教育費割合 (%)
land_price = df['C5401'].astype(float)           # 住宅地価格（千円/m²）
teacher_ratio = df['E2401'] / df['E2501']        # 小学教員/児童比
ict_share  = df['L322107'] / df['L3221']         # 交通通信費割合
univ_ratio = df['E6302'] / df['A1302'] * 100    # 大学生/15-64歳人口比 (%)
aging_rate = df['A1303'] / df['A1101'] * 100    # 高齢化率 (%)
passport_pc = df['G5105'] / df['A1101'] * 10000 # 国際旅券発行数per万人

feat = pd.DataFrame({
    '教育費割合(%)':          edu_cost,
    '住宅地価格(千円/m²)':    land_price,
    '小学教員/児童比':        teacher_ratio,
    '交通通信費割合':         ict_share,
    '大学生/15-64歳比(%)':    univ_ratio,
    '高齢化率(%)':            aging_rate,
    '旅券発行per万人':        passport_pc,
}, index=df.index)

FEAT_NAMES = list(feat.columns)
print(f"  変数数: {len(FEAT_NAMES)}")
print(f"  目的変数 大学進学率: 平均={y_univ.mean():.1f}%, 範囲=[{y_univ.min():.1f}, {y_univ.max():.1f}]%")

# ── 記述統計 ─────────────────────────────────────────────────────────────────
print("\n■ 記述統計")
desc = feat.describe().T[['mean','std','min','max']]
desc.insert(0, '大学進学率', y_univ.describe()[['mean','std','min','max']].iloc[0] if False else None)
print(pd.concat([
    pd.DataFrame({'mean': [y_univ.mean()], 'std': [y_univ.std()],
                  'min':  [y_univ.min()],  'max': [y_univ.max()]},
                 index=['大学進学率(%)']),
    desc
]).round(3).to_string())

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step3. Pearson 相関分析：大学進学率 vs 各変数
# ═══════════════════════════════════════════════════════════════════════════════
print("\n■ Pearson 相関分析（大学進学率 vs 各変数）")

corr_results = []
y_arr = y_univ.values.astype(float)
for col in FEAT_NAMES:
    x_arr = feat[col].values.astype(float)
    r, p  = stats.pearsonr(x_arr, y_arr)
    corr_results.append({'変数': col, '相関係数 r': r, 'p値': p,
                         '有意 (p<0.05)': '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else 'n.s.'))})

corr_df = pd.DataFrame(corr_results).sort_values('相関係数 r', ascending=False)
print(corr_df.round(4).to_string(index=False))

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step4. 重回帰分析：大学進学率 ~ 教育費割合 + 住宅地価格 + 教員比
# ═══════════════════════════════════════════════════════════════════════════════
print("\n■ 重回帰分析（OLS）")

X_reg = sm.add_constant(feat[['教育費割合(%)', '住宅地価格(千円/m²)', '小学教員/児童比']].astype(float))
ols_model = sm.OLS(y_arr, X_reg).fit()
print(ols_model.summary())

# ── 教育費割合 vs 大学進学率 の単回帰（scatter 用）─────────────────────────
x_edu = edu_cost.values.astype(float)
slope, intercept, r_val, p_val, se = stats.linregress(x_edu, y_arr)
print(f"\n  教育費割合 vs 大学進学率: r={r_val:.3f}, p={p_val:.4f}, slope={slope:.3f}")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step5. 都市農村格差：住宅地価格によるグループ分け（3群）
# ═══════════════════════════════════════════════════════════════════════════════
print("\n■ 都市農村格差（住宅地価格 3群）")

tertile_labels = ['低価格帯\n（農村型）', '中価格帯\n（中間型）', '高価格帯\n（都市型）']
land_arr = land_price.values.astype(float)
t33 = np.percentile(land_arr, 33.33)
t67 = np.percentile(land_arr, 66.67)
groups = []
for v in land_arr:
    if v <= t33:
        groups.append(0)
    elif v <= t67:
        groups.append(1)
    else:
        groups.append(2)
groups = np.array(groups)

for i, label in enumerate(tertile_labels):
    subset = y_arr[groups == i]
    print(f"  {label.replace(chr(10),'')}: N={len(subset)}, 平均={subset.mean():.1f}%, SD={subset.std():.1f}%")

# ANOVA
f_stat, anova_p = stats.f_oneway(y_arr[groups == 0], y_arr[groups == 1], y_arr[groups == 2])
print(f"  一元配置ANOVA: F={f_stat:.3f}, p={anova_p:.4f}")

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

# ────────────────────────────────────────────────────────────────────────────
# 図1：各変数と大学進学率の相関係数棒グラフ
# ────────────────────────────────────────────────────────────────────────────
print("\n図1: 相関係数棒グラフを作成中...")

corr_sorted = corr_df.sort_values('相関係数 r')
colors1 = []
for _, row in corr_sorted.iterrows():
    if row['有意 (p<0.05)'] != 'n.s.':
        colors1.append('#C62828' if row['相関係数 r'] >= 0 else '#1565C0')
    else:
        colors1.append('#BDBDBD')

fig1, ax1 = plt.subplots(figsize=(9, 5))
bars = ax1.barh(corr_sorted['変数'], corr_sorted['相関係数 r'],
                color=colors1, edgecolor='white', alpha=0.88, height=0.6)
ax1.axvline(0, color='black', linewidth=1.0)
ax1.set_xlabel('Pearson 相関係数 r', fontsize=12)
ax1.set_title('各変数と大学進学率（英語能力 proxy）との相関係数\n（赤：正・有意,  青：負・有意,  灰：非有意）',
              fontsize=12, fontweight='bold')
ax1.grid(axis='x', alpha=0.3)
ax1.set_xlim(-1.05, 1.05)

for bar, (_, row) in zip(ax1.patches, corr_sorted.iterrows()):
    r = row['相関係数 r']
    sig = row['有意 (p<0.05)']
    x_pos = r + 0.03 * np.sign(r) if r != 0 else 0.03
    ha = 'left' if r >= 0 else 'right'
    label = f'{r:.3f}{sig if sig != "n.s." else ""}'
    ax1.text(x_pos, bar.get_y() + bar.get_height() / 2, label,
             va='center', ha=ha, fontsize=9, color='#222222')

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2023_U5_3_fig1_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("  -> 2023_U5_3_fig1_corr.png 保存完了")

# ────────────────────────────────────────────────────────────────────────────
# 図2：教育費割合 vs 大学進学率（回帰線 + 都道府県ラベル）
# ────────────────────────────────────────────────────────────────────────────
print("図2: 散布図（教育費割合 vs 大学進学率）を作成中...")

x_plot = np.linspace(x_edu.min() - 0.1, x_edu.max() + 0.1, 200)
y_fit  = intercept + slope * x_plot

fig2, ax2 = plt.subplots(figsize=(10, 7))
ax2.scatter(x_edu, y_arr, color='#1565C0', s=55, alpha=0.75, zorder=3,
            edgecolors='white', linewidths=0.5)
ax2.plot(x_plot, y_fit, color='#C62828', linewidth=2.0, linestyle='-',
         label=f'回帰線 (r={r_val:.3f}, p={p_val:.4f})')

# 都道府県ラベル（全47都道府県）
for i, pref in enumerate(PREFS):
    short = pref.replace('県', '').replace('府', '').replace('都', '').replace('道', '')
    ax2.text(x_edu[i] + 0.01, y_arr[i] + 0.3, short,
             fontsize=6.5, color='#444444', va='bottom', ha='left')

ax2.set_xlabel('教育費割合 (%, 消費支出に占める教育費)', fontsize=12)
ax2.set_ylabel('大学進学率 (%, 英語能力 proxy)', fontsize=12)
ax2.set_title('教育費割合と大学進学率の関係（47都道府県, 2022年度）\n'
              'データ出所: SSDSE-B-2026（統計数理研究所）',
              fontsize=12, fontweight='bold')
ax2.legend(fontsize=11)
ax2.grid(True, alpha=0.2)
ax2.text(0.97, 0.05,
         f'y = {slope:.2f}x + {intercept:.1f}\nr = {r_val:.3f}  p = {p_val:.4f}',
         transform=ax2.transAxes, fontsize=10, va='bottom', ha='right',
         bbox=dict(boxstyle='round', facecolor='#FFF9C4', alpha=0.85))

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2023_U5_3_fig2_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("  -> 2023_U5_3_fig2_scatter.png 保存完了")

# ────────────────────────────────────────────────────────────────────────────
# 図3：47都道府県の大学進学率 水平棒グラフ（ランキング順）
# ────────────────────────────────────────────────────────────────────────────
print("図3: 47都道府県 大学進学率ランキングを作成中...")

rank_df = pd.DataFrame({'都道府県': PREFS, '大学進学率': y_arr}).sort_values('大学進学率')

# 上位5・下位5を強調
top5_set  = set(rank_df['都道府県'].iloc[-5:])
bot5_set  = set(rank_df['都道府県'].iloc[:5])
avg_val   = y_arr.mean()

colors3 = []
for pref in rank_df['都道府県']:
    if pref in top5_set:
        colors3.append('#C62828')
    elif pref in bot5_set:
        colors3.append('#1565C0')
    else:
        colors3.append('#90A4AE')

fig3, ax3 = plt.subplots(figsize=(9, 14))
bars3 = ax3.barh(rank_df['都道府県'], rank_df['大学進学率'],
                 color=colors3, edgecolor='white', alpha=0.88, height=0.75)
ax3.axvline(avg_val, color='black', linestyle='--', linewidth=1.5,
            label=f'全国平均 {avg_val:.1f}%')
ax3.set_xlabel('大学進学率 (%)', fontsize=12)
ax3.set_title('47都道府県 大学進学率ランキング（2022年度）\n'
              '赤：上位5県, 青：下位5県',
              fontsize=12, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(axis='x', alpha=0.3)
ax3.set_xlim(0, rank_df['大学進学率'].max() * 1.12)

for bar, val in zip(ax3.patches, rank_df['大学進学率']):
    ax3.text(val + 0.3, bar.get_y() + bar.get_height() / 2,
             f'{val:.1f}%', va='center', ha='left', fontsize=7.5, color='#333333')

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2023_U5_3_fig3_ranking.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("  -> 2023_U5_3_fig3_ranking.png 保存完了")

# ────────────────────────────────────────────────────────────────────────────
# 図4：住宅地価格グループ別 大学進学率の箱ひげ図
# ────────────────────────────────────────────────────────────────────────────
print("図4: 箱ひげ図（住宅地価格グループ別）を作成中...")

group_data = [y_arr[groups == i] for i in range(3)]
group_colors = ['#42A5F5', '#66BB6A', '#EF5350']

fig4, ax4 = plt.subplots(figsize=(8, 6))
bp = ax4.boxplot(group_data, patch_artist=True, notch=False,
                 medianprops=dict(color='black', linewidth=2),
                 whiskerprops=dict(linewidth=1.5),
                 capprops=dict(linewidth=1.5),
                 flierprops=dict(marker='o', markersize=5, markerfacecolor='gray', alpha=0.6))

for patch, color in zip(bp['boxes'], group_colors):
    patch.set_facecolor(color)
    patch.set_alpha(0.65)

# 個別データ点を重ねる（jitter）
for i, gd in enumerate(group_data):
    jitter = (np.arange(len(gd)) - len(gd) / 2) * 0.05
    ax4.scatter(np.full(len(gd), i + 1) + jitter, gd,
                color=group_colors[i], s=30, alpha=0.7, zorder=3, edgecolors='white', linewidths=0.5)

ax4.set_xticklabels(tertile_labels, fontsize=11)
ax4.set_ylabel('大学進学率 (%)', fontsize=12)
ax4.set_title(f'住宅地価格グループ別 大学進学率の分布\n'
              f'（一元配置ANOVA: F={f_stat:.2f}, p={anova_p:.4f}）',
              fontsize=12, fontweight='bold')
ax4.grid(axis='y', alpha=0.3)

# グループ平均ラベル
for i, gd in enumerate(group_data):
    ax4.text(i + 1, gd.max() + 0.5, f'平均\n{gd.mean():.1f}%',
             ha='center', va='bottom', fontsize=9.5, fontweight='bold',
             color=group_colors[i])

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2023_U5_3_fig4_box.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("  -> 2023_U5_3_fig4_box.png 保存完了")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ 主要知見サマリー
# ═══════════════════════════════════════════════════════════════════════════════
print("\n" + "=" * 60)
print("✓ 全図の生成完了（4枚）")
print("=" * 60)
print("\n【主要知見】")
print(f"  データ         : SSDSE-B-2026.csv（2022年度・47都道府県）")
print(f"  目的変数       : 大学進学率（英語能力 proxy）")
print(f"                   全国平均={y_arr.mean():.1f}%, 最大={y_arr.max():.1f}%, 最小={y_arr.min():.1f}%")
print(f"  相関 Top変数   : {corr_df.iloc[0]['変数']} (r={corr_df.iloc[0]['相関係数 r']:.3f})")
print(f"  重回帰 R²      : {ols_model.rsquared:.3f}")
print(f"  都市農村格差   : 低価格帯{y_arr[groups==0].mean():.1f}% vs 高価格帯{y_arr[groups==2].mean():.1f}%")
print(f"  ANOVA          : F={f_stat:.3f}, p={anova_p:.4f}")
top3 = pd.Series(y_arr, index=PREFS).nlargest(3)
bot3 = pd.Series(y_arr, index=PREFS).nsmallest(3)
print(f"  大学進学率 上位3: {list(top3.index)}")
print(f"  大学進学率 下位3: {list(bot3.index)}")