"""
教育用再現コード: 2025年 統計データ分析コンペティション 審査員奨励賞 [高校生の部]
=================================================================
論文タイトル：5つの視点に基づく不登校の原因究明と対策
著者：早稲田大学系属 早稲田実業学校

【データ出典】
  従属変数（不登校率）:
    文部科学省「児童生徒の問題行動・不登校等生徒指導上の諸課題に関する調査」
    平成29年度（2017年度） 表4-14 都道府県別 不登校児童生徒数
    出典: https://www.e-stat.go.jp/ 政府統計の総合窓口

  教育支援センター設置数:
    文部科学省 上記調査 表4-15 都道府県別教育支援センター状況（2017年度）

  説明変数:
    SSDSE-D-2023.csv（社会生活基本調査 2021年度）: 睡眠時間, 通学時間, 仕事時間
    SSDSE-B-2026.csv（2017年度）: 出生率, 教育費（二人以上世帯）
    SSDSE-E-2026.csv: 1人当たり県民所得
    SSDSE-F-2023v3.csv: 年間日照時間

【分析概要】
  5カテゴリの説明変数で都道府県別不登校率（小中計, 2017年度）を解析
    a.子ども: 睡眠時間, 趣味・娯楽時間（インターネット利用の代替指標）
    b.親: 父親仕事時間, 合計特殊出生率
    c.学校: 通学時間, 教育費
    d.教育支援機関: 教育支援センター数（人口万対）
    e.その他: 1人当たり県民所得, 年間日照時間

  Step1. カテゴリ別相関分析
  Step2. 重回帰分析（全変数）
  Step3. 外れ値分析（IQR法）
=================================================================
"""

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


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy import stats
from matplotlib.patches import Patch
import warnings
warnings.filterwarnings('ignore')
import os

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

BASE_DIR = os.path.join(_script_dir, '..')
FIG_DIR  = os.path.join(BASE_DIR, 'html', 'figures')
DATA_DIR = os.path.join(BASE_DIR, 'data', 'raw')

# ================================================================
# ■ 実データ埋め込み（文部科学省 平成29年度調査）
#   出典: 問題行動・不登校等調査 表4-14, 4-15 都道府県別
# ================================================================
# 不登校率（小中学校計, 1,000人当たり） — 2017年度 MEXT調査より
FUTOKO_RATE = {
    '北海道': 14.9, '青森県': 13.9, '岩手県': 11.2, '宮城県': 19.1,
    '秋田県': 10.8, '山形県': 12.1, '福島県': 13.2, '茨城県': 14.8,
    '栃木県': 16.8, '群馬県': 14.1, '埼玉県': 11.8, '千葉県': 13.3,
    '東京都': 14.5, '神奈川県': 17.6, '新潟県': 13.7, '富山県': 11.4,
    '石川県': 15.3, '福井県': 11.7, '山梨県': 15.2, '長野県': 15.3,
    '岐阜県': 15.4, '静岡県': 17.4, '愛知県': 16.7, '三重県': 14.9,
    '滋賀県': 13.8, '京都府': 13.7, '大阪府': 16.0, '兵庫県': 15.3,
    '奈良県': 13.0, '和歌山県': 13.4, '鳥取県': 14.4, '島根県': 16.8,
    '岡山県': 13.0, '広島県': 13.3, '山口県': 12.6, '徳島県': 11.5,
    '香川県': 13.4, '愛媛県': 11.4, '高知県': 17.7, '福岡県': 13.6,
    '佐賀県': 14.3, '長崎県': 13.3, '熊本県': 13.2, '大分県': 15.0,
    '宮崎県': 12.0, '鹿児島県': 12.4, '沖縄県': 17.3
}

