"""
2021年度 統計データ分析コンペティション 審査員奨励賞（高校生の部）
「都道府県別環境指標と住民生活の質の関係分析」
教育用再現コード — 実データ（SSDSE-B-2026）使用
"""

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


import os
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

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

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

# ── データ読み込み ──────────────────────────────────────────────
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().reset_index(drop=True)
print(f"2022年データ: {len(df)}都道府県")

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

df['地域'] = df['都道府県'].map(region_map)
df['地域色'] = df['地域'].map(region_colors)

# ── 変数の数値変換 ──────────────────────────────────────────────
COLS = {
    'ごみ排出量':    '1人1日当たりの排出量',
    'リサイクル率':  'ごみのリサイクル率',
    '消費支出':      '消費支出（二人以上の世帯）',
    '教養娯楽費':    '教養娯楽費（二人以上の世帯）',
    '出生率':        '合計特殊出生率',
    '気温':          '年平均気温',
    '降水量':        '降水量（年間）',
    '病院数':        '一般病院数',
}

for key, col in COLS.items():
    df[key] = pd.to_numeric(df[col], errors='coerce')

# 人口1万人あたり病院数
df['人口'] = pd.to_numeric(df['総人口'], errors='coerce')
df['病院数_per10k'] = df['病院数'] / df['人口'] * 10000

# 分析用データ（欠損除去）
ENV_VARS  = ['ごみ排出量', 'リサイクル率', '気温', '降水量']
LIFE_VARS = ['消費支出', '教養娯楽費', '出生率']
ALL_VARS  = ENV_VARS + LIFE_VARS
df_ana = df[['都道府県', '地域', '地域色'] + ALL_VARS].dropna().copy()
print(f"分析対象: {len(df_ana)}都道府県")

# ── 統計サマリー出力 ────────────────────────────────────────────
print("\n=== 基本統計量（2022年） ===")
for v in ALL_VARS:
    s = df_ana[v]
    print(f"  {v}: 平均={s.mean():.2f}, 標準偏差={s.std():.2f}, "
          f"最小={s.min():.2f}, 最大={s.max():.2f}")

# ── 相関分析 ────────────────────────────────────────────────────
print("\n=== 消費支出との相関（Pearson r）===")
for v in ENV_VARS + ['教養娯楽費', '出生率']:
    r, p = stats.pearsonr(df_ana[v], df_ana['消費支出'])
    sig = '**' if p < 0.01 else ('*' if p < 0.05 else 'ns')
    print(f"  {v}: r={r:.4f}, p={p:.4f} {sig}")

# ── OLS重回帰 ───────────────────────────────────────────────────
X_cols = ['ごみ排出量', 'リサイクル率', '気温', '降水量', '教養娯楽費', '出生率']
y_col  = '消費支出'

X_raw = df_ana[X_cols].copy()
y_raw = df_ana[y_col].copy()

# 標準化
X_std = (X_raw - X_raw.mean()) / X_raw.std()
y_std = (y_raw - y_raw.mean()) / y_raw.std()

X_ols = sm.add_constant(X_std)
model = sm.OLS(y_std, X_ols).fit()
print("\n=== OLS重回帰結果（標準化係数）===")
print(f"  R² = {model.rsquared:.4f}")
print(f"  Adj. R² = {model.rsquared_adj:.4f}")
print(f"  F-stat = {model.fvalue:.4f}, p = {model.f_pvalue:.4f}")
for name, coef, pval in zip(X_cols, model.params[1:], model.pvalues[1:]):
    sig = '**' if pval < 0.01 else ('*' if pval < 0.05 else 'ns')
    print(f"  {name}: β={coef:.4f}, p={pval:.4f} {sig}")

r2_val      = model.rsquared
adj_r2_val  = model.rsquared_adj
f_val       = model.fvalue
f_pval      = model.f_pvalue
coef_dict   = dict(zip(X_cols, model.params[1:]))
pval_dict   = dict(zip(X_cols, model.pvalues[1:]))

# ── 地域別ごみ排出量統計 ────────────────────────────────────────
print("\n=== 地域別ごみ排出量（1人1日あたり g）===")
for region in region_colors:
    sub = df_ana[df_ana['地域'] == region]['ごみ排出量']
    if len(sub) > 0:
        print(f"  {region}: 平均={sub.mean():.1f}g, n={len(sub)}")

# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
# 図1: ごみ排出量 vs 消費支出（地域色分け・都道府県ラベル・回帰線）
# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
fig, ax = plt.subplots(figsize=(10, 7))

for region, color in region_colors.items():
    mask = df_ana['地域'] == region
    sub  = df_ana[mask]
    ax.scatter(sub['ごみ排出量'], sub['消費支出'] / 10000,
               color=color, label=region, s=60, zorder=3, alpha=0.85, edgecolors='white', linewidth=0.5)

# 都道府県ラベル
for _, row in df_ana.iterrows():
    ax.annotate(row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', ''),
                (row['ごみ排出量'], row['消費支出'] / 10000),
                fontsize=7, ha='left', va='bottom',
                xytext=(3, 2), textcoords='offset points', color='#333333')

