地域産業構造と経済成長
特化係数・主成分分析による産業集積の分析

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

PREFIX = '2021_U3'

def save_fig(label):
    path = os.path.join(FIG_DIR, f'{PREFIX}_{label}.png')
    plt.savefig(path, bbox_inches='tight', dpi=150)
    plt.close()
    print(f'  → {path} 保存完了')

print('=== データ読み込み ===')
df_raw = pd.read_csv(DATA_B, encoding='cp932', header=1)

# 都道府県行のみ抽出（地域コードが R+5桁、かつ全国コードR00000を除外）
mask = (df_raw['地域コード'].str.match(r'^R\d{5}$', na=False) &
        (df_raw['地域コード'] != 'R00000'))
df_all = df_raw[mask].copy()

# 2022年データ（消費支出が揃っている最新年）
df = df_all[df_all['年度'] == 2022].copy().reset_index(drop=True)
print(f'2022年データ: {len(df)}都道府県')

# 数値変換
num_cols = ['総人口', '旅館営業施設数（ホテルを含む）', '大学学生数',
            '一般病院数', '一般診療所数', '着工新設住宅戸数',
            '1人1日当たりの排出量', '消費支出（二人以上の世帯）']
for c in num_cols:
    df[c] = pd.to_numeric(df[c], errors='coerce')

# 欠損値確認
print('欠損値確認:')
for c in num_cols:
    n_na = df[c].isna().sum()
    if n_na > 0:
        print(f'  {c}: {n_na}件')

print('\n=== 特化係数代理指標の計算 ===')

pop = df['総人口']

# 観光特化: 旅館施設数 / 総人口 × 10000
df['観光特化'] = df['旅館営業施設数（ホテルを含む）'] / pop * 10000

# 農業特化代理: ごみ排出量の逆数（農村部はごみが少ない）
df['農業特化代理'] = 1000 / df['1人1日当たりの排出量']

# 教育特化: 大学学生数 / 総人口 × 1000
df['教育特化'] = df['大学学生数'] / pop * 1000

# 医療特化: (病院数+診療所数) / 総人口 × 10000
df['医療特化'] = (df['一般病院数'] + df['一般診療所数']) / pop * 10000

# 建設特化代理: 新設住宅戸数 / 総人口 × 10000
df['建設特化代理'] = df['着工新設住宅戸数'] / pop * 10000

# 消費支出（log変換）
df['消費支出_log'] = np.log(df['消費支出（二人以上の世帯）'])

# 地域ブロック
df['地域'] = df['都道府県'].map(REGION_MAP)

# 変数リスト
VARS = ['観光特化', '農業特化代理', '教育特化', '医療特化', '建設特化代理', '消費支出_log']
VAR_LABELS = ['観光\n特化', '農業\n特化代理', '教育\n特化', '医療\n特化', '建設\n特化代理', '消費支出\n(log)']

# 欠損除去
df_clean = df[['都道府県', '地域'] + VARS + ['消費支出（二人以上の世帯）']].dropna().reset_index(drop=True)
print(f'分析対象: {len(df_clean)}都道府県')

# 統計サマリー
print('\n--- 代理指標サマリー ---')
print(df_clean[VARS].describe().round(4).to_string())

lq_df = df_clean[VARS].copy()
for v in VARS:
    nat_mean = lq_df[v].mean()
    lq_df[v] = lq_df[v] / nat_mean  # LQ: 地域値 / 全国平均

lq_df['地域'] = df_clean['地域'].values
lq_df['都道府県'] = df_clean['都道府県'].values

print('\n--- LQ値サマリー（全国平均=1.0）---')
print(lq_df[VARS].describe().round(3).to_string())

print('\n=== Fig1: レーダーチャート ===')

region_lq = lq_df.groupby('地域')[VARS].mean()

N_vars = len(VARS)
angles = np.linspace(0, 2 * np.pi, N_vars, endpoint=False).tolist()
angles += angles[:1]  # close the polygon

fig, axes = plt.subplots(2, 3, figsize=(14, 9),
                         subplot_kw=dict(polar=True))
fig.suptitle('地域ブロック別 産業特化指標（LQ代理）レーダーチャート\n（2022年, 全国平均=1.0）',
             fontsize=14, fontweight='bold', y=1.01)

axes_flat = axes.flatten()