# 教育支援センター設置数 — 2017年度 MEXT調査 表4-15より
SHIEN_CENTER = {
    '北海道': 55, '青森県': 12, '岩手県': 23, '宮城県': 33, '秋田県': 14,
    '山形県': 23, '福島県': 22, '茨城県': 50, '栃木県': 28, '群馬県': 35,
    '埼玉県': 66, '千葉県': 58, '東京都': 78, '神奈川県': 63, '新潟県': 39,
    '富山県': 16, '石川県': 19, '福井県': 19, '山梨県': 14, '長野県': 63,
    '岐阜県': 40, '静岡県': 41, '愛知県': 67, '三重県': 21, '滋賀県': 26,
    '京都府': 20, '大阪府': 39, '兵庫県': 53, '奈良県': 12, '和歌山県': 15,
    '鳥取県': 14, '島根県': 13, '岡山県': 27, '広島県': 28, '山口県': 20,
    '徳島県': 12, '香川県': 17, '愛媛県': 15, '高知県': 24, '福岡県': 42,
    '佐賀県': 19, '長崎県': 13, '熊本県': 30, '大分県': 20, '宮崎県': 21,
    '鹿児島県': 26, '沖縄県': 16
}

PREFS = list(FUTOKO_RATE.keys())  # 47都道府県（順序確定）

# ================================================================
# ■ 説明変数の読み込み
# ================================================================
# SSDSE-D: 社会生活基本調査 2021年（時間はすべて分/日）
df_d_all = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-D-2023.csv'), encoding='cp932', header=1)
df_d_total = df_d_all[(df_d_all['男女の別'] == '0_総数') &
                      (df_d_all['地域コード'] != 'R00000')].copy()
df_d_male  = df_d_all[(df_d_all['男女の別'] == '1_男') &
                      (df_d_all['地域コード'] != 'R00000')].copy()

sleep_map    = dict(zip(df_d_total['都道府県'], pd.to_numeric(df_d_total['睡眠'],             errors='coerce')))
leisure_map  = dict(zip(df_d_total['都道府県'], pd.to_numeric(df_d_total['趣味・娯楽'],        errors='coerce')))
commute_map  = dict(zip(df_d_total['都道府県'], pd.to_numeric(df_d_total['通勤・通学'],         errors='coerce')))
work_m_map   = dict(zip(df_d_male['都道府県'],  pd.to_numeric(df_d_male['仕事'],               errors='coerce')))

# SSDSE-B 2017: 出生率, 教育費
ssdse_b = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'), encoding='cp932', header=1)
df_b17 = ssdse_b[
    (ssdse_b['年度'] == 2017) &
    ssdse_b['地域コード'].str.match(r'^R\d{5}$', na=False)
].reset_index(drop=True)
tfr_map  = dict(zip(df_b17['都道府県'], pd.to_numeric(df_b17['合計特殊出生率'],          errors='coerce')))
edu_map  = dict(zip(df_b17['都道府県'], pd.to_numeric(df_b17['教育費（二人以上の世帯）'], errors='coerce')))
pop_map  = dict(zip(df_b17['都道府県'], pd.to_numeric(df_b17['総人口'],                  errors='coerce')))

# SSDSE-E: 1人当たり県民所得
ssdse_e_raw = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-E-2026.csv'), encoding='cp932', header=1)
ssdse_e = ssdse_e_raw.iloc[1:].copy()
ssdse_e.columns = ssdse_e_raw.iloc[0].values
ssdse_e = ssdse_e[ssdse_e['都道府県'] != '全国'].reset_index(drop=True)
income_map = dict(zip(
    ssdse_e['都道府県'],
    pd.to_numeric(ssdse_e['1人当たり県民所得（平成27年基準）'], errors='coerce')
))

# SSDSE-F: 年間日照時間（都道府県ごとに観測地点平均）
df_f = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-F-2023v3.csv'), encoding='cp932', header=1)
annual_f = (df_f.groupby(['都道府県', '市'])['日照時間の合計']
            .sum().reset_index()
            .groupby('都道府県')['日照時間の合計'].mean())
sunshine_map = annual_f.to_dict()