# 回帰線
x_plot = np.linspace(df_ana['ごみ排出量'].min(), df_ana['ごみ排出量'].max(), 100)
slope, intercept, r_val, p_val_reg, _ = stats.linregress(df_ana['ごみ排出量'], df_ana['消費支出'] / 10000)
ax.plot(x_plot, intercept + slope * x_plot, color='#333333', linewidth=1.5,
        linestyle='--', label=f'回帰線 r={r_val:.3f}', zorder=2)

ax.set_xlabel('1人1日あたりごみ排出量（g）', fontsize=12)
ax.set_ylabel('消費支出（万円/月、二人以上の世帯）', fontsize=12)
ax.set_title('都道府県別：ごみ排出量と消費支出の関係（2022年）', fontsize=13, fontweight='bold')
ax.legend(fontsize=9, loc='upper right', framealpha=0.9)
ax.grid(True, linestyle=':', alpha=0.5)
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2021_H5_6_fig1.png'), bbox_inches='tight')
plt.close()
print("\n図1保存完了")

# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
# 図2: 相関行列ヒートマップ
# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
HM_VARS = ['ごみ排出量', 'リサイクル率', '気温', '降水量', '消費支出', '教養娯楽費', '出生率']
HM_LABELS = ['ごみ排出量', 'リサイクル率', '年平均気温', '年間降水量', '消費支出', '教養娯楽費', '合計特殊\n出生率']

corr_mat = df_ana[HM_VARS].corr()
pval_mat = pd.DataFrame(np.ones((len(HM_VARS), len(HM_VARS))), index=HM_VARS, columns=HM_VARS)
for i, c1 in enumerate(HM_VARS):
    for j, c2 in enumerate(HM_VARS):
        if i != j:
            _, p = stats.pearsonr(df_ana[c1], df_ana[c2])
            pval_mat.loc[c1, c2] = p

fig, ax = plt.subplots(figsize=(9, 7))
im = ax.imshow(corr_mat.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax, shrink=0.8, label='Pearson r')

n = len(HM_VARS)
ax.set_xticks(range(n))
ax.set_yticks(range(n))
ax.set_xticklabels(HM_LABELS, fontsize=9, rotation=35, ha='right')
ax.set_yticklabels(HM_LABELS, fontsize=9)

for i in range(n):
    for j in range(n):
        r_ij = corr_mat.values[i, j]
        p_ij = pval_mat.values[i, j]
        sig_str = '**' if p_ij < 0.01 else ('*' if p_ij < 0.05 else '')
        txt = f'{r_ij:.2f}{sig_str}'
        ax.text(j, i, txt, ha='center', va='center', fontsize=8,
                color='white' if abs(r_ij) > 0.5 else 'black')

ax.set_title('環境・生活変数の相関行列（2022年、47都道府県）\n* p<0.05, ** p<0.01', fontsize=12, fontweight='bold')
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2021_H5_6_fig2.png'), bbox_inches='tight')
plt.close()
print("図2保存完了")

# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
# 図3: 標準化偏回帰係数プロット
# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
fig, ax = plt.subplots(figsize=(8, 5))

labels_jp = {
    'ごみ排出量': 'ごみ排出量\n（環境負荷↑）',
    'リサイクル率': 'リサイクル率\n（環境意識↑）',
    '気温': '年平均気温',
    '降水量': '年間降水量',
    '教養娯楽費': '教養娯楽費\n（余暇支出）',
    '出生率': '合計特殊出生率',
}

coefs  = [coef_dict[v] for v in X_cols]
pvals  = [pval_dict[v] for v in X_cols]
labels = [labels_jp[v] for v in X_cols]
colors = ['#C62828' if c > 0 else '#1565C0' for c in coefs]
alphas = [1.0 if p < 0.05 else 0.45 for p in pvals]

bars = ax.barh(range(len(X_cols)), coefs, color=colors, alpha=0.75, height=0.6)
for i, (bar, alpha) in enumerate(zip(bars, alphas)):
    bar.set_alpha(alpha)

# 信頼区間（95%CI）
ci = model.conf_int(alpha=0.05)
ci_lo = ci.iloc[1:, 0].values
ci_hi = ci.iloc[1:, 1].values
for i in range(len(X_cols)):
    ax.plot([ci_lo[i], ci_hi[i]], [i, i], color='black', linewidth=1.5, zorder=5)
    ax.plot([ci_lo[i]], [i], '|', color='black', markersize=6)
    ax.plot([ci_hi[i]], [i], '|', color='black', markersize=6)

# p値注釈
for i, (c, p) in enumerate(zip(coefs, pvals)):
    sig_str = '**' if p < 0.01 else ('*' if p < 0.05 else 'ns')
    offset = 0.02 if c >= 0 else -0.02
    ha = 'left' if c >= 0 else 'right'
    ax.text(c + offset, i, f' {sig_str}', ha=ha, va='center', fontsize=9,
            color='#C62828' if sig_str != 'ns' else '#888888')

