"""
教育用再現コード: 2021年度 統計データ分析コンペティション 審査員奨励賞（大学生の部）
=================================================================================
論文タイトル：女性就業率と賃金格差の地域要因分析
受賞区分  ：審査員奨励賞 ［大学生の部］

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県別パネルデータ）
  対象  ：全47都道府県、2022年断面データ

  本研究は「女性就業率」の地域差を生む要因を都道府県データで分析する。
  保育所整備・高齢化率・求人倍率・合計特殊出生率・消費支出を説明変数とした
  重回帰分析により決定要因を特定し、Ward法クラスタリングで都道府県を類型化する。

  1. 保育所密度 vs 女性就業率 散布図（地域別色分け、都道府県ラベル、回帰線）
  2. 標準化偏回帰係数プロット（95%信頼区間付き横棒グラフ）
  3. Ward法デンドログラム（47都道府県のクラスタリング）
  4. クラスター別主要指標の箱ひげ図

【変数定義】
  女性就業率   = 就職件数（一般）/ 女性15歳以上人口 × 100
                （新規就業活動率として地域労働市場の流動性を示す指標）
  保育所密度   = 保育所等数 / 総人口 × 1000（人口千人あたり）
  高齢化率     = 65歳以上人口 / 総人口 × 100
  求人倍率     = 月間有効求人数（一般）/ 月間有効求職者数（一般）
  合計特殊出生率 = データから直接取得
  消費支出_万  = 消費支出（二人以上の世帯）/ 10000（万円）

【推定手法】
  statsmodels OLS: 標準化偏回帰係数（z-score標準化後の係数）
  AgglomerativeClustering: Ward法、4クラスター

【使用データ】
  SSDSE-B-2026.csv（統計数理研究所 教育用標準データセット）
  https://www.ism.ac.jp/editsection/ssdse/

【図の出力】
  html/figures/2021_U5_2_fig1.png  ... 保育所密度 vs 女性就業率 散布図
  html/figures/2021_U5_2_fig2.png  ... 標準化偏回帰係数プロット
  html/figures/2021_U5_2_fig3.png  ... Ward法デンドログラム
  html/figures/2021_U5_2_fig4.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/2021_U5_2_shorei.py
# ============================================================


import os
import warnings
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import statsmodels.api as sm
from scipy import stats
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import AgglomerativeClustering
from scipy.cluster.hierarchy import dendrogram, linkage

warnings.filterwarnings('ignore')

# ── パス設定 ──────────────────────────────────────────────────────────────────
FIG_DIR = 'html/figures'
DATA_B  = '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,
})

# ── 地域区分マップ ─────────────────────────────────────────────────────────────
region_map = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北',
    '福島県': '北海道・東北',
    '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東',
    '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国',
    '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄',
}
region_colors = {
    '北海道・東北': '#4e9af1',
    '関東':       '#e05c5c',
    '中部':       '#f0a500',
    '近畿':       '#5cb85c',
    '中国・四国':  '#9b59b6',
    '九州・沖縄':  '#f39c12',
}

# ── データ読み込み ─────────────────────────────────────────────────────────────
print("データ読み込み中...")
df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}', na=False)].copy()
df_b['年度'] = df_b['年度'].astype(int)

# 2022年断面データ
df = df_b[df_b['年度'] == 2022].copy()
print(f"2022年断面データ: {len(df)}都道府県")

# ── 変数計算 ─────────────────────────────────────────────────────────────────
df['女性15歳以上人口'] = df['総人口（女）'] - df['15歳未満人口（女）']
df['女性就業率']   = df['就職件数（一般）'] / df['女性15歳以上人口'] * 100
df['保育所密度']   = df['保育所等数'] / df['総人口'] * 1000
df['高齢化率']     = df['65歳以上人口'] / df['総人口'] * 100
df['求人倍率']     = df['月間有効求人数（一般）'] / df['月間有効求職者数（一般）']
df['出生率']       = df['合計特殊出生率'].astype(float)
df['消費支出_万']  = df['消費支出（二人以上の世帯）'] / 10000
df['地域区分']     = df['都道府県'].map(region_map)

# 分析変数リスト
EXOG_VARS = ['保育所密度', '高齢化率', '求人倍率', '出生率', '消費支出_万']
EXOG_LABELS = {
    '保育所密度': '保育所密度\n（人口千人あたり）',
    '高齢化率':   '高齢化率 (%)',
    '求人倍率':   '有効求人倍率',
    '出生率':     '合計特殊出生率',
    '消費支出_万': '消費支出（万円）',
}

# ── 重回帰分析（標準化係数） ───────────────────────────────────────────────────
print("\n重回帰分析（OLS）...")
X_raw = df[EXOG_VARS].copy().astype(float)
y = df['女性就業率'].astype(float)

# 標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_raw)
X_scaled_df = pd.DataFrame(X_scaled, columns=EXOG_VARS, index=df.index)
X_sm = sm.add_constant(X_scaled_df)

model = sm.OLS(y, X_sm).fit()
print(model.summary())

# 係数・CI抽出（定数項を除く）
coef   = model.params[EXOG_VARS].values
ci_low = model.conf_int().loc[EXOG_VARS, 0].values
ci_up  = model.conf_int().loc[EXOG_VARS, 1].values
pvals  = model.pvalues[EXOG_VARS].values
r2     = model.rsquared
r2_adj = model.rsquared_adj
f_stat = model.fvalue
f_pval = model.f_pvalue

# ── クラスタリング ─────────────────────────────────────────────────────────────
print("\nクラスタリング（Ward法, 4クラスター）...")
cluster_vars = ['女性就業率', '保育所密度', '高齢化率', '求人倍率', '消費支出_万']
X_cl = df[cluster_vars].astype(float)
X_cl_scaled = StandardScaler().fit_transform(X_cl)

# Ward法階層クラスタリング（デンドログラム用）
Z = linkage(X_cl_scaled, method='ward')

# 4クラスターに分割
agg = AgglomerativeClustering(n_clusters=4, linkage='ward')
df['クラスター'] = agg.fit_predict(X_cl_scaled) + 1  # 1〜4に変換

print("クラスター別都道府県数:")
print(df['クラスター'].value_counts().sort_index())
print("\nクラスター別女性就業率平均:")
print(df.groupby('クラスター')['女性就業率'].mean().round(3))

# ── 統計量出力 ────────────────────────────────────────────────────────────────
print("\n=== 記述統計 ===")
desc_vars = ['女性就業率', '保育所密度', '高齢化率', '求人倍率', '出生率', '消費支出_万']
print(df[desc_vars].describe().round(3))

print(f"\n=== 重回帰分析結果 ===")
print(f"R² = {r2:.4f}, 調整済みR² = {r2_adj:.4f}")
print(f"F統計量 = {f_stat:.2f}, p値 = {f_pval:.2e}")
for v, c, p in zip(EXOG_VARS, coef, pvals):
    sig = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else 'n.s.'))
    print(f"  {v}: β={c:.4f}, p={p:.4f} {sig}")

# ───────────────────────────────────────────────────────────────────────────────
# 図1: 保育所密度 vs 女性就業率 散布図（地域色分け・都道府県ラベル・回帰線）
# ───────────────────────────────────────────────────────────────────────────────
print("\n図1: 散布図作成中...")
fig1, ax = plt.subplots(figsize=(11, 8))

for region, color in region_colors.items():
    mask = df['地域区分'] == region
    ax.scatter(
        df.loc[mask, '保育所密度'],
        df.loc[mask, '女性就業率'],
        c=color, label=region, s=60, alpha=0.85, zorder=3,
        edgecolors='white', linewidths=0.5,
    )

# 都道府県ラベル
for _, row in df.iterrows():
    ax.annotate(
        row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', ''),
        xy=(row['保育所密度'], row['女性就業率']),
        fontsize=6.5, ha='center', va='bottom',
        xytext=(0, 4), textcoords='offset points', color='#333',
    )

# 回帰直線
x_vals = df['保育所密度'].astype(float)
y_vals = df['女性就業率'].astype(float)
slope, intercept, r_val, p_val, _ = stats.linregress(x_vals, y_vals)
x_line = np.linspace(x_vals.min(), x_vals.max(), 100)
ax.plot(x_line, slope * x_line + intercept, color='#333', linewidth=1.5,
        linestyle='--', alpha=0.7, label=f'回帰線 (r={r_val:.3f}, p<0.001)')

ax.set_xlabel('保育所密度（人口千人あたり）', fontsize=12)
ax.set_ylabel('女性就業率（%）', fontsize=12)
ax.set_title('保育所密度と女性就業率の関係（2022年 都道府県別）', fontsize=13, fontweight='bold')

legend_patches = [mpatches.Patch(color=c, label=r) for r, c in region_colors.items()]
legend_patches.append(
    plt.Line2D([0], [0], color='#333', linestyle='--', linewidth=1.5,
               label=f'回帰線 (r={r_val:.3f})')
)
ax.legend(handles=legend_patches, loc='upper right', fontsize=9, framealpha=0.9)
ax.grid(True, alpha=0.3, linestyle='-')

fig1.tight_layout()
out1 = os.path.join(FIG_DIR, '2021_U5_2_fig1.png')
fig1.savefig(out1, dpi=150, bbox_inches='tight')
plt.close(fig1)
print(f"  -> {out1}")

# ───────────────────────────────────────────────────────────────────────────────
# 図2: 標準化偏回帰係数プロット（95%CI付き横棒グラフ）
# ───────────────────────────────────────────────────────────────────────────────
print("図2: 係数プロット作成中...")
fig2, ax = plt.subplots(figsize=(9, 5.5))

var_labels = [EXOG_LABELS[v] for v in EXOG_VARS]
colors_bar = ['#e05c5c' if p < 0.05 else '#aaaaaa' for p in pvals]

y_pos = np.arange(len(EXOG_VARS))
ax.barh(y_pos, coef, color=colors_bar, alpha=0.8, height=0.5, zorder=3)

# 95%CI エラーバー
xerr_low  = coef - ci_low
xerr_high = ci_up - coef
ax.errorbar(
    coef, y_pos,
    xerr=[xerr_low, xerr_high],
    fmt='none', ecolor='#333', elinewidth=1.5, capsize=5, zorder=4,
)

ax.axvline(0, color='#333', linewidth=1.0, linestyle='-')
ax.set_yticks(y_pos)
ax.set_yticklabels(var_labels, fontsize=10)
ax.set_xlabel('標準化偏回帰係数（β）', fontsize=11)
ax.set_title(f'重回帰分析：女性就業率の決定要因（2022年, N=47)\nR²={r2:.3f}, 調整済みR²={r2_adj:.3f}',
             fontsize=12, fontweight='bold')

# p値ラベル
for i, (c, p) in enumerate(zip(coef, pvals)):
    sig = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else ''))
    x_text = ci_up[i] + 0.01 if c >= 0 else ci_low[i] - 0.01
    ha_text = 'left' if c >= 0 else 'right'
    label_text = f'p={p:.3f}{sig}' if sig else f'p={p:.2f} n.s.'
    ax.text(x_text, i, label_text, va='center', ha=ha_text, fontsize=8.5, color='#333')

sig_patch = mpatches.Patch(color='#e05c5c', alpha=0.8, label='有意（p<0.05）')
ns_patch  = mpatches.Patch(color='#aaaaaa', alpha=0.8, label='非有意')
ax.legend(handles=[sig_patch, ns_patch], loc='lower right', fontsize=9)
ax.grid(True, axis='x', alpha=0.3)

fig2.tight_layout()
out2 = os.path.join(FIG_DIR, '2021_U5_2_fig2.png')
fig2.savefig(out2, dpi=150, bbox_inches='tight')
plt.close(fig2)
print(f"  -> {out2}")

# ───────────────────────────────────────────────────────────────────────────────
# 図3: Ward法デンドログラム
# ───────────────────────────────────────────────────────────────────────────────
print("図3: デンドログラム作成中...")
fig3, ax = plt.subplots(figsize=(16, 6))

pref_names = df['都道府県'].str.replace('県', '').str.replace('府', '').str.replace('都', '').str.replace('道', '北海').values

dendrogram(
    Z,
    labels=pref_names,
    ax=ax,
    leaf_rotation=90,
    leaf_font_size=8,
    color_threshold=Z[-3, 2],  # 4クラスターに対応する閾値
    above_threshold_color='#999999',
)

ax.set_title('都道府県クラスタリング（Ward法、4クラスター）\n変数: 女性就業率・保育所密度・高齢化率・求人倍率・消費支出',
             fontsize=12, fontweight='bold')
ax.set_ylabel('クラスター間距離（Ward法）', fontsize=11)
ax.set_xlabel('都道府県', fontsize=11)
ax.tick_params(axis='x', labelsize=7.5)

# 4クラスター境界線
threshold = Z[-3, 2]
ax.axhline(threshold, color='red', linewidth=1.2, linestyle='--', alpha=0.7,
           label=f'4クラスター切断（距離={threshold:.2f}）')
ax.legend(fontsize=9)

fig3.tight_layout()
out3 = os.path.join(FIG_DIR, '2021_U5_2_fig3.png')
fig3.savefig(out3, dpi=150, bbox_inches='tight')
plt.close(fig3)
print(f"  -> {out3}")

# ───────────────────────────────────────────────────────────────────────────────
# 図4: クラスター別の主要指標分布（箱ひげ図）
# ───────────────────────────────────────────────────────────────────────────────
print("図4: 箱ひげ図作成中...")
box_vars = ['女性就業率', '保育所密度', '高齢化率', '求人倍率', '消費支出_万']
box_titles = ['女性就業率（%）', '保育所密度\n（人口千人あたり）', '高齢化率（%）',
              '有効求人倍率', '消費支出（万円）']

cluster_palette = {1: '#4e9af1', 2: '#e05c5c', 3: '#f0a500', 4: '#5cb85c'}

fig4, axes = plt.subplots(1, 5, figsize=(16, 5.5))

for ax, var, title in zip(axes, box_vars, box_titles):
    data_by_cl = [
        df.loc[df['クラスター'] == cl, var].astype(float).values
        for cl in range(1, 5)
    ]
    bp = ax.boxplot(
        data_by_cl,
        patch_artist=True,
        medianprops=dict(color='white', linewidth=2),
        whiskerprops=dict(linewidth=1.2),
        capprops=dict(linewidth=1.2),
        flierprops=dict(marker='o', markersize=4, linestyle='none',
                        markeredgecolor='#555', markerfacecolor='none'),
    )
    for patch, cl in zip(bp['boxes'], range(1, 5)):
        patch.set_facecolor(cluster_palette[cl])
        patch.set_alpha(0.8)

    ax.set_xticklabels([f'C{i}' for i in range(1, 5)], fontsize=9)
    ax.set_title(title, fontsize=9.5, fontweight='bold')
    ax.grid(True, axis='y', alpha=0.3)

fig4.suptitle('クラスター別主要指標の分布（Ward法 4クラスター, 2022年）',
              fontsize=12, fontweight='bold', y=1.01)

legend_patches4 = [
    mpatches.Patch(color=cluster_palette[cl], alpha=0.8, label=f'クラスター{cl}')
    for cl in range(1, 5)
]
fig4.legend(handles=legend_patches4, loc='lower center', ncol=4, fontsize=9,
            bbox_to_anchor=(0.5, -0.08))

fig4.tight_layout()
out4 = os.path.join(FIG_DIR, '2021_U5_2_fig4.png')
fig4.savefig(out4, dpi=150, bbox_inches='tight')
plt.close(fig4)
print(f"  -> {out4}")

# ── 最終確認 ─────────────────────────────────────────────────────────────────
print("\n=== 生成ファイル確認 ===")
for fp in [out1, out2, out3, out4]:
    size = os.path.getsize(fp) / 1024
    print(f"  {os.path.basename(fp)}: {size:.1f} KB")

print("\n=== クラスター別都道府県一覧 ===")
for cl in range(1, 5):
    prefs = df.loc[df['クラスター'] == cl, '都道府県'].tolist()
    avg_emp = df.loc[df['クラスター'] == cl, '女性就業率'].mean()
    avg_age = df.loc[df['クラスター'] == cl, '高齢化率'].mean()
    avg_nur = df.loc[df['クラスター'] == cl, '保育所密度'].mean()
    print(f"  クラスター{cl} (n={len(prefs)}): {', '.join(prefs)}")
    print(f"    女性就業率={avg_emp:.3f}%, 高齢化率={avg_age:.1f}%, 保育所密度={avg_nur:.3f}")

print("\n完了!")