# ================================================================
# ■ データフレームの構築
# ================================================================
df = pd.DataFrame({'都道府県': PREFS})
df['不登校率']         = df['都道府県'].map(FUTOKO_RATE)
df['睡眠時間']         = df['都道府県'].map(sleep_map)
df['趣味娯楽時間']     = df['都道府県'].map(leisure_map)   # インターネット利用の代替指標
df['父親仕事時間']     = df['都道府県'].map(work_m_map)
df['出生率']           = df['都道府県'].map(tfr_map)
df['通学時間']         = df['都道府県'].map(commute_map)
df['教育費']           = df['都道府県'].map(edu_map)        # 円/月
df['支援センター万対'] = df['都道府県'].map(SHIEN_CENTER) / df['都道府県'].map(pop_map) * 10000
df['県民所得']         = df['都道府県'].map(income_map)     # 万円
df['日照時間']         = df['都道府県'].map(sunshine_map)   # 時間/年

df = df.dropna().reset_index(drop=True)
N = len(df)

VAR_NAMES_ALL = ['睡眠時間', '趣味娯楽時間', '父親仕事時間', '出生率',
                 '通学時間', '教育費', '支援センター万対', '県民所得', '日照時間']
CATEGORIES = {
    'a.子ども（生活習慣）': ['睡眠時間', '趣味娯楽時間'],
    'b.親（家庭環境）':     ['父親仕事時間', '出生率'],
    'c.学校（環境）':       ['通学時間', '教育費'],
    'd.教育支援機関':       ['支援センター万対'],
    'e.その他':             ['県民所得', '日照時間'],
}
CAT_COLORS = {
    'a.子ども（生活習慣）': '#C62828',
    'b.親（家庭環境）':     '#E65100',
    'c.学校（環境）':       '#1565C0',
    'd.教育支援機関':       '#2E7D32',
    'e.その他':             '#6A1B9A'
}

y = df['不登校率'].values
X_raw = df[VAR_NAMES_ALL].values

print("=" * 60)
print(f"■ 記述統計（47都道府県, 2017年度 MEXT調査）N = {N}")
print("=" * 60)
print(df[['不登校率'] + VAR_NAMES_ALL].describe().round(2))

# ================================================================
# ■ Step1. カテゴリ別相関分析
# ================================================================
print("\n" + "=" * 60)
print("■ Step1. カテゴリ別相関分析（目的変数：不登校率）")
print("=" * 60)
for cat, vlist in CATEGORIES.items():
    print(f"\n  【{cat}】")
    for vname in vlist:
        xv = df[vname].values
        r, p = stats.pearsonr(xv, y)
        sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
        print(f"    {vname:<16} r={r:>7.4f}  p={p:.4f}  {sig}")

# ================================================================
# ■ Step2. 重回帰分析
# ================================================================
print("\n" + "=" * 60)
print("■ Step2. 重回帰分析（全説明変数）")
print("=" * 60)
X_ols = sm.add_constant(X_raw)
model = sm.OLS(y, X_ols).fit()
print(model.summary2())

# ================================================================
# ■ Step3. IQR 外れ値の確認
# ================================================================
print("\n" + "=" * 60)
print("■ Step3. 外れ値（不登校率 IQR×1.5）")
print("=" * 60)
Q1, Q3 = np.percentile(y, 25), np.percentile(y, 75)
IQR_val = Q3 - Q1
upper, lower = Q3 + 1.5 * IQR_val, Q1 - 1.5 * IQR_val
outliers = df[(df['不登校率'] > upper) | (df['不登校率'] < lower)][['都道府県', '不登校率']]
print(f"  Q1={Q1:.2f}, Q3={Q3:.2f}, IQR={IQR_val:.2f}")
print(f"  上限={upper:.2f}, 下限={lower:.2f}")
print(f"  外れ値都道府県: {outliers['都道府県'].tolist()}")
print(outliers.to_string(index=False))

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

# ──────────────────────────────────────────────────────────
# 図1: 都道府県別不登校率 + IQR箱ひげ図
# ──────────────────────────────────────────────────────────
fig1, axes1 = plt.subplots(1, 2, figsize=(14, 6))
fig1.suptitle('都道府県別不登校率（小中学校計, 2017年度）\n出典: 文部科学省 問題行動等調査', fontsize=13, fontweight='bold')

