市区町村ごとの失業率の要因分析

import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

from scipy import stats as scipy_stats
from scipy.stats import gaussian_kde
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor

# ── パス設定 ─────────────────────────────────────────────────────────────────
FIG_DIR  = 'html/figures'
DATA_DIR = 'data/raw'
os.makedirs(FIG_DIR, exist_ok=True)

plt.rcParams.update({
    'font.family':        'Hiragino Sans',
    'axes.unicode_minus': False,
    'figure.dpi':         150,
    'axes.spines.top':    False,
    'axes.spines.right':  False,
})

print("=" * 65)
print("■ データ読み込み（SSDSE-B-2026 実データのみ）")
print("=" * 65)

df_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'),
    encoding='cp932', header=1
)

# 47都道府県のみ（地域コード R + 5桁数字、かつ全国合計 R00000 を除く）
df_raw = df_raw[df_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()
df_raw['年度'] = pd.to_numeric(df_raw['年度'], errors='coerce')

# 2022年データを使用（最新の完全データ年）
YEAR = 2022
df = df_raw[df_raw['年度'] == YEAR].copy().reset_index(drop=True)
print(f"SSDSE-B {YEAR}年: {len(df)}都道府県")

# 数値変換
NUM_COLS = [
    '月間有効求職者数（一般）', '月間有効求人数（一般）',
    '総人口', '15～64歳人口', '15～64歳人口（男）', '65歳以上人口',
    '消費支出（二人以上の世帯）',
    '標準価格（平均価格）（住宅地）',
    '合計特殊出生率',
    '高等学校卒業者数', '高等学校卒業者のうち進学者数',
    '年平均気温',
    '延べ宿泊者数',
]
for c in NUM_COLS:
    df[c] = pd.to_numeric(df[c], errors='coerce')

# ── 変数の構築 ────────────────────────────────────────────────────────────────
# 目的変数：求職圧力指数（Jobseeker Pressure Index）
df['求職圧力指数'] = (
    df['月間有効求職者数（一般）'] /
    (df['月間有効求職者数（一般）'] + df['月間有効求人数（一般）'])
    * 100
)

# 説明変数
df['高齢化率']   = df['65歳以上人口'] / df['総人口'] * 100
df['男性率']     = df['15～64歳人口（男）'] / df['15～64歳人口'] * 100
df['消費支出']   = df['消費支出（二人以上の世帯）']          # 円/月
df['住宅地価格'] = df['標準価格（平均価格）（住宅地）']        # 円/m²
df['TFR']        = df['合計特殊出生率']
df['大学進学率'] = df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数'] * 100
df['年平均気温'] = df['年平均気温']
df['宿泊者per_capita'] = df['延べ宿泊者数'] / df['総人口']

# 分析変数名
Y_NAME = '求職圧力指数'
X_NAMES = ['高齢化率', '男性率', '消費支出', '住宅地価格', 'TFR', '大学進学率', '年平均気温', '宿泊者per_capita']

# 表示用変数名（短縮）
X_LABELS = {
    '高齢化率':        '高齢化率\n（%）',
    '男性率':          '男性率\n（生産年齢，%）',
    '消費支出':        '消費支出\n（円/月）',
    '住宅地価格':      '住宅地価格\n（円/m²）',
    'TFR':             '合計特殊\n出生率',
    '大学進学率':      '大学進学率\n（%）',
    '年平均気温':      '年平均気温\n（℃）',
    '宿泊者per_capita':'宿泊者\nper capita',
}

# 欠損除外
VARS_ALL = [Y_NAME] + X_NAMES + ['都道府県']
df_ana = df[VARS_ALL].dropna().reset_index(drop=True)
N = len(df_ana)
prefs = df_ana['都道府県'].tolist()

y = df_ana[Y_NAME].values
X = df_ana[X_NAMES].values

print(f"\n分析対象: {N}都道府県")
print(f"\n{Y_NAME} 記述統計:")
print(f"  平均   = {y.mean():.2f}%")
print(f"  標準偏差 = {y.std():.2f}%")
print(f"  最小   = {y.min():.2f}%  ({prefs[y.argmin()]})")
print(f"  最大   = {y.max():.2f}%  ({prefs[y.argmax()]})")

print("\n図1: 分布図（ヒストグラム + KDE）を作成中...")

fig1, ax1 = plt.subplots(figsize=(9, 5))

# ヒストグラム
n_bins = 10
counts, bin_edges, _ = ax1.hist(
    y, bins=n_bins, color='#1565C0', alpha=0.65, edgecolor='white',
    label='度数', density=True
)

# KDE（カーネル密度推定）
kde = gaussian_kde(y, bw_method='scott')
x_kde = np.linspace(y.min() - 1, y.max() + 1, 300)
ax1.plot(x_kde, kde(x_kde), color='#E65100', linewidth=2.5, label='KDE（カーネル密度）')

# 平均・中央値の縦線
ax1.axvline(y.mean(),   color='#1565C0', linestyle='--', linewidth=1.8,
            label=f'平均 = {y.mean():.2f}%')
ax1.axvline(np.median(y), color='#2E7D32', linestyle=':', linewidth=1.8,
            label=f'中央値 = {np.median(y):.2f}%')

ax1.set_xlabel('求職圧力指数（%）', fontsize=12)
ax1.set_ylabel('密度', fontsize=12)
ax1.set_title(
    f'求職圧力指数の分布 — 47都道府県（{YEAR}年）\n'
    f'（月間有効求職者数 / (求職者数 + 求人数) × 100）\nデータ：SSDSE-B-2026 実データ',
    fontsize=12, fontweight='bold'
)
ax1.legend(fontsize=10)
ax1.grid(axis='y', alpha=0.3)

# 上位・下位都道府県ラベル
highlight = ['北海道', '沖縄県', '東京都', '愛知県', '秋田県']
for pref in highlight:
    if pref in prefs:
        idx = prefs.index(pref)
        kde_val = float(kde(np.array([y[idx]]))[0])
        ax1.annotate(
            pref.replace('道','').replace('府','').replace('都','').replace('県',''),
            xy=(y[idx], 0.005),
            xytext=(y[idx], kde_val * 0.55),
            fontsize=8, color='#C62828',
            arrowprops=dict(arrowstyle='->', color='#C62828', lw=1.0),
            ha='center'
        )

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

print("図2: 相関係数ヒートマップを作成中...")

fig2, ax2 = plt.subplots(figsize=(10, 9))
corr_vals = corr_matrix.values
n_all = len(corr_matrix.columns)
col_labels_short = [X_LABELS.get(c, c) if c != Y_NAME else '求職圧力\n指数' for c in corr_matrix.columns]

im = ax2.imshow(corr_vals, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
cbar = plt.colorbar(im, ax=ax2, fraction=0.046, pad=0.04)
cbar.set_label('Pearson相関係数', fontsize=10)

ax2.set_xticks(range(n_all))
ax2.set_yticks(range(n_all))
ax2.set_xticklabels(col_labels_short, fontsize=8, rotation=35, ha='right')
ax2.set_yticklabels(col_labels_short, fontsize=8)
ax2.set_title(
    f'Pearson相関係数行列（n={N}都道府県）\n* p<0.05  ** p<0.01  *** p<0.001\nデータ：SSDSE-B-2026',
    fontsize=12, fontweight='bold'
)

for i in range(n_all):
    for j in range(n_all):
        val = corr_vals[i, j]
        sig_str = ''
        if i != j:
            try:
                _, pv = scipy_stats.pearsonr(
                    df_corr.iloc[:, i].dropna(),
                    df_corr.iloc[:, j].dropna()
                )
                sig_str = '***' if pv < 0.001 else '**' if pv < 0.01 else '*' if pv < 0.05 else ''
            except Exception:
                pass
        text_color = 'white' if abs(val) > 0.6 else 'black'
        ax2.text(j, i, f'{val:.2f}{sig_str}', ha='center', va='center',
                 fontsize=7.5, fontweight='bold', color=text_color)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2023_U5_5_fig2_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("  -> 2023_U5_5_fig2_corr.png 保存完了")

print("図3: 標準化係数棒グラフを作成中...")

# 95% CI（±1.96×SE）
ci95 = 1.96 * std_bse

# 有意性フラグで色分け
bar_colors = []
for pv in std_pvals:
    if pv < 0.01:
        bar_colors.append('#C62828')   # 深紅：高有意
    elif pv < 0.05:
        bar_colors.append('#E65100')   # オレンジ：有意
    else:
        bar_colors.append('#90CAF9')   # 薄青：非有意

short_labels = [X_LABELS.get(xn, xn) for xn in X_NAMES]

# 標準化係数の降順ソート
sort_idx = np.argsort(std_coefs)
sorted_labels  = [short_labels[i] for i in sort_idx]
sorted_coefs   = std_coefs[sort_idx]
sorted_ci      = ci95[sort_idx]
sorted_pvals   = std_pvals[sort_idx]
sorted_colors  = [bar_colors[i] for i in sort_idx]

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

bars3 = ax3.barh(
    np.arange(len(X_NAMES)), sorted_coefs,
    color=sorted_colors, alpha=0.85, edgecolor='white',
    xerr=sorted_ci, error_kw=dict(ecolor='#333', capsize=4, linewidth=1.2)
)
ax3.axvline(0, color='black', linewidth=1.0)
ax3.set_yticks(np.arange(len(X_NAMES)))
ax3.set_yticklabels(sorted_labels, fontsize=10)
ax3.set_xlabel('標準化回帰係数（β）', fontsize=12)
ax3.set_title(
    f'OLS 重回帰：標準化係数（β）と95%信頼区間\n'
    f'目的変数：求職圧力指数　R²={ols_result.rsquared:.3f}  AdjR²={ols_result.rsquared_adj:.3f}  n={N}',
    fontsize=12, fontweight='bold'
)
ax3.grid(axis='x', alpha=0.3)
ax3.invert_yaxis()

# 有意性ラベル
for i, (bar, pv) in enumerate(zip(bars3, sorted_pvals)):
    w = bar.get_width()
    sig_str = '***' if pv < 0.001 else '**' if pv < 0.01 else '*' if pv < 0.05 else 'ns'
    offset = sorted_ci[i] + 0.02
    ax3.text(w + offset if w >= 0 else w - offset - 0.05,
             bar.get_y() + bar.get_height() / 2,
             sig_str, va='center', fontsize=10,
             color='#C62828' if sig_str != 'ns' else '#999')

# 凡例
from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor='#C62828', alpha=0.85, label='p < 0.01（高度有意）'),
    Patch(facecolor='#E65100', alpha=0.85, label='p < 0.05（有意）'),
    Patch(facecolor='#90CAF9', alpha=0.85, label='p ≥ 0.05（非有意）'),
]
ax3.legend(handles=legend_elements, fontsize=9, loc='lower right')

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

