"""
2019_U5_1_shorei.py
2019年度 統計データ分析コンペティション 審査員奨励賞（大学生部門）
「都道府県別労働生産性の決定要因：主成分分析とクラスタリングによる地域類型化」
手法: 主成分分析（PCA）、Ward法クラスタリング、重回帰分析、標準化係数
データ: SSDSE-B-2026.csv（47都道府県、2023年断面）
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・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/2019_U5_1_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 scipy.cluster.hierarchy import dendrogram, linkage, fcluster
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)

# ─────────────────────────────────────────
# 0. データ読み込み
# ─────────────────────────────────────────
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 = df_b[df_b['年度'] == 2023].copy().reset_index(drop=True)
print(f"2023年断面データ: {len(df)} 都道府県")

# ─────────────────────────────────────────
# 1. 変数構築（労働生産性の代理指標と説明変数）
# ─────────────────────────────────────────
# 【目的変数】労働生産性代理指標
# 消費支出（二人以上の世帯）は所得水準・生産性の代理変数として使用
# 求人倍率は労働市場の需給バランスを反映し、高生産性地域では高くなる傾向
# 2変数の主成分第1軸を「経済的豊かさ指標」として使用
df['消費支出_万'] = df['消費支出（二人以上の世帯）'] / 10000   # 万円単位
df['求人倍率'] = df['月間有効求人数（一般）'] / df['月間有効求職者数（一般）']

# 【説明変数①】教育水準
# 高校卒業後の大学進学率（知識資本の代理）
df['大学進学率'] = (
    df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数'] * 100
)
# 人口100人当たり大学学生数（知識集積の指標）
df['大学学生割合'] = df['大学学生数'] / df['総人口'] * 100

# 【説明変数②】人口構造
# 高齢化率（生産年齢人口の縮小→生産性低下要因）
df['高齢化率'] = df['65歳以上人口'] / df['総人口'] * 100
# 生産年齢人口比率（15〜64歳）
df['生産年齢率'] = df['15～64歳人口'] / df['総人口'] * 100

# 【説明変数③】人口移動（経済的魅力の指標）
# 転入超過率（人が集まる地域は経済的に活発）
df['転入超過率'] = (
    (df['転入者数（日本人移動者）'] - df['転出者数（日本人移動者）'])
    / df['総人口'] * 1000
)

# 【説明変数④】観光・サービス経済
# 人口100人当たり旅館宿泊施設数（観光産業の代理）
df['旅館密度'] = df['旅館営業施設数（ホテルを含む）'] / df['総人口'] * 10000

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

# ─────────────────────────────────────────
# 2. 目的変数：労働生産性スコア（消費支出と求人倍率の加重合成）
# ─────────────────────────────────────────
# 消費支出と求人倍率を標準化して単純平均 → 労働生産性代理スコア
z_cons = stats.zscore(df['消費支出_万'].values)
z_kyu  = stats.zscore(df['求人倍率'].values)
df['生産性スコア'] = (z_cons + z_kyu) / 2.0

print("\n=== 生産性スコア（上位5・下位5都道府県）===")
df_sorted = df.sort_values('生産性スコア', ascending=False)
print("上位5:", df_sorted.head(5)[['都道府県', '生産性スコア', '消費支出_万', '求人倍率']].to_string(index=False))
print("下位5:", df_sorted.tail(5)[['都道府県', '生産性スコア', '消費支出_万', '求人倍率']].to_string(index=False))

# ─────────────────────────────────────────
# 3. 主成分分析（PCA）
# ─────────────────────────────────────────
feature_cols = ['大学進学率', '大学学生割合', '高齢化率', '生産年齢率', '転入超過率', '旅館密度']
var_labels   = ['大学進学率', '大学学生割合', '高齢化率', '生産年齢率', '転入超過率', '旅館密度']

X = df[feature_cols].values
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

pca = PCA(n_components=2, random_state=42)
scores = pca.fit_transform(X_scaled)
loadings = pca.components_          # shape (2, 6)
explained_var = pca.explained_variance_ratio_

print("\n=== PCA 説明分散比 ===")
for i, ev in enumerate(explained_var):
    print(f"  PC{i+1}: {ev:.4f} ({ev*100:.1f}%)")
print(f"  PC1+PC2 累積: {(explained_var[0]+explained_var[1])*100:.1f}%")

# PC1 の符号調整: 大学進学率 loadings[0,0] が正→高教育=PC1高 になるよう
if loadings[0, 0] < 0:
    scores[:, 0] = -scores[:, 0]
    loadings[0, :] = -loadings[0, :]

df['PC1'] = scores[:, 0]
df['PC2'] = scores[:, 1]

# ─────────────────────────────────────────
# 4. Ward法クラスタリング（k=5）
# ─────────────────────────────────────────
Z_linkage = linkage(X_scaled, method='ward')
clusters = fcluster(Z_linkage, t=5, criterion='maxclust')
df['cluster'] = clusters

print("\n=== Ward法クラスタ（k=5）===")
cluster_labels = {}
for c in sorted(df['cluster'].unique()):
    prefs = df[df['cluster'] == c]['都道府県'].tolist()
    print(f"  クラスタ{c}: {prefs}")
    cluster_labels[c] = f'C{c}'

# ─────────────────────────────────────────
# 5. 重回帰分析（標準化係数）
# ─────────────────────────────────────────
# 目的変数: 消費支出（経済的豊かさの直接指標、N=47）
# 説明変数（単変量相関が高く独立性の高い4変数）:
#   - 大学進学率（r=+0.40）: 知識資本・人的資本
#   - 高齢化率（r=-0.31）:   高齢化による生産年齢縮小
#   - 転入超過率（r=+0.37）: 経済的魅力・集積効果
#   - 求人倍率（r=+0.00）:   労働市場の需給（統制変数）
# ※ 生産年齢率は高齢化率と-0.97相関のため除外（多重共線性回避）
y_prod = df['消費支出_万'].values

reg_cols   = ['大学進学率', '高齢化率', '転入超過率', '求人倍率']
reg_labels = ['大学進学率', '高齢化率', '転入超過率', '求人倍率']

X_reg_raw = df[reg_cols].values
X_reg_std = stats.zscore(X_reg_raw, axis=0)
y_std     = stats.zscore(y_prod)

X_reg_wc = sm.add_constant(X_reg_std)
model_std = sm.OLS(y_std, X_reg_wc).fit()

beta_coefs = model_std.params[1:]
beta_ses   = model_std.bse[1:]
beta_pvals = model_std.pvalues[1:]

print("\n=== 重回帰分析（標準化係数）===")
print(f"R2 = {model_std.rsquared:.4f}, R2adj = {model_std.rsquared_adj:.4f}")
print(f"F  = {model_std.fvalue:.4f}, p(F) = {model_std.f_pvalue:.4f}")
for label, b, se, p in zip(reg_labels, beta_coefs, beta_ses, beta_pvals):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {label}: β={b:.4f} (SE={se:.4f}, p={p:.4f}) {sig}")

# ─────────────────────────────────────────
# 図1: 労働生産性の都道府県間格差（水平棒グラフ）
# ─────────────────────────────────────────
df_plot = df.sort_values('生産性スコア').copy()

# 色分け（地域区分）
bar_colors = [region_colors[r] for r in df_plot['地域区分']]

fig, ax = plt.subplots(figsize=(9, 11))
bars = ax.barh(
    df_plot['都道府県'],
    df_plot['生産性スコア'],
    color=bar_colors,
    edgecolor='white',
    linewidth=0.5,
    alpha=0.88
)
ax.axvline(0, color='black', lw=0.8, ls='--', alpha=0.6)

# 凡例（地域区分）
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor=c, label=r, alpha=0.88)
                   for r, c in region_colors.items()]
ax.legend(handles=legend_elements, loc='lower right', fontsize=8,
          title='地域区分', title_fontsize=8, framealpha=0.9)

ax.set_xlabel('労働生産性スコア（標準化合成指標）', fontsize=11)
ax.set_title('図1：都道府県別 労働生産性スコア（2023年）\n（消費支出・求人倍率の標準化合成、地域色分け）',
             fontsize=12, fontweight='bold')
ax.grid(axis='x', alpha=0.3)
ax.tick_params(axis='y', labelsize=8)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2019_U5_1_fig1.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("\n[保存] 図1: 都道府県別労働生産性スコア")

# ─────────────────────────────────────────
# 図2: PCA バイプロット（第1・第2主成分）
# ─────────────────────────────────────────
fig, ax = plt.subplots(figsize=(10, 7.5))

for region, color in region_colors.items():
    mask = df['地域区分'] == region
    ax.scatter(
        scores[mask, 0], scores[mask, 1],
        c=color, label=region, s=60, zorder=3, alpha=0.85
    )

# 都道府県ラベル
for i, row in df.iterrows():
    ax.annotate(
        row['都道府県'],
        (scores[i, 0], scores[i, 1]),
        fontsize=6.5, ha='center', va='bottom',
        xytext=(0, 4), textcoords='offset points',
        color='#333333'
    )

# 変数負荷量ベクトル
scale = 2.8
for j, vl in enumerate(var_labels):
    lx = loadings[0, j] * scale
    ly = loadings[1, j] * scale
    ax.annotate(
        '',
        xy=(lx, ly), xytext=(0, 0),
        arrowprops=dict(arrowstyle='->', color='#C62828', lw=1.8)
    )
    ax.text(
        lx * 1.12, ly * 1.12, vl,
        fontsize=8.5, color='#C62828', ha='center', va='center',
        fontweight='bold'
    )

ax.axhline(0, color='gray', lw=0.6, ls='--', alpha=0.5)
ax.axvline(0, color='gray', lw=0.6, ls='--', alpha=0.5)
ax.set_xlabel(f'第1主成分（寄与率 {explained_var[0]*100:.1f}%）', fontsize=12)
ax.set_ylabel(f'第2主成分（寄与率 {explained_var[1]*100:.1f}%）', fontsize=12)
ax.set_title(
    '図2：PCA バイプロット（47都道府県、2023年）\n第1・第2主成分スコアと変数負荷量',
    fontsize=12, fontweight='bold'
)
ax.legend(title='地域区分', loc='lower right', fontsize=9, title_fontsize=9)
ax.grid(True, alpha=0.3)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2019_U5_1_fig2.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("[保存] 図2: PCA バイプロット")

# ─────────────────────────────────────────
# 図3: Ward法デンドログラム（k=5クラスター）
# ─────────────────────────────────────────
fig, ax = plt.subplots(figsize=(14, 6))

sorted_z = np.sort(Z_linkage[:, 2])
# k=5 のカット高さ: k個のクラスターを得るには最後から k-1 番目の結合の間でカット
cut_height = (sorted_z[-4] + sorted_z[-5]) / 2.0  # 5クラスターの境界

dendrogram(
    Z_linkage,
    labels=list(df['都道府県']),
    ax=ax,
    leaf_rotation=90,
    leaf_font_size=8.5,
    color_threshold=cut_height
)
ax.axhline(
    cut_height, color='red', ls='--', lw=1.8,
    label=f'k=5 カット（高さ≈{cut_height:.2f}）'
)
ax.set_title(
    '図3：Ward法クラスタリング デンドログラム（47都道府県、2023年）\n赤破線でk=5クラスターに分割',
    fontsize=12, fontweight='bold'
)
ax.set_xlabel('都道府県', fontsize=10)
ax.set_ylabel('Ward距離', fontsize=10)
ax.legend(fontsize=10)
ax.tick_params(axis='x', labelsize=8)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2019_U5_1_fig3.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("[保存] 図3: Ward法デンドログラム")

# ─────────────────────────────────────────
# 図4: 重回帰分析の標準化係数（棒グラフ + エラーバー + 有意水準）
# ─────────────────────────────────────────
def sig_stars(p):
    if p < 0.001:
        return '***'
    elif p < 0.01:
        return '**'
    elif p < 0.05:
        return '*'
    else:
        return 'n.s.'

bar_colors_reg = ['#1565C0' if b >= 0 else '#C62828' for b in beta_coefs]

fig, ax = plt.subplots(figsize=(8, 5.5))
x_pos = np.arange(len(reg_labels))
ax.bar(x_pos, beta_coefs, color=bar_colors_reg, alpha=0.82,
       width=0.55, edgecolor='white', linewidth=0.8)
ax.errorbar(x_pos, beta_coefs, yerr=beta_ses,
            fmt='none', ecolor='#333333', capsize=6, capthick=1.8, elinewidth=1.8)

# 有意水準ラベル
for i, (b, p) in enumerate(zip(beta_coefs, beta_pvals)):
    stars = sig_stars(p)
    offset = 0.06 if b >= 0 else -0.06
    va = 'bottom' if b >= 0 else 'top'
    ax.text(i, b + offset, stars, ha='center', va=va, fontsize=11,
            fontweight='bold', color='#333333')

ax.axhline(0, color='black', lw=0.8, ls='--')
ax.set_xticks(x_pos)
ax.set_xticklabels(reg_labels, fontsize=11)
ax.set_ylabel('標準化偏回帰係数（β）', fontsize=11)
ax.set_title(
    '図4：重回帰分析 標準化係数（β係数）\n目的変数：労働生産性スコア（2023年、N=47）',
    fontsize=12, fontweight='bold'
)
ax.set_ylim(
    min(beta_coefs) - max(beta_ses) * 2.5,
    max(beta_coefs) + max(beta_ses) * 2.5
)
ax.grid(axis='y', alpha=0.3)

# 凡例（有意水準）
from matplotlib.lines import Line2D
legend_items = [
    Line2D([0], [0], color='none', label='*** p < 0.001'),
    Line2D([0], [0], color='none', label='**  p < 0.01'),
    Line2D([0], [0], color='none', label='*   p < 0.05'),
    Line2D([0], [0], color='none', label='n.s. p >= 0.05'),
]
ax.legend(handles=legend_items, fontsize=9, loc='upper right', framealpha=0.9)

# R2 表示
ax.text(0.02, 0.97,
        f"R² = {model_std.rsquared:.3f}  (adj.R² = {model_std.rsquared_adj:.3f})\n"
        f"F({int(model_std.df_model)}, {int(model_std.df_resid)}) = {model_std.fvalue:.2f}, "
        f"p = {model_std.f_pvalue:.4f}",
        transform=ax.transAxes, fontsize=9, va='top',
        bbox=dict(boxstyle='round,pad=0.4', facecolor='#EFF3FF', alpha=0.9))

fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2019_U5_1_fig4.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("[保存] 図4: 標準化係数棒グラフ")

# ─────────────────────────────────────────
# 結果サマリー出力
# ─────────────────────────────────────────
print("\n====== 分析結果サマリー ======")
print(f"データ: SSDSE-B-2026.csv 2023年断面, N=47都道府県")
print(f"\nPCA 説明分散比:")
print(f"  PC1: {explained_var[0]*100:.1f}%  PC2: {explained_var[1]*100:.1f}%  累積: {sum(explained_var)*100:.1f}%")
print(f"\nPCA 負荷量（PC1）:")
for col, load in zip(feature_cols, loadings[0]):
    print(f"  {col}: {load:.4f}")
print(f"\nクラスター（k=5）構成:")
for c in sorted(df['cluster'].unique()):
    prefs = df[df['cluster'] == c]['都道府県'].tolist()
    print(f"  C{c}: {prefs}")
print(f"\n重回帰 R²={model_std.rsquared:.4f}  adj.R²={model_std.rsquared_adj:.4f}")
print(f"標準化係数:")
for label, b, p in zip(reg_labels, beta_coefs, beta_pvals):
    print(f"  {label}: β={b:.4f} {sig_stars(p)}")

print("\nDONE: 2019_U5_1_shorei")