都道府県別労働生産性の決定要因
主成分分析とクラスタリングによる地域類型化

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)

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))

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}'

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")

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: 都道府県別労働生産性スコア")

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 バイプロット")

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法デンドログラム")