製造業・産業構造と地域雇用・所得の関係分析

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

print("SSDSE-E を読み込み中...")
df_e_raw = pd.read_csv(DATA_E, encoding='cp932', header=0)

# 行0: 年度, 行1: 変数名, 行2以降: データ
# 都道府県データ（R01000〜R47000）と全国（R00000）を分離
df_e_all = df_e_raw[df_e_raw['SSDSE-E-2026'] == 'R00000'].copy()  # 全国
df_e_pref = df_e_raw[
    df_e_raw['SSDSE-E-2026'].str.match(r'^R\d{5}$', na=False) &
    (df_e_raw['SSDSE-E-2026'] != 'R00000')
].copy()

# 数値変換対象列
e_num_cols = [
    'C2208',    # 従業者数（民営）計
    'C220836',  # 建設業
    'C220837',  # 製造業
    'C220839',  # 情報通信業
    'C220841',  # 卸売業・小売業
    'C220847',  # 宿泊業・飲食サービス業
    'C220848',  # 生活関連サービス業・娯楽業
    'C220850',  # 医療・福祉
    'C122101',  # 1人当たり県民所得（万円）
    'A1101',    # 総人口
    'A1303',    # 65歳以上人口
    'B1101',    # 総面積
]
for col in e_num_cols:
    df_e_pref[col] = pd.to_numeric(df_e_pref[col], errors='coerce')
    df_e_all[col]  = pd.to_numeric(df_e_all[col],  errors='coerce')

# 都道府県名を列として整理
df_e_pref = df_e_pref.rename(columns={'SSDSE-E-2026': '地域コード', 'Prefecture': '都道府県'})
df_e_pref = df_e_pref.reset_index(drop=True)

total_workers = df_e_pref['C2208']

df_e_pref['従業者割合_製造業']       = df_e_pref['C220837'] / total_workers * 100
df_e_pref['従業者割合_建設業']       = df_e_pref['C220836'] / total_workers * 100
df_e_pref['従業者割合_情報通信業']   = df_e_pref['C220839'] / total_workers * 100
df_e_pref['従業者割合_卸売小売業']   = df_e_pref['C220841'] / total_workers * 100
df_e_pref['従業者割合_宿泊飲食業']   = df_e_pref['C220847'] / total_workers * 100
df_e_pref['従業者割合_生活娯楽業']   = df_e_pref['C220848'] / total_workers * 100
df_e_pref['従業者割合_医療福祉業']   = df_e_pref['C220850'] / total_workers * 100

nat_total = float(df_e_all['C2208'].values[0])
nat_mfg   = float(df_e_all['C220837'].values[0])
nat_const = float(df_e_all['C220836'].values[0])
nat_ict   = float(df_e_all['C220839'].values[0])
nat_retail= float(df_e_all['C220841'].values[0])
nat_med   = float(df_e_all['C220850'].values[0])

nat_mfg_ratio   = nat_mfg   / nat_total
nat_const_ratio = nat_const / nat_total
nat_ict_ratio   = nat_ict   / nat_total
nat_retail_ratio= nat_retail/ nat_total
nat_med_ratio   = nat_med   / nat_total

df_e_pref['LQ_製造業']     = (df_e_pref['C220837'] / total_workers) / nat_mfg_ratio
df_e_pref['LQ_建設業']     = (df_e_pref['C220836'] / total_workers) / nat_const_ratio
df_e_pref['LQ_情報通信業'] = (df_e_pref['C220839'] / total_workers) / nat_ict_ratio
df_e_pref['LQ_卸売小売業'] = (df_e_pref['C220841'] / total_workers) / nat_retail_ratio
df_e_pref['LQ_医療福祉業'] = (df_e_pref['C220850'] / total_workers) / nat_med_ratio

# 高齢化率・人口密度
df_e_pref['高齢化率'] = df_e_pref['A1303'] / df_e_pref['A1101'] * 100
df_e_pref['人口密度'] = df_e_pref['A1101'] / df_e_pref['B1101']  # 人/km²

print(f"SSDSE-E: {len(df_e_pref)} 都道府県")

print("SSDSE-B を読み込み中...")
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_b2022 = df_b[df_b['年度'] == 2022].copy()