# 棒グラフ（上位10）
sorted_idx = np.argsort(y)[::-1]
top10_prefs = df['都道府県'].values[sorted_idx[:10]]
top10_vals  = y[sorted_idx[:10]]
ax1a = axes1[0]
colors_bar = ['#C62828' if i < 3 else '#FF8F00' if i < 7 else '#1565C0' for i in range(10)]
ax1a.barh(range(10), top10_vals[::-1], color=colors_bar[::-1], alpha=0.8, edgecolor='white')
ax1a.set_yticks(range(10))
ax1a.set_yticklabels(top10_prefs[::-1], fontsize=10)
ax1a.set_xlabel('不登校率（小中計, 千対）', fontsize=11)
ax1a.set_title('不登校率 上位10都道府県', fontsize=11, fontweight='bold')
ax1a.axvline(y.mean(), color='gray', linestyle='--', linewidth=1.2, label=f'全国平均={y.mean():.2f}‰')
ax1a.legend(fontsize=9)
ax1a.grid(axis='x', alpha=0.3)

# IQR箱ひげ図
ax1b = axes1[1]
bp1 = ax1b.boxplot(y, patch_artist=True, vert=True,
                   boxprops=dict(facecolor='#BBDEFB', color='#1565C0'),
                   medianprops=dict(color='#C62828', linewidth=2),
                   flierprops=dict(marker='o', markerfacecolor='#FF8F00', markersize=7))
ax1b.axhline(upper, color='#C62828', linestyle='--', linewidth=1.2,
             label=f'IQR上限={upper:.2f}')
ax1b.axhline(lower, color='#1565C0', linestyle='--', linewidth=1.2,
             label=f'IQR下限={lower:.2f}')
for _, row in outliers.iterrows():
    ax1b.annotate(row['都道府県'], (1.02, row['不登校率']), fontsize=8)
ax1b.set_ylabel('不登校率（千対）', fontsize=11)
ax1b.set_title('不登校率の箱ひげ図（IQR×1.5）', fontsize=11, fontweight='bold')
ax1b.legend(fontsize=9)
ax1b.grid(axis='y', alpha=0.3)
ax1b.set_xlim(0.5, 1.5)

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2025_H5_4_fig1_pref.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2025_H5_4_fig1_pref.png")

# ──────────────────────────────────────────────────────────
# 図2: カテゴリ別相関係数の棒グラフ
# ──────────────────────────────────────────────────────────
fig2, ax2 = plt.subplots(figsize=(9, 5))
corrs = [stats.pearsonr(df[v].values, y)[0] for v in VAR_NAMES_ALL]
pvals2 = [stats.pearsonr(df[v].values, y)[1] for v in VAR_NAMES_ALL]

bar_colors2 = []
for v in VAR_NAMES_ALL:
    for cat, vlist in CATEGORIES.items():
        if v in vlist:
            bar_colors2.append(CAT_COLORS[cat])
            break

y_pos2 = np.arange(len(VAR_NAMES_ALL))
var_labels = ['睡眠時間', '趣味娯楽時間\n（ネット代替）', '父親仕事時間', '出生率',
              '通学時間', '教育費', '支援センター\n（万対）', '県民所得', '日照時間']
ax2.barh(y_pos2, corrs, color=bar_colors2, alpha=0.8, edgecolor='white', height=0.6)
for i, (patch, p) in enumerate(zip(ax2.patches, pvals2)):
    if p >= 0.05:
        patch.set_alpha(0.3)
ax2.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax2.set_yticks(y_pos2)
ax2.set_yticklabels(var_labels, fontsize=10)
ax2.set_xlabel('ピアソン相関係数（目的変数：不登校率）', fontsize=11)
ax2.set_title('5カテゴリ別の相関係数（薄色=非有意 p≥0.05）\nデータ: MEXT 2017 + SSDSE-B/D/E/F', fontsize=12, fontweight='bold')
ax2.invert_yaxis()
ax2.grid(axis='x', alpha=0.3)
legend_els = [Patch(color=c, alpha=0.8, label=cat) for cat, c in CAT_COLORS.items()]
ax2.legend(handles=legend_els, fontsize=8, loc='lower right')
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2025_H5_4_fig2_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2025_H5_4_fig2_corr.png")