print("図4: 散布図（消費支出 vs 求職圧力指数）を作成中...")

x_scatter = df_ana['消費支出'].values
y_scatter = y

r_s, p_s = scipy_stats.pearsonr(x_scatter, y_scatter)
coef_fit  = np.polyfit(x_scatter, y_scatter, 1)
x_fit     = np.linspace(x_scatter.min(), x_scatter.max(), 200)
y_fit     = np.polyval(coef_fit, x_fit)

fig4, ax4 = plt.subplots(figsize=(10, 7))

scatter = ax4.scatter(
    x_scatter / 1000, y_scatter,    # 千円単位に変換
    c='#1565C0', alpha=0.75, s=70,
    edgecolors='white', linewidth=0.6, zorder=3
)
ax4.plot(x_fit / 1000, y_fit, color='#E65100', linewidth=2.2,
         linestyle='--', label=f'回帰直線 (r={r_s:.3f})', zorder=2)

# 都道府県ラベル（代表的な都道府県）
LABEL_PREFS = ['北海道', '東京都', '大阪府', '愛知県', '沖縄県',
               '秋田県', '青森県', '鹿児島県', '福岡県', '神奈川県']
for i, pref in enumerate(prefs):
    if pref in LABEL_PREFS:
        short = pref.replace('都','').replace('道','').replace('府','').replace('県','')
        ax4.annotate(
            short, (x_scatter[i] / 1000, y_scatter[i]),
            textcoords='offset points', xytext=(5, 3),
            fontsize=8.5, color='#333',
            bbox=dict(boxstyle='round,pad=0.2', facecolor='white', alpha=0.7, edgecolor='none')
        )