# 数値変換
b_num_cols = [
    '消費支出（二人以上の世帯）',
    '月間有効求人数（一般）',
    '就職件数（一般）',
    '転入者数（日本人移動者）',
    '転出者数（日本人移動者）',
    '総人口',
    '65歳以上人口',
    '15～64歳人口',
]
for col in b_num_cols:
    if col in df_b2022.columns:
        df_b2022[col] = pd.to_numeric(df_b2022[col], errors='coerce')

# 有効求人倍率（月間有効求人数 / 就職件数（分母として近似））
# 実際はハローワーク統計で算出されるが、代理指標として使用
df_b2022['有効求人倍率_代理'] = (
    df_b2022['月間有効求人数（一般）'] / df_b2022['就職件数（一般）']
)

# 転入率（転入/総人口×1000）
df_b2022['転入率'] = df_b2022['転入者数（日本人移動者）'] / df_b2022['総人口'] * 1000

print(f"SSDSE-B 2022年: {len(df_b2022)} 都道府県")

df = pd.merge(
    df_e_pref[['都道府県', '地域コード',
               '従業者割合_製造業', '従業者割合_建設業', '従業者割合_情報通信業',
               '従業者割合_卸売小売業', '従業者割合_宿泊飲食業', '従業者割合_医療福祉業',
               'LQ_製造業', 'LQ_建設業', 'LQ_情報通信業', 'LQ_卸売小売業', 'LQ_医療福祉業',
               'C122101', '高齢化率', '人口密度']],
    df_b2022[['都道府県', '消費支出（二人以上の世帯）',
              '有効求人倍率_代理', '転入率']],
    on='都道府県',
    how='inner'
)

df['地域区分'] = df['都道府県'].map(region_map)
df['地域色']  = df['地域区分'].map(region_colors)

# 1人当たり県民所得を万円単位に（元データは千円単位）
df['県民所得_万円'] = df['C122101'] / 10  # 千円→万円

print(f"結合後: {len(df)} 都道府県")
print()

print()
print("図1 生成中...")

# 製造業割合で上位10・下位10を選択
df_sorted_mfg = df.sort_values('従業者割合_製造業', ascending=False)
target_prefs = pd.concat([df_sorted_mfg.head(10), df_sorted_mfg.tail(10)])
target_prefs = target_prefs.reset_index(drop=True)

# 積み上げ棒グラフの業種
bar_cols = [
    ('従業者割合_製造業',     '製造業',     '#e05c5c'),
    ('従業者割合_建設業',     '建設業',     '#f0a500'),
    ('従業者割合_情報通信業', '情報通信業', '#4e9af1'),
    ('従業者割合_卸売小売業', '卸売・小売業', '#5cb85c'),
    ('従業者割合_宿泊飲食業', '宿泊・飲食業', '#9b59b6'),
    ('従業者割合_医療福祉業', '医療・福祉業', '#f39c12'),
]

# 残りをその他として計算
listed_cols = [c[0] for c in bar_cols]
target_prefs['その他'] = 100 - target_prefs[listed_cols].sum(axis=1)

fig1, ax1 = plt.subplots(figsize=(14, 7))

x = np.arange(len(target_prefs))
bottom = np.zeros(len(target_prefs))

for col, label, color in bar_cols:
    vals = target_prefs[col].values
    ax1.bar(x, vals, bottom=bottom, label=label, color=color, width=0.72, edgecolor='white', linewidth=0.5)
    bottom += vals

# その他
ax1.bar(x, target_prefs['その他'].values, bottom=bottom,
        label='その他', color='#aaaaaa', width=0.72, edgecolor='white', linewidth=0.5)

# 上位10と下位10の境界線
ax1.axvline(x=9.5, color='black', linestyle='--', linewidth=1.2, alpha=0.6)
ax1.text(4.5, 102, '製造業割合 上位10県', ha='center', fontsize=11, fontweight='bold', color='#333333')
ax1.text(14.5, 102, '製造業割合 下位10県', ha='center', fontsize=11, fontweight='bold', color='#666666')

# ラベル設定
ax1.set_xticks(x)
pref_labels = target_prefs['都道府県'].str.replace('県', '').str.replace('府', '').str.replace('都', '').str.replace('道', '')
ax1.set_xticklabels(pref_labels, rotation=45, ha='right', fontsize=9)
ax1.set_ylabel('産業別従業者割合（%）', fontsize=12)
ax1.set_title('都道府県別産業構造（製造業割合 上位・下位各10都道府県）\n─ 民営事業所 従業者数ベース（2021年） ─',
              fontsize=13, fontweight='bold', pad=12)