ax.set_yticks(range(len(X_cols)))
ax.set_yticklabels(labels, fontsize=9)
ax.axvline(0, color='black', linewidth=0.8)
ax.set_xlabel('標準化偏回帰係数（β）', fontsize=11)
ax.set_title(f'消費支出を目的変数とするOLS重回帰\n標準化偏回帰係数（R²={r2_val:.3f}、Adj.R²={adj_r2_val:.3f}）\n* p<0.05, ** p<0.01, ns=非有意', fontsize=11, fontweight='bold')
ax.grid(True, axis='x', linestyle=':', alpha=0.5)
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2021_H5_6_fig3.png'), bbox_inches='tight')
plt.close()
print("図3保存完了")

# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
# 図4: 都道府県別ごみ排出量ランキング（棒グラフ・地域色分け）
# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
df_rank = df_ana[['都道府県', '地域', 'ごみ排出量']].sort_values('ごみ排出量', ascending=True).reset_index(drop=True)

fig, ax = plt.subplots(figsize=(10, 12))
bar_colors = [region_colors[r] for r in df_rank['地域']]
ax.barh(range(len(df_rank)), df_rank['ごみ排出量'], color=bar_colors, alpha=0.85, height=0.7, edgecolor='white', linewidth=0.3)

ax.set_yticks(range(len(df_rank)))
pref_labels = df_rank['都道府県'].str.replace('県', '').str.replace('府', '').str.replace('都', '').str.replace('道', '')
ax.set_yticklabels(pref_labels, fontsize=8.5)

# 全国平均線
avg_gomi = df_ana['ごみ排出量'].mean()
ax.axvline(avg_gomi, color='#333333', linewidth=1.5, linestyle='--', label=f'全国平均 {avg_gomi:.0f}g')

# 数値ラベル
for i, val in enumerate(df_rank['ごみ排出量']):
    ax.text(val + 5, i, f'{val:.0f}', va='center', fontsize=7.5, color='#333')

ax.set_xlabel('1人1日あたりごみ排出量（g）', fontsize=11)
ax.set_title('都道府県別 1人1日あたりごみ排出量（2022年）\n（地域色分け）', fontsize=12, fontweight='bold')
ax.legend(fontsize=9, loc='lower right')

# 地域凡例
from matplotlib.patches import Patch
legend_patches = [Patch(facecolor=c, label=r) for r, c in region_colors.items()]
ax.legend(handles=legend_patches, fontsize=8.5, loc='lower right', framealpha=0.9, title='地域', title_fontsize=9)

ax.grid(True, axis='x', linestyle=':', alpha=0.4)
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2021_H5_6_fig4.png'), bbox_inches='tight')
plt.close()
print("図4保存完了")

# ── 最終サマリー出力（HTML用） ──────────────────────────────────
print("\n=== HTML記載用 統計値サマリー ===")
gomi_r, gomi_p = stats.pearsonr(df_ana['ごみ排出量'], df_ana['消費支出'])
recycl_r, recycl_p = stats.pearsonr(df_ana['リサイクル率'], df_ana['消費支出'])
edu_r, edu_p = stats.pearsonr(df_ana['教養娯楽費'], df_ana['消費支出'])
birth_r, birth_p = stats.pearsonr(df_ana['出生率'], df_ana['消費支出'])
print(f"  ごみ排出量 vs 消費支出: r={gomi_r:.3f}, p={gomi_p:.4f}")
print(f"  リサイクル率 vs 消費支出: r={recycl_r:.3f}, p={recycl_p:.4f}")
print(f"  教養娯楽費 vs 消費支出: r={edu_r:.3f}, p={edu_p:.4f}")
print(f"  出生率 vs 消費支出: r={birth_r:.3f}, p={birth_p:.4f}")
print(f"  OLS R²={r2_val:.3f}, Adj.R²={adj_r2_val:.3f}")
print(f"  F({model.df_model:.0f},{model.df_resid:.0f})={f_val:.2f}, p={f_pval:.4f}")
print(f"  全国平均 ごみ排出量: {df_ana['ごみ排出量'].mean():.1f}g")
print(f"  全国平均 消費支出: {df_ana['消費支出'].mean():.0f}円/月")
print(f"  全国平均 合計特殊出生率: {df_ana['出生率'].mean():.3f}")
# 最大・最小ごみ排出量都道府県
idx_max = df_ana['ごみ排出量'].idxmax()
idx_min = df_ana['ごみ排出量'].idxmin()
print(f"  ごみ排出量 最大: {df_ana.loc[idx_max,'都道府県']} {df_ana.loc[idx_max,'ごみ排出量']:.0f}g")
print(f"  ごみ排出量 最小: {df_ana.loc[idx_min,'都道府県']} {df_ana.loc[idx_min,'ごみ排出量']:.0f}g")

print("\n全ての図の生成が完了しました。")