ax4.set_xlabel('消費支出（二人以上の世帯）[千円/月]', fontsize=12)
ax4.set_ylabel('求職圧力指数（%）', fontsize=12)
ax4.set_title(
    f'消費支出と求職圧力指数の関係（都道府県別, {YEAR}年）\n'
    f'r = {r_s:.3f}  {"p < 0.05" if p_s < 0.05 else f"p = {p_s:.3f}"}  '
    f'回帰係数 = {coef_fit[0]*1000:.4f} (%/千円)',
    fontsize=12, fontweight='bold'
)
ax4.legend(fontsize=11)
ax4.grid(True, alpha=0.25)

# テキストボックス：解釈
ax4.text(
    0.03, 0.95,
    f'r = {r_s:.3f}\n{"p < 0.05" if p_s < 0.05 else "p = {:.3f}".format(p_s)}\n'
    f'{"消費支出が高いほど\n求職圧力が低い傾向" if r_s < 0 else "消費支出が高いほど\n求職圧力が高い傾向"}',
    transform=ax4.transAxes, fontsize=10, va='top',
    bbox=dict(boxstyle='round', facecolor='#E3F2FD', alpha=0.85)
)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2023_U5_5_fig4_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("  -> 2023_U5_5_fig4_scatter.png 保存完了")