ax1.legend(loc='upper right', fontsize=9, ncol=2)
ax1.set_ylim(0, 115)
ax1.set_xlim(-0.6, len(target_prefs) - 0.4)
ax1.yaxis.grid(True, alpha=0.3, linestyle='--')
ax1.set_axisbelow(True)

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2020_H5_3_fig1.png'), bbox_inches='tight')
plt.close(fig1)
print("  -> 2020_H5_3_fig1.png 保存完了")

print("図2 生成中...")

fig2, ax2 = plt.subplots(figsize=(10, 7))

for region, color in region_colors.items():
    mask = df['地域区分'] == region
    ax2.scatter(
        df.loc[mask, '従業者割合_製造業'],
        df.loc[mask, '消費支出（二人以上の世帯）'],
        c=color, label=region, s=70, alpha=0.85, edgecolors='white', linewidth=0.6, zorder=3
    )

# 都道府県ラベル
for _, row in df.iterrows():
    pref_short = row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', '')
    ax2.annotate(
        pref_short,
        xy=(row['従業者割合_製造業'], row['消費支出（二人以上の世帯）']),
        xytext=(3, 3), textcoords='offset points',
        fontsize=7, color='#444444', zorder=4
    )

# 回帰直線
x_mfg = df['従業者割合_製造業'].values
y_csp  = df['消費支出（二人以上の世帯）'].values
mask_valid = ~(np.isnan(x_mfg) | np.isnan(y_csp))
slope, intercept, r_val, p_val, std_err = stats.linregress(x_mfg[mask_valid], y_csp[mask_valid])
x_line = np.linspace(x_mfg[mask_valid].min(), x_mfg[mask_valid].max(), 100)
ax2.plot(x_line, slope * x_line + intercept, 'k--', linewidth=1.5, alpha=0.7, zorder=2, label=f'回帰直線 (r={r_val:.3f})')

ax2.set_xlabel('製造業従業者割合（%）', fontsize=12)
ax2.set_ylabel('消費支出（円/月・二人以上世帯）', fontsize=12)
ax2.set_title(f'製造業割合と消費支出の関係（47都道府県、2021/2022年）\nr={r_val:.3f}, p={p_val:.4f}',
              fontsize=13, fontweight='bold')
ax2.legend(fontsize=9, loc='upper left', ncol=2)
ax2.yaxis.get_major_formatter().set_useOffset(False)
ax2.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, _: f'{val/10000:.0f}万'))
ax2.xaxis.grid(True, alpha=0.3, linestyle='--')
ax2.yaxis.grid(True, alpha=0.3, linestyle='--')
ax2.set_axisbelow(True)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2020_H5_3_fig2.png'), bbox_inches='tight')
plt.close(fig2)
print(f"  -> 2020_H5_3_fig2.png 保存完了 (r={r_val:.3f}, p={p_val:.4f})")

print("図3 生成中...")

heat_vars = {
    '製造業\n従業者割合': '従業者割合_製造業',
    '建設業\n従業者割合': '従業者割合_建設業',
    '情報通信業\n従業者割合': '従業者割合_情報通信業',
    '医療福祉業\n従業者割合': '従業者割合_医療福祉業',
    '製造業\nLQ': 'LQ_製造業',
    '高齢化率': '高齢化率',
    '人口密度': '人口密度',
    '有効求人\n倍率(代理)': '有効求人倍率_代理',
    '転入率': '転入率',
    '消費支出': '消費支出（二人以上の世帯）',
}

heat_df = df[[v for v in heat_vars.values()]].copy()
heat_df.columns = list(heat_vars.keys())
corr_matrix = heat_df.corr()

fig3, ax3 = plt.subplots(figsize=(10, 8))

