"""
教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [大学生・一般の部]
=================================================================
論文タイトル：デジタルデバイドの地域格差：交通・通信費を指標とした情報格差の分析
受賞：審査員奨励賞（大学生・一般の部）

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県別パネルデータ, 2012〜2023年度）
  対象：全47都道府県 × 12年（2012〜2023）

  代理変数：交通・通信費（二人以上の世帯）→ IT・情報インフラへの支出の代理指標

  分析の流れ
  1. 時系列：交通・通信費の地域別推移
  2. OLS回帰：交通・通信費の決定要因（高齢化率・大学数・消費支出・人口密度）
  3. PCA：情報格差関連指標の主成分分析
  4. Ward法クラスタリング：都道府県のデジタルデバイドグループ分類

【目的変数（IT格差代理）】
  交通・通信費（二人以上の世帯）= IT・通信インフラ支出の近似

【説明変数】
  高齢化率 = 65歳以上人口 / 総人口 × 100  （高齢者↑ → IT利用↓）
  大学数              （高等教育機関の集積 → IT親和性↑）
  消費支出_log       （所得水準の代理）
  保育所密度          （若い世帯の多さ → IT利用↑）

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）

【データサイエンス学習ポイント】
  1. 代理変数の設計（直接測定できない「IT格差」の近似方法）
  2. 主成分分析（PCA）による多次元指標の次元削減
  3. Ward法クラスタリングの手順と解釈
  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/2022_U5_15_shorei.py
# ============================================================


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from scipy.cluster.hierarchy import dendrogram, linkage
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
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

import os
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)
df_b = df_b.sort_values(['都道府県', '年度']).reset_index(drop=True)

# 変数生成
df_b['通信費_千円'] = df_b['交通・通信費（二人以上の世帯）'] / 1000
df_b['通信費_log'] = np.log(df_b['交通・通信費（二人以上の世帯）'].clip(lower=1))
df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100
df_b['消費支出_log'] = np.log(df_b['消費支出（二人以上の世帯）'].clip(lower=1))
df_b['保育所密度'] = df_b['保育所等数'] / df_b['総人口'] * 10000
df_b['教育費_千円'] = df_b['教育費（二人以上の世帯）'] / 1000

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

# Fig1: 交通・通信費の時系列推移（地域別）
fig, ax = plt.subplots(figsize=(10, 5))
yearly = df_b.groupby(['年度', '地域'])['通信費_千円'].mean().reset_index()
for reg, grp in yearly.groupby('地域'):
    ax.plot(grp['年度'], grp['通信費_千円'], marker='o', markersize=4,
            label=reg, color=region_colors.get(reg, 'gray'))
ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('交通・通信費（千円/月）', fontsize=12)
ax.set_title('地域別 交通・通信費の推移（2012〜2023年）\n（IT格差の代理指標）', fontsize=14, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_15_fig1_ts.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig1 saved")

# Fig2: OLS回帰 係数プロット
df_2022 = df_b[df_b['年度'] == 2022].copy()
xvars = ['高齢化率', '大学数', '消費支出_log', '保育所密度']
df_reg = df_2022[['通信費_log'] + xvars].dropna()
X = sm.add_constant(df_reg[xvars])
res = sm.OLS(df_reg['通信費_log'], X).fit()
print(res.summary())
coefs = res.params.drop('const')
ses = res.bse.drop('const')
pvals = res.pvalues.drop('const')
fig, ax = plt.subplots(figsize=(8, 5))
colors_c = ['#e05c5c' if p < 0.05 else '#888888' for p in pvals]
ax.barh(range(len(coefs)), coefs, xerr=1.96 * ses, color=colors_c, alpha=0.8,
        error_kw={'elinewidth': 1.5, 'capsize': 4})
ax.set_yticks(range(len(coefs)))
ax.set_yticklabels(coefs.index, fontsize=10)
ax.axvline(0, color='black', linewidth=0.8)
ax.set_xlabel('OLS回帰係数', fontsize=12)
ax.set_title(f'交通・通信費の決定要因（OLS, N=47）\nR²={res.rsquared:.3f}（赤=p<0.05）', fontsize=12, fontweight='bold')
ax.grid(axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_15_fig2_ols.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig2 saved")

# Fig3: PCA 主成分分析
pca_vars = ['通信費_千円', '高齢化率', '大学数', '消費支出_log', '教育費_千円', '保育所密度']
df_pca = df_2022[['都道府県', '地域'] + pca_vars].dropna()
scaler = StandardScaler()
X_sc = scaler.fit_transform(df_pca[pca_vars])
pca = PCA(n_components=2)
scores = pca.fit_transform(X_sc)
explained = pca.explained_variance_ratio_

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
# スコアプロット
for reg in df_pca['地域'].unique():
    mask = df_pca['地域'] == reg
    ax1.scatter(scores[mask, 0], scores[mask, 1],
                color=region_colors.get(reg, 'gray'), label=reg, s=60, alpha=0.8)
for i, pref in enumerate(df_pca['都道府県']):
    ax1.annotate(pref[:2], (scores[i, 0], scores[i, 1]),
                 fontsize=6, alpha=0.6, xytext=(2, 2), textcoords='offset points')
ax1.set_xlabel(f'PC1 ({explained[0]*100:.1f}%)', fontsize=11)
ax1.set_ylabel(f'PC2 ({explained[1]*100:.1f}%)', fontsize=11)
ax1.set_title('PCAスコアプロット（2022年, N=47）', fontsize=12, fontweight='bold')
ax1.legend(fontsize=8, ncol=2)
ax1.grid(alpha=0.3)
ax1.axhline(0, color='black', linewidth=0.5)
ax1.axvline(0, color='black', linewidth=0.5)
# 寄与率
ax2.bar(range(1, len(pca_vars) + 1),
        PCA(n_components=len(pca_vars)).fit(X_sc).explained_variance_ratio_ * 100,
        color='#4e9af1', alpha=0.8)
ax2.set_xlabel('主成分番号', fontsize=11)
ax2.set_ylabel('寄与率（%）', fontsize=11)
ax2.set_title('主成分の寄与率（スクリープロット）', fontsize=12, fontweight='bold')
ax2.grid(axis='y', alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_15_fig3_pca.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig3 saved")

# Fig4: Ward法クラスタリング デンドログラム
from scipy.cluster.hierarchy import dendrogram, linkage
Z = linkage(X_sc, method='ward')
fig, ax = plt.subplots(figsize=(12, 6))
dendrogram(Z, labels=df_pca['都道府県'].tolist(), ax=ax,
           color_threshold=Z[-3, 2], leaf_rotation=90, leaf_font_size=8)
ax.set_title('Ward法クラスタリング：デジタルデバイド指標による都道府県分類（2022年）',
             fontsize=13, fontweight='bold')
ax.set_xlabel('都道府県', fontsize=11)
ax.set_ylabel('距離（Ward法）', fontsize=11)
ax.grid(axis='y', alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_15_fig4_cluster.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig4 saved")
print("All done!")