df_corr = df_ana[X_NAMES + [Y_NAME]].copy()
corr_matrix = df_corr.corr()

print("\n【説明変数と求職圧力指数の相関】")
for xn in X_NAMES:
    r, p = scipy_stats.pearsonr(df_ana[xn], y)
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'ns'
    print(f"  {xn:<16} r = {r:+.3f}  p = {p:.4f}  {sig}")

X_reg = sm.add_constant(X)
vif_vals = [variance_inflation_factor(X_reg, i + 1) for i in range(len(X_NAMES))]
vif_df = pd.DataFrame({'変数': X_NAMES, 'VIF': vif_vals})

print("\n【VIF（分散拡大係数）】")
for i, (xn, v) in enumerate(zip(X_NAMES, vif_vals)):
    flag = '  ★多重共線性の疑い' if v > 5 else ''
    print(f"  {xn:<16} VIF = {v:.2f}{flag}")

ols_result = sm.OLS(y, X_reg).fit()
print("\n【OLS 重回帰結果（非標準化）】")
print(f"  R² = {ols_result.rsquared:.3f}  AdjR² = {ols_result.rsquared_adj:.3f}")
print(f"  F統計量 = {ols_result.fvalue:.2f}  p = {ols_result.f_pvalue:.4f}")
print(ols_result.summary().tables[1])

# 標準化 OLS（Zスコア変換）
y_z = (y - y.mean()) / y.std()
X_z = (X - X.mean(axis=0)) / X.std(axis=0)
X_z_const = sm.add_constant(X_z)
ols_std = sm.OLS(y_z, X_z_const).fit()

std_coefs = ols_std.params[1:]        # 定数項を除く
std_pvals = ols_std.pvalues[1:]
std_bse   = ols_std.bse[1:]           # 標準誤差（95%CI用）

print("\n【標準化回帰係数】")
for xn, sc, pv in zip(X_NAMES, std_coefs, std_pvals):
    sig = '***' if pv < 0.001 else '**' if pv < 0.01 else '*' if pv < 0.05 else 'ns'
    print(f"  {xn:<16} β = {sc:+.3f}  p = {pv:.4f}  {sig}")

print("\n" + "=" * 65)
print("✓ 全図の生成完了（4枚）")
print("=" * 65)
print(f"\n保存先: {os.path.abspath(FIG_DIR)}")
print("  2023_U5_5_fig1_dist.png    - 求職圧力指数の分布（ヒストグラム+KDE）")
print("  2023_U5_5_fig2_corr.png    - 相関係数ヒートマップ")
print("  2023_U5_5_fig3_coef.png    - OLS標準化係数棒グラフ（エラーバー付き）")
print("  2023_U5_5_fig4_scatter.png - 消費支出 vs 求職圧力指数散布図")
print()
print("【主要知見】")
print(f"  分析対象: {N}都道府県（SSDSE-B {YEAR}年, 実データ）")
print(f"  求職圧力指数: 平均 {y.mean():.2f}%  SD={y.std():.2f}%")
print(f"  OLS  R² = {ols_result.rsquared:.3f}  AdjR² = {ols_result.rsquared_adj:.3f}")
print(f"  VIF 最大値: {max(vif_vals):.2f}  最小値: {min(vif_vals):.2f}")
sig_vars = [X_NAMES[i] for i, p in enumerate(std_pvals) if p < 0.05]
print(f"  有意な説明変数（p<0.05）: {sig_vars if sig_vars else '—（なし）'}")
print(f"  r（消費支出 vs 求職圧力）: {r_s:.3f}  p={p_s:.4f}")
print()
print("  ※ 合成データ・乱数は一切使用していません")
print("  ※ データ出典: 統計センター SSDSE-B-2026")