# ──────────────────────────────────────────────────────────
# 図3: 重回帰係数プロット（小中不登校率）
# ──────────────────────────────────────────────────────────
fig3, ax3 = plt.subplots(figsize=(9, 5))
coefs = model.params[1:]
pvals3 = model.pvalues[1:]
ses3   = model.bse[1:]
bar_cols3 = ['#C62828' if p < 0.05 else '#1565C0' if p < 0.1 else '#9E9E9E' for p in pvals3]
ax3.barh(y_pos2, coefs, color=bar_cols3, alpha=0.75, edgecolor='white', height=0.6)
ax3.errorbar(coefs, y_pos2, xerr=1.96 * ses3, fmt='none', color='#333', capsize=3, linewidth=1.0)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(y_pos2)
ax3.set_yticklabels(var_labels, fontsize=10)
ax3.set_xlabel('回帰係数（±1.96SE）', fontsize=11)
ax3.set_title(f'重回帰分析の係数（不登校率)\nR²={model.rsquared:.3f}, N={N}都道府県', fontsize=12, fontweight='bold')
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)
legend_els3 = [Patch(color='#C62828', alpha=0.75, label='p<0.05'),
               Patch(color='#1565C0', alpha=0.75, label='p<0.10'),
               Patch(color='#9E9E9E', alpha=0.75, label='非有意')]
ax3.legend(handles=legend_els3, fontsize=9)
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2025_H5_4_fig3_coef.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2025_H5_4_fig3_coef.png")

# ──────────────────────────────────────────────────────────
# 図4: 主要変数の散布図（睡眠時間 / 日照時間）
# ──────────────────────────────────────────────────────────
fig4, axes4 = plt.subplots(1, 2, figsize=(13, 5))
fig4.suptitle('生活環境変数と不登校率の関係（2017年度）', fontsize=13, fontweight='bold')

for ax, (vname, xlabel) in zip(axes4,
        [('睡眠時間', '睡眠時間（分/日, SSDSE-D 2021年）'),
         ('日照時間', '年間日照時間（時間, SSDSE-F）')]):
    xv = df[vname].values
    r, p = stats.pearsonr(xv, y)
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'

    # カラー: 外れ値（高不登校率県）を強調
    is_out = (y > upper) | (y < lower)
    ax.scatter(xv[~is_out], y[~is_out], color='#1565C0', s=55, alpha=0.65,
               edgecolors='white', linewidth=0.5, label='通常', zorder=3)
    ax.scatter(xv[is_out], y[is_out], color='#C62828', s=80, alpha=0.85,
               edgecolors='white', linewidth=0.5, label='外れ値', zorder=4)
    for i, pref in enumerate(df['都道府県']):
        if is_out[i]:
            ax.annotate(pref, (xv[i], y[i]), fontsize=7,
                        xytext=(3, 3), textcoords='offset points')

    z = np.polyfit(xv, y, 1)
    xs = np.linspace(xv.min() - 1, xv.max() + 1, 100)
    ax.plot(xs, np.poly1d(z)(xs), 'k--', linewidth=1.5, alpha=0.6)
    ax.set_xlabel(xlabel, fontsize=10)
    ax.set_ylabel('不登校率（千対）', fontsize=10)
    ax.set_title(f'{vname}と不登校率  r={r:.3f} {sig}', fontsize=11, fontweight='bold')
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2025_H5_4_fig4_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("図4保存: 2025_H5_4_fig4_scatter.png")

print("\n全図の生成完了（4枚）")
print("\nデータ出典:")
print("  従属変数: 文部科学省 問題行動等調査 平成29年度（e-Stat: 00400304）")
print("  教育支援センター: 文部科学省 上記調査 表4-15")
print("  SSDSE-B-2026（2017年度）: 出生率, 教育費（統計センター）")
print("  SSDSE-D-2023（2021年度）: 睡眠時間, 通学時間, 仕事時間（統計センター）")
print("  SSDSE-E-2026: 1人当たり県民所得（統計センター）")
print("  SSDSE-F-2023v3: 年間日照時間（統計センター）")