for i, region in enumerate(REGION_ORDER[:6]):
    ax = axes_flat[i]
    if region not in region_lq.index:
        ax.set_visible(False)
        continue
    vals = region_lq.loc[region, VARS].tolist()
    vals += vals[:1]

    color = REGION_COLORS[region]
    ax.plot(angles, vals, color=color, linewidth=2.0)
    ax.fill(angles, vals, color=color, alpha=0.25)

    # 参照線（LQ=1.0）
    ax.plot(angles, [1.0] * len(angles), color='gray',
            linewidth=0.8, linestyle='--', alpha=0.6)

    ax.set_xticks(angles[:-1])
    ax.set_xticklabels(VAR_LABELS, fontsize=8)
    ax.set_ylim(0, None)
    ax.set_title(region, fontsize=11, fontweight='bold',
                 color=color, pad=12)
    ax.tick_params(axis='y', labelsize=7)

plt.tight_layout()
save_fig('fig1_radar')

print('\n=== Fig2: PCAスコアプロット ===')

fig, ax = plt.subplots(figsize=(11, 8))

for region in REGION_ORDER:
    mask_r = df_clean['地域'].values == region
    if mask_r.sum() == 0:
        continue
    ax.scatter(scores[mask_r, 0], scores[mask_r, 1],
               color=REGION_COLORS[region], s=80, alpha=0.85,
               label=region, zorder=3)
    for idx in np.where(mask_r)[0]:
        pref = df_clean['都道府県'].iloc[idx]
        ax.annotate(pref, (scores[idx, 0], scores[idx, 1]),
                    fontsize=7, xytext=(4, 3),
                    textcoords='offset points',
                    color=REGION_COLORS[region])

# 原点に十字線
ax.axhline(0, color='gray', linewidth=0.6, linestyle='--', alpha=0.5)
ax.axvline(0, color='gray', linewidth=0.6, linestyle='--', alpha=0.5)

ax.set_xlabel(f'PC1（寄与率 {ev_ratio[0]*100:.1f}%）', fontsize=12)
ax.set_ylabel(f'PC2（寄与率 {ev_ratio[1]*100:.1f}%）', fontsize=12)
ax.set_title('主成分分析（PCA）スコアプロット\n産業特化指標の都道府県別分布（2022年）',
             fontsize=13, fontweight='bold')
ax.legend(title='地域ブロック', bbox_to_anchor=(1.01, 1), loc='upper left',
          fontsize=9, title_fontsize=9)

# 因子負荷量の矢印（biplot）
scale_factor = 2.5
for j, v in enumerate(VARS):
    lx, ly = loadings[0, j] * scale_factor, loadings[1, j] * scale_factor
    ax.annotate('', xy=(lx, ly), xytext=(0, 0),
                arrowprops=dict(arrowstyle='->', color='#E53935', lw=1.5))
    ax.text(lx * 1.08, ly * 1.08, v.replace('代理', '代理\n').replace('特化', '特化\n'),
            fontsize=8, color='#E53935', ha='center')

ax.grid(True, alpha=0.3)
plt.tight_layout()
save_fig('fig2_pca')

print('\n=== Fig3: Ward法デンドログラム ===')

linked = linkage(X_scaled, method='ward')

fig, ax = plt.subplots(figsize=(14, 6))

dend = dendrogram(
    linked,
    labels=df_clean['都道府県'].values,
    ax=ax,
    leaf_rotation=90,
    leaf_font_size=8,
    color_threshold=0.7 * max(linked[:, 2]),
)

ax.set_title('Ward法階層クラスタリング デンドログラム\n産業特化指標による都道府県類型（2022年）',
             fontsize=13, fontweight='bold')
ax.set_xlabel('都道府県', fontsize=11)
ax.set_ylabel('非類似度（Ward距離）', fontsize=11)

# k=4 でカットする水準を示す補助線
# Ward距離の上位ギャップを探してカット点を示す
sorted_dist = np.sort(linked[:, 2])
gaps = np.diff(sorted_dist)
# 上位5ギャップのうち最大を採用（k=4カット）
n_clusters_target = 4
cutoff_height = (linked[-n_clusters_target + 1, 2] + linked[-n_clusters_target, 2]) / 2
ax.axhline(y=cutoff_height, color='red', linestyle='--', linewidth=1.5,
           label=f'k={n_clusters_target} クラスタ分割線')