import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.linear_model import ElasticNetCV
from sklearn.ensemble import RandomForestRegressor
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
import shap

# SSDSE-A（市区町村データ、原論文と同じデータ種別）
a = pd.read_csv('data/raw/SSDSE-A-2025.csv', encoding='cp932', header=2)
a = a[a['地域コード'].str.match(r'^R\d+', na=False)].copy()
num = lambda c: pd.to_numeric(a[c], errors='coerce')

# 目的変数: 失業率 = 完全失業者数 ÷ (就業者数 + 完全失業者数) × 100
labor = num('就業者数') + num('完全失業者数')
y_all = num('完全失業者数') / labor * 100

X_all = pd.DataFrame({
    '高齢化率': num('65歳以上人口') / num('総人口') * 100,
    '人口密度': num('総人口') / num('可住地面積'),
    '転入率': num('転入者数（日本人移動者）') / num('総人口') * 100,
    '単独世帯割合': num('単独世帯数') / num('一般世帯数') * 100,
    '第1次産業割合': num('第１次産業就業者数') / num('就業者数') * 100,
    '第2次産業割合': num('第２次産業就業者数') / num('就業者数') * 100,
    '第3次産業割合': num('第３次産業就業者数') / num('就業者数') * 100,
    '経常収支比率': num('経常収支比率（市町村財政）'),
    '小売店密度': num('小売店数') / num('総人口') * 10000,
    '医療福祉従業者割合': num('従業者数（民営）（医療、福祉）') / num('従業者数（民営）') * 100,
})
mask = ~(y_all.isna() | X_all.isna().any(axis=1)) & (labor > 0)
X_all, y_all = X_all[mask].reset_index(drop=True), y_all[mask].reset_index(drop=True)
names = a.loc[mask.values, ['都道府県', '市区町村']].reset_index(drop=True)
print(f'市区町村数: {len(X_all)}（欠損除外後）  失業率 平均 {y_all.mean():.2f}% / 最大 {y_all.max():.2f}%')

# ── (1) 重回帰（OLS） ──
Xz = pd.DataFrame(StandardScaler().fit_transform(X_all), columns=X_all.columns)
ols = sm.OLS(y_all, sm.add_constant(Xz)).fit()
print('\n=== (1) 重回帰（標準化係数、上位5） ===')
top = ols.params.drop('const').abs().sort_values(ascending=False).head(5)
for n in top.index:
    print(f'  {n:12s} {ols.params[n]:+.3f} (p={ols.pvalues[n]:.1e})')

# ── (2) Elastic Net（CV で正則化を選択） ──
en = ElasticNetCV(cv=5, l1_ratio=[.1, .5, .9], random_state=0).fit(Xz, y_all)
kept = [(n, c) for n, c in zip(X_all.columns, en.coef_) if abs(c) > 1e-10]
print(f'=== (2) Elastic Net === alpha={en.alpha_:.4f}, l1_ratio={en.l1_ratio_}, 残った変数 {len(kept)}/10')

# ── (3) RandomForest + SHAP ──
X_tr, X_te, y_tr, y_te = train_test_split(X_all, y_all, test_size=0.2, random_state=42)
rf = RandomForestRegressor(n_estimators=300, random_state=42).fit(X_tr, y_tr)
print(f'=== (3) RandomForest === テストR2 = {r2_score(y_te, rf.predict(X_te)):.3f}')
sv = shap.TreeExplainer(rf).shap_values(X_te)
imp = pd.Series(np.abs(sv).mean(axis=0), index=X_all.columns).sort_values(ascending=False)
print('  SHAP 上位3:', '、'.join(f'{n}({v:.3f})' for n, v in imp.head(3).items()))

# ── (4) PCA + k平均法 ──
pca = PCA(n_components=2)
pc = pca.fit_transform(Xz)
km = KMeans(n_clusters=4, random_state=0, n_init=10).fit(pc)
print(f'=== (4) PCA + k-means === PC1+PC2 寄与率 = {pca.explained_variance_ratio_.sum()*100:.1f}%')
for k in range(4):
    m = km.labels_ == k
    print(f'  クラスタ{k}: {m.sum():4d} 市区町村  平均失業率 {y_all[m].mean():.2f}%  例: '
          + '、'.join((names['市区町村'][m]).head(3)))