"""
2021年度 統計データ分析コンペティション 審査員奨励賞（大学生の部）
「インターネット利用率・ICT環境と地域経済指標の関係分析」
教育用再現スクリプト（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_U5_4_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
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

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)
assert len(df) == 47, f"47都道府県を期待したが {len(df)} 行"

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

# ── ICT関連プロキシ変数の構築 ─────────────────────────────────────────────────
# 大学進学率 = 高校卒業者のうち進学者 / 高校卒業者数
df['大学進学率'] = df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数']

# 大学数per万人（高等教育集積度）
df['大学数per万人'] = df['大学数'] / df['総人口'] * 10000

# 大学生比率（高等教育在学率）
df['大学生比率'] = df['大学学生数'] / df['総人口']

# 通信費支出割合（ICT消費の代理）
df['通信費比率'] = df['交通・通信費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）']

# 教育文化費比率（教育・文化投資の代理）
df['教育文化費比率'] = (
    df['教育費（二人以上の世帯）'] + df['教養娯楽費（二人以上の世帯）']
) / df['消費支出（二人以上の世帯）']

# 経済指標（目的変数候補）
df['消費水準'] = df['消費支出（二人以上の世帯）']  # 万円換算
df['通信費絶対額'] = df['交通・通信費（二人以上の世帯）']

# ICT分析に使う変数群
ict_vars = ['大学進学率', '大学数per万人', '大学生比率', '通信費比率', '教育文化費比率']
ict_labels = ['大学進学率', '大学数\n(per万人)', '大学生比率', '通信費比率', '教育文化費比率']

eco_vars = ['消費支出（二人以上の世帯）', '教育費（二人以上の世帯）',
            '教養娯楽費（二人以上の世帯）']

# ── PCA（ICT活用指数の構築） ─────────────────────────────────────────────────
X_ict = df[ict_vars].values
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_ict)

pca = PCA(n_components=2)
pca_scores = pca.fit_transform(X_scaled)
df['PC1'] = pca_scores[:, 0]
df['PC2'] = pca_scores[:, 1]

# ICT活用指数 = PC1（第一主成分、ICT集積度を主に反映）
# PC1が大きいほどICT先進
df['ICT指数'] = df['PC1']
# PC1の符号を確認（大学進学率などと同方向になるよう調整）
if pca.components_[0][0] < 0:  # 大学進学率の係数が負なら反転
    df['ICT指数'] = -df['ICT指数']
    pca_scores[:, 0] = -pca_scores[:, 0]

# 都道府県分類
q33 = df['ICT指数'].quantile(0.67)
q33_lo = df['ICT指数'].quantile(0.33)
df['ICT分類'] = pd.cut(
    df['ICT指数'],
    bins=[-np.inf, q33_lo, q33, np.inf],
    labels=['ICT後発県', '標準県', 'ICT先進県']
)

# ── 重回帰分析 ───────────────────────────────────────────────────────────────
# 目的変数: 消費水準（豊かさの代理）
# 説明変数: ICT関連5変数（標準化）
X_reg = X_scaled.copy()
y_reg = df['消費支出（二人以上の世帯）'].values

X_reg_const = sm.add_constant(X_reg)
ols_model = sm.OLS(y_reg, X_reg_const).fit()

print("=" * 60)
print("重回帰分析結果 (目的変数: 消費支出)")
print("=" * 60)
print(ols_model.summary())

coef_vals = ols_model.params[1:]       # 定数項除く
_ci_raw   = ols_model.conf_int()
conf_int  = (_ci_raw.values if hasattr(_ci_raw, 'values') else _ci_raw)[1:]   # ndarray (n_vars, 2)
pvals     = ols_model.pvalues[1:]
r2        = ols_model.rsquared
adj_r2    = ols_model.rsquared_adj

print(f"\nR² = {r2:.4f}, Adj-R² = {adj_r2:.4f}")

# ── 相関分析 ─────────────────────────────────────────────────────────────────
corr_cols = ict_vars + ['消費支出（二人以上の世帯）', '教育費（二人以上の世帯）', '教養娯楽費（二人以上の世帯）']
corr_labels = ['大学進学率', '大学数\nper万人', '大学生\n比率', '通信費\n比率', '教育文化\n費比率',
               '消費支出', '教育費', '教養娯\n楽費']
corr_matrix = df[corr_cols].corr(method='pearson')

# ── 上位・下位ランキング ──────────────────────────────────────────────────────
df_rank = df[['都道府県', 'ICT指数', 'ICT分類', '地域']].copy()
df_rank = df_rank.sort_values('ICT指数', ascending=False).reset_index(drop=True)
top10 = df_rank.head(10)
bot10 = df_rank.tail(10)

print("\nICT先進県 TOP10:")
print(top10[['都道府県', 'ICT指数']].to_string(index=False))
print("\nICT後発県 BOTTOM10:")
print(bot10[['都道府県', 'ICT指数']].to_string(index=False))

# PCAの寄与率
var_expl = pca.explained_variance_ratio_
print(f"\nPC1寄与率: {var_expl[0]:.3f}, PC2寄与率: {var_expl[1]:.3f}")
print(f"累積寄与率: {sum(var_expl):.3f}")

# ══════════════════════════════════════════════════════════════════════════════
# 図1: 相関ヒートマップ
# ══════════════════════════════════════════════════════════════════════════════
fig1, ax1 = plt.subplots(figsize=(9, 7))

cmap = plt.get_cmap('RdBu_r')
im = ax1.imshow(corr_matrix.values, cmap=cmap, vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax1, label='Pearson相関係数')

n = len(corr_labels)
ax1.set_xticks(range(n))
ax1.set_yticks(range(n))
ax1.set_xticklabels(corr_labels, fontsize=10)
ax1.set_yticklabels(corr_labels, fontsize=10)
ax1.tick_params(axis='x', rotation=0)

# 相関係数の数値を表示
for i in range(n):
    for j in range(n):
        val = corr_matrix.values[i, j]
        color = 'white' if abs(val) > 0.6 else 'black'
        ax1.text(j, i, f'{val:.2f}', ha='center', va='center',
                 fontsize=9, color=color, fontweight='bold')

ax1.set_title('図1: ICT関連変数と経済指標の相関行列（2022年・47都道府県）\nPearson相関係数',
              fontsize=13, pad=14)
ax1.set_xlim(-0.5, n - 0.5)
ax1.set_ylim(n - 0.5, -0.5)
fig1.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2021_U5_4_fig1.png'), bbox_inches='tight')
plt.close(fig1)
print("図1 保存完了")

# ══════════════════════════════════════════════════════════════════════════════
# 図2: PCAバイプロット (PC1 vs PC2、地域色分け・ラベル・変数ベクトル)
# ══════════════════════════════════════════════════════════════════════════════
fig2, ax2 = plt.subplots(figsize=(11, 9))

for region, color in region_colors.items():
    mask = df['地域'] == region
    ax2.scatter(
        pca_scores[mask, 0], pca_scores[mask, 1],
        c=color, label=region, s=70, zorder=3, edgecolors='white', linewidths=0.5
    )

# 都道府県ラベル
for i, row in df.iterrows():
    ax2.annotate(
        row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', ''),
        (pca_scores[i, 0], pca_scores[i, 1]),
        fontsize=7.5, ha='center', va='bottom', xytext=(0, 5),
        textcoords='offset points', color='#222'
    )

# 変数ベクトル（バイプロット）
components = pca.components_  # shape (2, n_vars)
scale = 3.0  # ベクトルスケール
for k, label in enumerate(ict_vars):
    vx = components[0, k] * scale
    vy = components[1, k] * scale
    ax2.annotate(
        '', xy=(vx, vy), xytext=(0, 0),
        arrowprops=dict(arrowstyle='->', color='#333', lw=1.8)
    )
    ax2.text(vx * 1.12, vy * 1.12, ict_labels[k],
             fontsize=9, color='#333', ha='center', va='center',
             bbox=dict(boxstyle='round,pad=0.2', fc='#ffffcc', ec='#999', alpha=0.85))

ax2.axhline(0, color='#bbb', lw=0.8, ls='--')
ax2.axvline(0, color='#bbb', lw=0.8, ls='--')
ax2.set_xlabel(f'PC1（寄与率 {var_expl[0]*100:.1f}%） ← ICT集積度', fontsize=12)
ax2.set_ylabel(f'PC2（寄与率 {var_expl[1]*100:.1f}%）', fontsize=12)
ax2.set_title(f'図2: PCAバイプロット（47都道府県・2022年）\n累積寄与率 {sum(var_expl)*100:.1f}%',
              fontsize=13, pad=12)
ax2.legend(loc='lower right', fontsize=10, framealpha=0.85)
fig2.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2021_U5_4_fig2.png'), bbox_inches='tight')
plt.close(fig2)
print("図2 保存完了")

# ══════════════════════════════════════════════════════════════════════════════
# 図3: 重回帰係数プロット（標準化偏回帰係数 + 95%CI）
# ══════════════════════════════════════════════════════════════════════════════
fig3, ax3 = plt.subplots(figsize=(8, 5))

colors_coef = ['#e05c5c' if p < 0.05 else '#aaa' for p in pvals]
y_pos = range(len(ict_labels))

for i, (coef, ci, color, label) in enumerate(zip(coef_vals, conf_int, colors_coef, ict_labels)):
    ax3.barh(i, coef, color=color, height=0.5, zorder=3)
    ax3.plot([ci[0], ci[1]], [i, i], 'k-', lw=2, zorder=4)
    ax3.plot([ci[0]], [i], 'k|', markersize=8, lw=2)
    ax3.plot([ci[1]], [i], 'k|', markersize=8, lw=2)

ax3.set_yticks(list(y_pos))
ax3.set_yticklabels([l.replace('\n', '') for l in ict_labels], fontsize=11)
ax3.axvline(0, color='black', lw=1.2)
ax3.set_xlabel('標準化偏回帰係数', fontsize=11)
ax3.set_title(f'図3: ICT指標が消費水準に与える影響\n標準化偏回帰係数と95%信頼区間  (R²={r2:.3f})', fontsize=12)
ax3.set_xlim(min(conf_int[:, 0]) * 1.3, max(conf_int[:, 1]) * 1.3)

# 有意性注記
from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor='#e05c5c', label='p < 0.05（有意）'),
    Patch(facecolor='#aaa',    label='p ≥ 0.05（非有意）')
]
ax3.legend(handles=legend_elements, fontsize=9, loc='lower right')
ax3.grid(axis='x', alpha=0.3, zorder=0)
fig3.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2021_U5_4_fig3.png'), bbox_inches='tight')
plt.close(fig3)
print("図3 保存完了")

# ══════════════════════════════════════════════════════════════════════════════
# 図4: ICTスコアランキング（上位10・下位10）
# ══════════════════════════════════════════════════════════════════════════════
top10_sorted = top10.sort_values('ICT指数', ascending=True)
bot10_sorted = bot10.sort_values('ICT指数', ascending=True)

fig4, axes = plt.subplots(1, 2, figsize=(13, 6), sharey=False)

# 上位10（赤系）
ax_top = axes[0]
colors_top = [region_colors[r] for r in top10_sorted['地域']]
bars = ax_top.barh(top10_sorted['都道府県'], top10_sorted['ICT指数'],
                   color=colors_top, edgecolor='white', height=0.7)
ax_top.axvline(0, color='black', lw=0.8)
ax_top.set_xlabel('ICT指数（PC1スコア）', fontsize=10)
ax_top.set_title('ICT先進県 上位10', fontsize=12, color='#c0392b', fontweight='bold')
for bar, val in zip(bars, top10_sorted['ICT指数']):
    ax_top.text(val + 0.05, bar.get_y() + bar.get_height() / 2,
                f'{val:.2f}', va='center', fontsize=9)
ax_top.tick_params(axis='y', labelsize=10)
ax_top.grid(axis='x', alpha=0.3)

# 下位10（青系）
ax_bot = axes[1]
bot10_rev = bot10.sort_values('ICT指数', ascending=False)
colors_bot = [region_colors[r] for r in bot10_rev['地域']]
bars2 = ax_bot.barh(bot10_rev['都道府県'], bot10_rev['ICT指数'],
                    color=colors_bot, edgecolor='white', height=0.7)
ax_bot.axvline(0, color='black', lw=0.8)
ax_bot.set_xlabel('ICT指数（PC1スコア）', fontsize=10)
ax_bot.set_title('ICT後発県 下位10', fontsize=12, color='#2980b9', fontweight='bold')
for bar, val in zip(bars2, bot10_rev['ICT指数']):
    offset = 0.05 if val >= 0 else -0.05
    ha = 'left' if val >= 0 else 'right'
    ax_bot.text(val + offset, bar.get_y() + bar.get_height() / 2,
                f'{val:.2f}', va='center', fontsize=9, ha=ha)
ax_bot.tick_params(axis='y', labelsize=10)
ax_bot.grid(axis='x', alpha=0.3)

# 凡例
from matplotlib.patches import Patch as MPatch
handles = [MPatch(facecolor=c, label=r) for r, c in region_colors.items()]
fig4.legend(handles=handles, loc='lower center', ncol=3, fontsize=9,
            bbox_to_anchor=(0.5, -0.04), framealpha=0.9)

fig4.suptitle('図4: ICT活用指数スコア ランキング（2022年・47都道府県）', fontsize=13, y=1.01)
fig4.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2021_U5_4_fig4.png'), bbox_inches='tight')
plt.close(fig4)
print("図4 保存完了")

# ── 最終サマリー出力 ─────────────────────────────────────────────────────────
print("\n" + "=" * 60)
print("最終統計サマリー（HTMLへの埋め込み用）")
print("=" * 60)
print(f"分析対象: 47都道府県（2022年断面）")
print(f"PC1寄与率: {var_expl[0]*100:.1f}%")
print(f"PC2寄与率: {var_expl[1]*100:.1f}%")
print(f"累積寄与率 (PC1+PC2): {sum(var_expl)*100:.1f}%")
print(f"重回帰 R²: {r2:.3f}")
print(f"重回帰 Adj-R²: {adj_r2:.3f}")

for i, (var, coef, pval) in enumerate(zip(ict_vars, coef_vals, pvals)):
    sig = '***' if pval < 0.001 else '**' if pval < 0.01 else '*' if pval < 0.05 else 'ns'
    print(f"  {var}: β={coef:.4f}, p={pval:.4f} {sig}")

print(f"\nICT先進県 TOP5: {', '.join(top10['都道府県'].head(5))}")
print(f"ICT後発県 BOT5: {', '.join(bot10['都道府県'].tail(5).tolist()[::-1])}")

# 相関係数（代表的なもの）
r_prog_cons, p_prog_cons = stats.pearsonr(df['大学進学率'], df['消費支出（二人以上の世帯）'])
r_univ_cons, p_univ_cons = stats.pearsonr(df['大学数per万人'], df['消費支出（二人以上の世帯）'])
r_comm_cons, p_comm_cons = stats.pearsonr(df['通信費比率'], df['消費支出（二人以上の世帯）'])
print(f"\n主要相関係数:")
print(f"  大学進学率 × 消費支出: r={r_prog_cons:.3f}, p={p_prog_cons:.4f}")
print(f"  大学数per万人 × 消費支出: r={r_univ_cons:.3f}, p={p_univ_cons:.4f}")
print(f"  通信費比率 × 消費支出: r={r_comm_cons:.3f}, p={p_comm_cons:.4f}")

print("\n全図の保存完了:")
for i in range(1, 5):
    path = os.path.join(FIG_DIR, f'2021_U5_4_fig{i}.png')
    print(f"  fig{i}: {os.path.exists(path)} → {path}")