ax.legend(fontsize=9)

plt.tight_layout()
save_fig('fig3_dendro')

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from scipy.cluster.hierarchy import dendrogram, linkage
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

import os
FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

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

REGION_COLORS = {
    '北海道':     '#E53935',
    '東北':       '#FB8C00',
    '関東':       '#43A047',
    '中部':       '#1E88E5',
    '近畿':       '#8E24AA',
    '中国・四国': '#00ACC1',
    '九州・沖縄': '#6D4C41',
}

REGION_ORDER = ['北海道', '東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄']

print('\n=== PCA ===')

X_raw = df_clean[VARS].values
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_raw)

pca = PCA(n_components=min(len(VARS), len(df_clean)))
scores = pca.fit_transform(X_scaled)
loadings = pca.components_  # shape: (n_components, n_vars)
ev_ratio = pca.explained_variance_ratio_

print('寄与率:')
for i, ev in enumerate(ev_ratio[:4]):
    print(f'  PC{i+1}: {ev*100:.1f}%  (累積: {ev_ratio[:i+1].sum()*100:.1f}%)')

# 因子負荷量
print('\n--- 因子負荷量（PC1, PC2）---')
for j, v in enumerate(VARS):
    print(f'  {v:12s}: PC1={loadings[0,j]:+.3f}, PC2={loadings[1,j]:+.3f}')

print('\n=== Fig4: 観光特化 vs 消費支出 散布図 ===')

x_plot = df_clean['観光特化'].values
y_plot = df_clean['消費支出（二人以上の世帯）'].values / 1000  # 千円単位

# 相関係数
r_val, p_val = stats.pearsonr(x_plot, y_plot)
print(f'相関係数: r={r_val:.4f}, p={p_val:.4f}')

# 回帰直線
slope, intercept, _, _, stderr = stats.linregress(x_plot, y_plot)
x_fit = np.linspace(x_plot.min(), x_plot.max(), 100)
y_fit = slope * x_fit + intercept

fig, ax = plt.subplots(figsize=(10, 7))

for region in REGION_ORDER:
    mask_r = df_clean['地域'].values == region
    if mask_r.sum() == 0:
        continue
    ax.scatter(x_plot[mask_r], y_plot[mask_r],
               color=REGION_COLORS[region], s=80, alpha=0.85,
               label=region, zorder=3)
    for idx in np.where(mask_r)[0]:
        pref = df_clean['都道府県'].iloc[idx]
        ax.annotate(pref, (x_plot[idx], y_plot[idx]),
                    fontsize=7, xytext=(4, 3),
                    textcoords='offset points',
                    color=REGION_COLORS[region])

# 回帰直線
line_style = '-' if p_val < 0.05 else '--'
ax.plot(x_fit, y_fit, color='black', linewidth=1.8,
        linestyle=line_style, label=f'回帰直線 (r={r_val:.3f}, p={p_val:.3f})')

ax.set_xlabel('観光特化指標（旅館施設数/万人）', fontsize=12)
ax.set_ylabel('消費支出（千円, 二人以上世帯）', fontsize=12)
ax.set_title('観光産業集積と消費支出水準の関係（2022年）\n都道府県別クロスセクション分析',
             fontsize=13, fontweight='bold')
ax.legend(title='地域ブロック', bbox_to_anchor=(1.01, 1), loc='upper left',
          fontsize=9, title_fontsize=9)
ax.grid(True, alpha=0.3)

# 相関注釈
sig_str = '**' if p_val < 0.01 else ('*' if p_val < 0.05 else 'n.s.')
ax.text(0.05, 0.95,
        f'r = {r_val:.3f} {sig_str}\np = {p_val:.3f}\nN = {len(df_clean)}',
        transform=ax.transAxes, fontsize=10,
        verticalalignment='top',
        bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))

plt.tight_layout()
save_fig('fig4_scatter')

print('\n=== 全図の保存完了 ===')
print('Fig1: 産業特化指標レーダーチャート（地域ブロック別）')
print('Fig2: PCAスコアプロット（都道府県ラベル付き biplot）')
print('Fig3: Ward法クラスタリング デンドログラム')
print('Fig4: 観光特化 vs 消費支出 散布図（2022年）')