n = len(corr_matrix)
# ヒートマップ描画
im = ax3.imshow(corr_matrix.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')

# カラーバー
cbar = fig3.colorbar(im, ax=ax3, fraction=0.046, pad=0.04)
cbar.set_label('Pearson相関係数', fontsize=10)

# 軸ラベル
labels = list(corr_matrix.columns)
ax3.set_xticks(range(n))
ax3.set_yticks(range(n))
ax3.set_xticklabels(labels, fontsize=9, rotation=45, ha='right')
ax3.set_yticklabels(labels, fontsize=9)

# セル内に数値を表示
for i in range(n):
    for j in range(n):
        val = corr_matrix.values[i, j]
        text_color = 'white' if abs(val) > 0.6 else 'black'
        ax3.text(j, i, f'{val:.2f}', ha='center', va='center',
                 fontsize=8.5, color=text_color, fontweight='bold' if abs(val) > 0.5 else 'normal')

ax3.set_title('産業構造・経済指標の相関ヒートマップ（47都道府県）\n─ Pearson相関係数 ─',
              fontsize=13, fontweight='bold', pad=14)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2020_H5_3_fig3.png'), bbox_inches='tight')
plt.close(fig3)
print("  -> 2020_H5_3_fig3.png 保存完了")

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
import warnings
warnings.filterwarnings('ignore')

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'
DATA_E  = 'data/raw/SSDSE-E-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

print("=" * 60)
print("産業別従業者割合 (2021年度) 統計サマリー")
print("=" * 60)
for col in ['従業者割合_製造業', '従業者割合_建設業', '従業者割合_情報通信業']:
    m = df[col].mean()
    s = df[col].std()
    print(f"  {col}: 平均={m:.1f}%, 標準偏差={s:.1f}%")

print()
print("=" * 60)
print("LQ（製造業特化係数）上位・下位5都道府県")
print("=" * 60)
lq_sorted = df[['都道府県', 'LQ_製造業', '従業者割合_製造業']].sort_values('LQ_製造業', ascending=False)
print("製造業LQ 上位5:")
print(lq_sorted.head(5).to_string(index=False))
print("製造業LQ 下位5:")
print(lq_sorted.tail(5).to_string(index=False))

print()
print("=" * 60)
print("消費支出との相関分析")
print("=" * 60)
corr_vars = ['従業者割合_製造業', '従業者割合_建設業', '従業者割合_情報通信業',
             'LQ_製造業', '高齢化率', '人口密度', '有効求人倍率_代理', '転入率']
for var in corr_vars:
    r, p = stats.pearsonr(df[var].dropna(), df.loc[df[var].notna(), '消費支出（二人以上の世帯）'])
    sig = "**" if p < 0.01 else "*" if p < 0.05 else "  "
    print(f"  {var:20s}: r={r:+.3f} (p={p:.4f}) {sig}")

print()
print("=" * 60)
print("OLS 重回帰分析（目的変数：消費支出）")
print("=" * 60)

reg_vars = ['従業者割合_製造業', '従業者割合_情報通信業',
            '高齢化率', '人口密度', '有効求人倍率_代理']

df_reg = df[['消費支出（二人以上の世帯）'] + reg_vars].dropna()
y = df_reg['消費支出（二人以上の世帯）'].values
X = df_reg[reg_vars].values

# 標準化
y_std = (y - y.mean()) / y.std()
X_std = (X - X.mean(axis=0)) / X.std(axis=0)

X_sm = sm.add_constant(X_std)
model = sm.OLS(y_std, X_sm).fit()
print(model.summary())

# 標準化偏回帰係数
coefs = model.params[1:]  # 定数項除く
pvals = model.pvalues[1:]
conf = pd.DataFrame(model.conf_int()).iloc[1:]

print()
print("標準化偏回帰係数:")
for name, c, p in zip(reg_vars, coefs, pvals):
    sig = "**" if p < 0.01 else "*" if p < 0.05 else "  "
    print(f"  {name:25s}: β={c:+.4f} (p={p:.4f}) {sig}")

print("図4 生成中...")

var_labels_jp = {
    '従業者割合_製造業': '製造業\n従業者割合',
    '従業者割合_情報通信業': '情報通信業\n従業者割合',
    '高齢化率': '高齢化率',
    '人口密度': '人口密度',
    '有効求人倍率_代理': '有効求人倍率\n（代理）',
}

coef_vals = np.array(model.params[1:])
pval_vals = np.array(model.pvalues[1:])
ci_low    = pd.DataFrame(model.conf_int()).iloc[1:, 0].values
ci_high   = pd.DataFrame(model.conf_int()).iloc[1:, 1].values

labels_jp = [var_labels_jp[v] for v in reg_vars]
colors_bar = ['#e05c5c' if c > 0 else '#4e9af1' for c in coef_vals]
sig_markers = ['**' if p < 0.01 else '*' if p < 0.05 else '' for p in pval_vals]

fig4, ax4 = plt.subplots(figsize=(9, 6))

y_pos = np.arange(len(labels_jp))
bars = ax4.barh(y_pos, coef_vals, color=colors_bar, alpha=0.82, edgecolor='white', linewidth=0.5, height=0.55)

# 信頼区間
for i, (lo, hi) in enumerate(zip(ci_low, ci_high)):
    ax4.plot([lo, hi], [i, i], 'k-', linewidth=2.5, alpha=0.7)
    ax4.plot([lo, lo], [i-0.08, i+0.08], 'k-', linewidth=1.5, alpha=0.7)
    ax4.plot([hi, hi], [i-0.08, i+0.08], 'k-', linewidth=1.5, alpha=0.7)

# 係数値と有意性マーカー
for i, (c, sig) in enumerate(zip(coef_vals, sig_markers)):
    offset = 0.02 if c >= 0 else -0.02
    align = 'left' if c >= 0 else 'right'
    ax4.text(c + offset, i, f'{c:+.3f}{sig}',
             va='center', ha=align, fontsize=10, fontweight='bold',
             color='#333333')

ax4.axvline(x=0, color='black', linewidth=1.2, alpha=0.7)
ax4.set_yticks(y_pos)
ax4.set_yticklabels(labels_jp, fontsize=10)
ax4.set_xlabel('標準化偏回帰係数（β）', fontsize=12)
ax4.set_title('重回帰分析：地域消費支出の決定要因\n（標準化偏回帰係数と95%信頼区間、n=47都道府県）',
              fontsize=13, fontweight='bold')
ax4.xaxis.grid(True, alpha=0.3, linestyle='--')
ax4.set_axisbelow(True)

# 凡例（有意性）
ax4.text(0.99, 0.02, '**p<0.01  *p<0.05  (95%CI)', transform=ax4.transAxes,
         ha='right', va='bottom', fontsize=9, color='#555555',
         bbox=dict(boxstyle='round,pad=0.3', fc='white', alpha=0.7))

# 統計値の注記
r2_adj = model.rsquared_adj
ax4.text(0.99, 0.98, f'修正済みR²={r2_adj:.3f}',
         transform=ax4.transAxes, ha='right', va='top',
         fontsize=10, color='#222222',
         bbox=dict(boxstyle='round,pad=0.3', fc='#fff9c4', ec='#f9a825', alpha=0.9))

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2020_H5_3_fig4.png'), bbox_inches='tight')
plt.close(fig4)
print(f"  -> 2020_H5_3_fig4.png 保存完了 (修正済みR²={r2_adj:.3f})")

print()
print("=" * 60)
print("最終統計サマリー（HTML作成用）")
print("=" * 60)

# 製造業LQが高い都道府県（製造業特化地域）
lq_high = df[df['LQ_製造業'] >= 1.5][['都道府県', 'LQ_製造業', '従業者割合_製造業', '消費支出（二人以上の世帯）']].sort_values('LQ_製造業', ascending=False)
print("製造業LQ≥1.5の都道府県:")
print(lq_high.to_string(index=False))

print()
mfg_csp_r, mfg_csp_p = stats.pearsonr(
    df['従業者割合_製造業'].dropna(),
    df.loc[df['従業者割合_製造業'].notna(), '消費支出（二人以上の世帯）']
)
print(f"製造業割合 vs 消費支出: r={mfg_csp_r:.3f}, p={mfg_csp_p:.4f}")

ict_csp_r, ict_csp_p = stats.pearsonr(
    df['従業者割合_情報通信業'].dropna(),
    df.loc[df['従業者割合_情報通信業'].notna(), '消費支出（二人以上の世帯）']
)
print(f"情報通信業割合 vs 消費支出: r={ict_csp_r:.3f}, p={ict_csp_p:.4f}")

aging_csp_r, aging_csp_p = stats.pearsonr(
    df['高齢化率'].dropna(),
    df.loc[df['高齢化率'].notna(), '消費支出（二人以上の世帯）']
)
print(f"高齢化率 vs 消費支出: r={aging_csp_r:.3f}, p={aging_csp_p:.4f}")

mfg_mean = df['従業者割合_製造業'].mean()
mfg_std  = df['従業者割合_製造業'].std()
csp_mean = df['消費支出（二人以上の世帯）'].mean()
csp_std  = df['消費支出（二人以上の世帯）'].std()

print()
print(f"製造業従業者割合 全国平均: {mfg_mean:.1f}% ± {mfg_std:.1f}%")
print(f"消費支出 全国平均: {csp_mean:.0f}円 ± {csp_std:.0f}円")

print()
print("全図生成完了。")