"""
教育用再現コード: 2020年度（令和2年度） 統計データ分析コンペティション 審査員奨励賞 [高校生の部]
=======================================================================
論文タイトル：都道府県別CO2排出量の決定要因と地球温暖化対策の分析
手法：相関分析、重回帰分析、特化係数（LQ）、地域比較

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県データ）2022年度断面
  目的変数：1人1日当たりのごみ排出量（CO2排出の代理変数）
  説明変数：
    - 光熱・水道費（エネルギー消費の代理）
    - 交通・通信費（自動車等移動由来排出の代理）
    - ごみのリサイクル率（廃棄物対策）
    - 年平均気温（冷暖房需要の代理）
    - 降水日数（気候特性）

【分析ステップ】
  Step1. 都道府県別ごみ排出量ランキング（棒グラフ）
  Step2. 散布図：光熱・水道費 vs 1人1日ごみ排出量（地域色分け）
  Step3. 相関ヒートマップ（環境・生活指標間）
  Step4. 重回帰分析：標準化偏回帰係数プロット

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）
  統計センターより公表の実データ

【データサイエンス学習ポイント】
  1. 特化係数（LQ）による地域産業構造の比較
  2. 代理変数（Proxy Variable）の限界と解釈
  3. 環境クズネッツ曲線（EKC）仮説
  4. 標準化偏回帰係数による変数重要度の比較
=======================================================================
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・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/2020_H5_5_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
import warnings
warnings.filterwarnings('ignore')

# ──────────────────────────────────────────────────────────────
# 共通設定
# ──────────────────────────────────────────────────────────────
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)

# ================================================================
# ■ 実データ読み込み（SSDSE-B: 都道府県データ）
# ================================================================
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)

# ================================================================
# ■ 分析変数の準備
# ================================================================
# 目的変数: 1人1日当たりごみ排出量（g/人・日）= CO2排出の代理変数
df['y_gomi_pp']    = df['1人1日当たりの排出量']

# 説明変数
df['x_konetsu']    = df['光熱・水道費（二人以上の世帯）']      # エネルギー消費代理
df['x_kotsu']      = df['交通・通信費（二人以上の世帯）']       # 自動車利用代理
df['x_recycle']    = df['ごみのリサイクル率']                   # 廃棄物対策
df['x_temp']       = df['年平均気温']                           # 気候（冷暖房需要）
df['x_rain']       = df['降水日数（年間）']                     # 降水日数

# 分析データセット
ana_cols = ['都道府県', '地域区分', '地域色',
            'y_gomi_pp', 'x_konetsu', 'x_kotsu', 'x_recycle', 'x_temp', 'x_rain']
dfa = df[ana_cols].dropna().copy()
print(f"[分析対象] {len(dfa)}都道府県（欠損除外後）")

# 基本統計量の表示
print("\n=== 基本統計量 ===")
desc_cols = {
    '1人1日ごみ排出量 (g/人・日)': 'y_gomi_pp',
    '光熱・水道費 (円/月)': 'x_konetsu',
    '交通・通信費 (円/月)': 'x_kotsu',
    'ごみリサイクル率 (%)': 'x_recycle',
    '年平均気温 (℃)': 'x_temp',
    '降水日数 (日/年)': 'x_rain',
}
for label, col in desc_cols.items():
    s = dfa[col]
    print(f"  {label}: 平均={s.mean():.1f}, 標準偏差={s.std():.1f}, "
          f"最小={s.min():.1f}, 最大={s.max():.1f}")

# ================================================================
# ■ 相関係数（Pearson）
# ================================================================
print("\n=== 相関分析（目的変数: 1人1日ごみ排出量） ===")
corr_vars = ['x_konetsu', 'x_kotsu', 'x_recycle', 'x_temp', 'x_rain']
corr_labels = {
    'x_konetsu': '光熱・水道費',
    'x_kotsu':   '交通・通信費',
    'x_recycle': 'リサイクル率',
    'x_temp':    '年平均気温',
    'x_rain':    '降水日数',
}
corr_results = {}
for col in corr_vars:
    r, p = stats.pearsonr(dfa['y_gomi_pp'], dfa[col])
    corr_results[col] = (r, p)
    sig = '**' if p < 0.01 else ('*' if p < 0.05 else 'n.s.')
    print(f"  {corr_labels[col]:12s}: r={r:.3f}, p={p:.4f} {sig}")

# ================================================================
# ■ 重回帰分析 (OLS)
# ================================================================
X_cols = corr_vars
X = dfa[X_cols].values
y = dfa['y_gomi_pp'].values

# 標準化
X_mean = X.mean(axis=0)
X_std  = X.std(axis=0, ddof=1)
y_mean = y.mean()
y_std  = y.std(ddof=1)

X_std_arr = (X - X_mean) / X_std
y_std_arr = (y - y_mean) / y_std

X_with_const = sm.add_constant(X_std_arr)
model = sm.OLS(y_std_arr, X_with_const).fit()

print("\n=== 重回帰分析結果 (標準化係数) ===")
print(f"  R² = {model.rsquared:.4f}, 調整済みR² = {model.rsquared_adj:.4f}")
print(f"  F値 = {model.fvalue:.3f}, p値 = {model.f_pvalue:.6f}")
print()
coef_labels = ['定数項'] + [corr_labels[c] for c in X_cols]
for i, (label, coef, pval) in enumerate(zip(coef_labels,
                                             model.params,
                                             model.pvalues)):
    sig = '**' if pval < 0.01 else ('*' if pval < 0.05 else 'n.s.')
    print(f"  {label:14s}: β={coef:+.4f}, p={pval:.4f} {sig}")

# ================================================================
# ■ 図1: 散布図（光熱・水道費 vs 1人1日ごみ排出量、地域色分け）
# ================================================================
fig1, ax1 = plt.subplots(figsize=(9, 6))

for region, color in region_colors.items():
    sub = dfa[dfa['地域区分'] == region]
    ax1.scatter(sub['x_konetsu'], sub['y_gomi_pp'],
                color=color, label=region, s=60, zorder=3, alpha=0.9)
    for _, row in sub.iterrows():
        pref = row['都道府県'].replace('県', '').replace('都', '').replace('府', '').replace('道', '')
        ax1.annotate(pref,
                     xy=(row['x_konetsu'], row['y_gomi_pp']),
                     xytext=(3, 3), textcoords='offset points',
                     fontsize=6.5, color='#444', zorder=4)

# 回帰直線
x_line = np.linspace(dfa['x_konetsu'].min(), dfa['x_konetsu'].max(), 100)
r_val, p_val = stats.pearsonr(dfa['x_konetsu'], dfa['y_gomi_pp'])
slope, intercept, *_ = stats.linregress(dfa['x_konetsu'], dfa['y_gomi_pp'])
y_line = intercept + slope * x_line
ax1.plot(x_line, y_line, color='#333', linewidth=1.5, linestyle='--',
         label=f'回帰直線 (r={r_val:.2f})')

ax1.set_xlabel('光熱・水道費（円/月）\n（エネルギー消費の代理変数）', fontsize=11)
ax1.set_ylabel('1人1日当たりごみ排出量（g/人・日）\n（CO2排出の代理変数）', fontsize=11)
ax1.set_title('光熱・水道費と1人1日ごみ排出量の散布図\n（2022年度 都道府県別）', fontsize=12, fontweight='bold')
ax1.legend(loc='upper left', fontsize=9, framealpha=0.8)
ax1.grid(True, alpha=0.3)

fig1.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2020_H5_5_fig1.png'), dpi=150, bbox_inches='tight')
plt.close(fig1)
print("\n[図1] 散布図 → 保存完了")

# ================================================================
# ■ 図2: 都道府県別ごみ排出量ランキング（棒グラフ）
# ================================================================
dfa_sorted = dfa.sort_values('y_gomi_pp', ascending=False).reset_index(drop=True)

fig2, ax2 = plt.subplots(figsize=(12, 7))
bars = ax2.bar(range(len(dfa_sorted)),
               dfa_sorted['y_gomi_pp'],
               color=dfa_sorted['地域色'].values,
               edgecolor='white', linewidth=0.5, zorder=3)

# 都道府県ラベル（X軸）
pref_labels = [p.replace('県', '').replace('都', '').replace('府', '').replace('道', '')
               for p in dfa_sorted['都道府県']]
ax2.set_xticks(range(len(dfa_sorted)))
ax2.set_xticklabels(pref_labels, rotation=90, fontsize=8)

# 凡例（地域区分）
from matplotlib.patches import Patch
legend_patches = [Patch(color=c, label=r) for r, c in region_colors.items()]
ax2.legend(handles=legend_patches, loc='upper right', fontsize=9, framealpha=0.8)

# 全国平均線
mean_val = dfa['y_gomi_pp'].mean()
ax2.axhline(mean_val, color='#333', linestyle='--', linewidth=1.5, zorder=4)
ax2.text(len(dfa_sorted) - 0.5, mean_val + 5,
         f'全国平均\n{mean_val:.0f}g', ha='right', fontsize=8.5, color='#333')

ax2.set_xlabel('都道府県', fontsize=11)
ax2.set_ylabel('1人1日当たりごみ排出量（g/人・日）', fontsize=11)
ax2.set_title('都道府県別 1人1日当たりごみ排出量ランキング\n（2022年度・CO2排出の代理変数）', fontsize=12, fontweight='bold')
ax2.set_xlim(-0.7, len(dfa_sorted) - 0.3)
ax2.set_ylim(700, max(dfa_sorted['y_gomi_pp']) * 1.12)
ax2.grid(True, axis='y', alpha=0.3, zorder=0)

fig2.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2020_H5_5_fig2.png'), dpi=150, bbox_inches='tight')
plt.close(fig2)
print("[図2] ランキング棒グラフ → 保存完了")

# ================================================================
# ■ 図3: 相関ヒートマップ
# ================================================================
heat_vars = ['y_gomi_pp', 'x_konetsu', 'x_kotsu', 'x_recycle', 'x_temp', 'x_rain']
heat_labels = ['1人1日\nごみ排出量', '光熱・\n水道費', '交通・\n通信費', 'リサイ\nクル率', '年平均\n気温', '降水\n日数']

corr_matrix = dfa[heat_vars].corr()

fig3, ax3 = plt.subplots(figsize=(7, 6))
n = len(heat_vars)
im = ax3.imshow(corr_matrix.values, cmap='RdBu_r', vmin=-1, vmax=1)
plt.colorbar(im, ax=ax3, shrink=0.85, label='Pearson r')

ax3.set_xticks(range(n))
ax3.set_yticks(range(n))
ax3.set_xticklabels(heat_labels, fontsize=9)
ax3.set_yticklabels(heat_labels, fontsize=9)

for i in range(n):
    for j in range(n):
        val = corr_matrix.values[i, j]
        r_ij, p_ij = stats.pearsonr(dfa[heat_vars[i]], dfa[heat_vars[j]])
        sig_mark = '**' if p_ij < 0.01 else ('*' if p_ij < 0.05 else '')
        txt_color = 'white' if abs(val) > 0.65 else 'black'
        ax3.text(j, i, f'{val:.2f}{sig_mark}', ha='center', va='center',
                 fontsize=8, color=txt_color, fontweight='bold' if sig_mark else 'normal')

ax3.set_title('環境・生活指標の相関ヒートマップ\n（Pearson r、*p<0.05, **p<0.01）', fontsize=11, fontweight='bold')
fig3.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2020_H5_5_fig3.png'), dpi=150, bbox_inches='tight')
plt.close(fig3)
print("[図3] 相関ヒートマップ → 保存完了")

# ================================================================
# ■ 図4: 標準化偏回帰係数プロット
# ================================================================
coef_names = [corr_labels[c] for c in X_cols]
coef_vals   = model.params[1:]   # 定数項除く
pvals_vals  = model.pvalues[1:]
conf_int    = np.array(model.conf_int())[1:]  # 信頼区間

# 係数の絶対値で降順ソート
order = np.argsort(np.abs(coef_vals))[::-1]
coef_vals_s  = coef_vals[order]
pvals_vals_s = pvals_vals[order]
names_s      = [coef_names[i] for i in order]
ci_s         = conf_int[order]

colors_bar = ['#e05c5c' if p < 0.05 else '#aaa' for p in pvals_vals_s]

fig4, ax4 = plt.subplots(figsize=(7, 5))
y_pos = range(len(names_s))
ax4.barh(y_pos, coef_vals_s, color=colors_bar, edgecolor='white', height=0.6, zorder=3)

# 信頼区間（95%CI）
for i, (ci_l, ci_u) in enumerate(ci_s):
    ax4.plot([ci_l, ci_u], [i, i], color='#333', linewidth=1.5, zorder=4)

ax4.set_yticks(list(y_pos))
ax4.set_yticklabels(names_s, fontsize=11)
ax4.axvline(0, color='black', linewidth=0.8)
ax4.set_xlabel('標準化偏回帰係数 β', fontsize=11)
ax4.set_title('重回帰分析：標準化偏回帰係数\n（赤：p<0.05有意、灰色：非有意）', fontsize=11, fontweight='bold')
ax4.grid(True, axis='x', alpha=0.3, zorder=0)

# p値ラベル
for i, (coef, pval) in enumerate(zip(coef_vals_s, pvals_vals_s)):
    sig = '**' if pval < 0.01 else ('*' if pval < 0.05 else 'n.s.')
    xpos = coef + (0.02 if coef >= 0 else -0.02)
    ha = 'left' if coef >= 0 else 'right'
    ax4.text(xpos, i, sig, va='center', ha=ha, fontsize=9, color='#333')

# 凡例
from matplotlib.patches import Patch
ax4.legend(handles=[
    Patch(color='#e05c5c', label='有意 (p<0.05)'),
    Patch(color='#aaa',    label='非有意')
], loc='lower right', fontsize=9)

fig4.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2020_H5_5_fig4.png'), dpi=150, bbox_inches='tight')
plt.close(fig4)
print("[図4] 標準化偏回帰係数プロット → 保存完了")

# ================================================================
# ■ 追加統計：地域別平均（地域比較用）
# ================================================================
print("\n=== 地域別 1人1日ごみ排出量 平均 ===")
region_stats = dfa.groupby('地域区分')['y_gomi_pp'].agg(['mean', 'std', 'count'])
region_stats.columns = ['平均', '標準偏差', '都道府県数']
print(region_stats.sort_values('平均', ascending=False).to_string())

# 上位・下位都道府県
print("\n=== 上位5都道府県（ごみ排出量多） ===")
top5 = dfa_sorted.head(5)[['都道府県', 'y_gomi_pp', '地域区分']]
top5.columns = ['都道府県', '1人1日排出量(g)', '地域区分']
print(top5.to_string(index=False))

print("\n=== 下位5都道府県（ごみ排出量少） ===")
bot5 = dfa_sorted.tail(5)[['都道府県', 'y_gomi_pp', '地域区分']]
bot5.columns = ['都道府県', '1人1日排出量(g)', '地域区分']
print(bot5.to_string(index=False))

# 特化係数（LQ: Location Quotient）の計算例
# LQ = (都道府県iの光熱水道費シェア) / (全国の光熱水道費シェア)
# ここでは光熱水道費の全国消費支出に占める割合の特化係数を計算
df_lq = df.copy()
national_share = df_lq['光熱・水道費（二人以上の世帯）'].mean() / df_lq['消費支出（二人以上の世帯）'].mean()
df_lq['LQ_konetsu'] = (df_lq['光熱・水道費（二人以上の世帯）'] / df_lq['消費支出（二人以上の世帯）']) / national_share
print("\n=== 光熱・水道費特化係数（LQ）上位10都道府県 ===")
lq_top = df_lq[['都道府県', 'LQ_konetsu']].sort_values('LQ_konetsu', ascending=False).head(10)
print(lq_top.to_string(index=False))

print("\n[完了] 全4枚の図を保存しました。")
print(f"  保存先: {FIG_DIR}")