女性就業率と賃金格差の地域要因分析

import os
import warnings
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import statsmodels.api as sm
from scipy import stats
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import AgglomerativeClustering
from scipy.cluster.hierarchy import dendrogram, linkage

warnings.filterwarnings('ignore')

# ── パス設定 ──────────────────────────────────────────────────────────────────
FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'
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,
})

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

# ── データ読み込み ─────────────────────────────────────────────────────────────
print("データ読み込み中...")
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)

# 2022年断面データ
df = df_b[df_b['年度'] == 2022].copy()
print(f"2022年断面データ: {len(df)}都道府県")

# ── 変数計算 ─────────────────────────────────────────────────────────────────
df['女性15歳以上人口'] = df['総人口（女）'] - df['15歳未満人口（女）']
df['女性就業率']   = df['就職件数（一般）'] / df['女性15歳以上人口'] * 100
df['保育所密度']   = df['保育所等数'] / df['総人口'] * 1000
df['高齢化率']     = df['65歳以上人口'] / df['総人口'] * 100
df['求人倍率']     = df['月間有効求人数（一般）'] / df['月間有効求職者数（一般）']
df['出生率']       = df['合計特殊出生率'].astype(float)
df['消費支出_万']  = df['消費支出（二人以上の世帯）'] / 10000
df['地域区分']     = df['都道府県'].map(region_map)

# 分析変数リスト
EXOG_VARS = ['保育所密度', '高齢化率', '求人倍率', '出生率', '消費支出_万']
EXOG_LABELS = {
    '保育所密度': '保育所密度\n（人口千人あたり）',
    '高齢化率':   '高齢化率 (%)',
    '求人倍率':   '有効求人倍率',
    '出生率':     '合計特殊出生率',
    '消費支出_万': '消費支出（万円）',
}

# ── 重回帰分析（標準化係数） ───────────────────────────────────────────────────
print("\n重回帰分析（OLS）...")
X_raw = df[EXOG_VARS].copy().astype(float)
y = df['女性就業率'].astype(float)

# 標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_raw)
X_scaled_df = pd.DataFrame(X_scaled, columns=EXOG_VARS, index=df.index)
X_sm = sm.add_constant(X_scaled_df)

model = sm.OLS(y, X_sm).fit()
print(model.summary())

# 係数・CI抽出（定数項を除く）
coef   = model.params[EXOG_VARS].values
ci_low = model.conf_int().loc[EXOG_VARS, 0].values
ci_up  = model.conf_int().loc[EXOG_VARS, 1].values
pvals  = model.pvalues[EXOG_VARS].values
r2     = model.rsquared
r2_adj = model.rsquared_adj
f_stat = model.fvalue
f_pval = model.f_pvalue

# ── クラスタリング ─────────────────────────────────────────────────────────────
print("\nクラスタリング（Ward法, 4クラスター）...")
cluster_vars = ['女性就業率', '保育所密度', '高齢化率', '求人倍率', '消費支出_万']
X_cl = df[cluster_vars].astype(float)
X_cl_scaled = StandardScaler().fit_transform(X_cl)

# Ward法階層クラスタリング（デンドログラム用）
Z = linkage(X_cl_scaled, method='ward')

# 4クラスターに分割
agg = AgglomerativeClustering(n_clusters=4, linkage='ward')
df['クラスター'] = agg.fit_predict(X_cl_scaled) + 1  # 1〜4に変換

print("クラスター別都道府県数:")
print(df['クラスター'].value_counts().sort_index())
print("\nクラスター別女性就業率平均:")
print(df.groupby('クラスター')['女性就業率'].mean().round(3))

# ── 統計量出力 ────────────────────────────────────────────────────────────────
print("\n=== 記述統計 ===")
desc_vars = ['女性就業率', '保育所密度', '高齢化率', '求人倍率', '出生率', '消費支出_万']
print(df[desc_vars].describe().round(3))

print(f"\n=== 重回帰分析結果 ===")
print(f"R² = {r2:.4f}, 調整済みR² = {r2_adj:.4f}")
print(f"F統計量 = {f_stat:.2f}, p値 = {f_pval:.2e}")
for v, c, p in zip(EXOG_VARS, coef, pvals):
    sig = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else 'n.s.'))
    print(f"  {v}: β={c:.4f}, p={p:.4f} {sig}")

print("\n図1: 散布図作成中...")
fig1, ax = plt.subplots(figsize=(11, 8))

for region, color in region_colors.items():
    mask = df['地域区分'] == region
    ax.scatter(
        df.loc[mask, '保育所密度'],
        df.loc[mask, '女性就業率'],
        c=color, label=region, s=60, alpha=0.85, zorder=3,
        edgecolors='white', linewidths=0.5,
    )

# 都道府県ラベル
for _, row in df.iterrows():
    ax.annotate(
        row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', ''),
        xy=(row['保育所密度'], row['女性就業率']),
        fontsize=6.5, ha='center', va='bottom',
        xytext=(0, 4), textcoords='offset points', color='#333',
    )

# 回帰直線
x_vals = df['保育所密度'].astype(float)
y_vals = df['女性就業率'].astype(float)
slope, intercept, r_val, p_val, _ = stats.linregress(x_vals, y_vals)
x_line = np.linspace(x_vals.min(), x_vals.max(), 100)
ax.plot(x_line, slope * x_line + intercept, color='#333', linewidth=1.5,
        linestyle='--', alpha=0.7, label=f'回帰線 (r={r_val:.3f}, p<0.001)')

ax.set_xlabel('保育所密度（人口千人あたり）', fontsize=12)
ax.set_ylabel('女性就業率（%）', fontsize=12)
ax.set_title('保育所密度と女性就業率の関係（2022年 都道府県別）', fontsize=13, fontweight='bold')

legend_patches = [mpatches.Patch(color=c, label=r) for r, c in region_colors.items()]
legend_patches.append(
    plt.Line2D([0], [0], color='#333', linestyle='--', linewidth=1.5,
               label=f'回帰線 (r={r_val:.3f})')
)
ax.legend(handles=legend_patches, loc='upper right', fontsize=9, framealpha=0.9)
ax.grid(True, alpha=0.3, linestyle='-')

fig1.tight_layout()
out1 = os.path.join(FIG_DIR, '2021_U5_2_fig1.png')
fig1.savefig(out1, dpi=150, bbox_inches='tight')
plt.close(fig1)
print(f"  -> {out1}")

print("図2: 係数プロット作成中...")
fig2, ax = plt.subplots(figsize=(9, 5.5))

var_labels = [EXOG_LABELS[v] for v in EXOG_VARS]
colors_bar = ['#e05c5c' if p < 0.05 else '#aaaaaa' for p in pvals]

y_pos = np.arange(len(EXOG_VARS))
ax.barh(y_pos, coef, color=colors_bar, alpha=0.8, height=0.5, zorder=3)

# 95%CI エラーバー
xerr_low  = coef - ci_low
xerr_high = ci_up - coef
ax.errorbar(
    coef, y_pos,
    xerr=[xerr_low, xerr_high],
    fmt='none', ecolor='#333', elinewidth=1.5, capsize=5, zorder=4,
)

ax.axvline(0, color='#333', linewidth=1.0, linestyle='-')
ax.set_yticks(y_pos)
ax.set_yticklabels(var_labels, fontsize=10)
ax.set_xlabel('標準化偏回帰係数（β）', fontsize=11)
ax.set_title(f'重回帰分析：女性就業率の決定要因（2022年, N=47)\nR²={r2:.3f}, 調整済みR²={r2_adj:.3f}',
             fontsize=12, fontweight='bold')

# p値ラベル
for i, (c, p) in enumerate(zip(coef, pvals)):
    sig = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else ''))
    x_text = ci_up[i] + 0.01 if c >= 0 else ci_low[i] - 0.01
    ha_text = 'left' if c >= 0 else 'right'
    label_text = f'p={p:.3f}{sig}' if sig else f'p={p:.2f} n.s.'
    ax.text(x_text, i, label_text, va='center', ha=ha_text, fontsize=8.5, color='#333')

sig_patch = mpatches.Patch(color='#e05c5c', alpha=0.8, label='有意（p<0.05）')
ns_patch  = mpatches.Patch(color='#aaaaaa', alpha=0.8, label='非有意')
ax.legend(handles=[sig_patch, ns_patch], loc='lower right', fontsize=9)
ax.grid(True, axis='x', alpha=0.3)

fig2.tight_layout()
out2 = os.path.join(FIG_DIR, '2021_U5_2_fig2.png')
fig2.savefig(out2, dpi=150, bbox_inches='tight')
plt.close(fig2)
print(f"  -> {out2}")

print("図3: デンドログラム作成中...")
fig3, ax = plt.subplots(figsize=(16, 6))

pref_names = df['都道府県'].str.replace('県', '').str.replace('府', '').str.replace('都', '').str.replace('道', '北海').values

dendrogram(
    Z,
    labels=pref_names,
    ax=ax,
    leaf_rotation=90,
    leaf_font_size=8,
    color_threshold=Z[-3, 2],  # 4クラスターに対応する閾値
    above_threshold_color='#999999',
)

ax.set_title('都道府県クラスタリング（Ward法、4クラスター）\n変数: 女性就業率・保育所密度・高齢化率・求人倍率・消費支出',
             fontsize=12, fontweight='bold')
ax.set_ylabel('クラスター間距離（Ward法）', fontsize=11)
ax.set_xlabel('都道府県', fontsize=11)
ax.tick_params(axis='x', labelsize=7.5)

# 4クラスター境界線
threshold = Z[-3, 2]
ax.axhline(threshold, color='red', linewidth=1.2, linestyle='--', alpha=0.7,
           label=f'4クラスター切断（距離={threshold:.2f}）')
ax.legend(fontsize=9)

fig3.tight_layout()
out3 = os.path.join(FIG_DIR, '2021_U5_2_fig3.png')
fig3.savefig(out3, dpi=150, bbox_inches='tight')
plt.close(fig3)
print(f"  -> {out3}")

print("図4: 箱ひげ図作成中...")
box_vars = ['女性就業率', '保育所密度', '高齢化率', '求人倍率', '消費支出_万']
box_titles = ['女性就業率（%）', '保育所密度\n（人口千人あたり）', '高齢化率（%）',
              '有効求人倍率', '消費支出（万円）']

cluster_palette = {1: '#4e9af1', 2: '#e05c5c', 3: '#f0a500', 4: '#5cb85c'}

fig4, axes = plt.subplots(1, 5, figsize=(16, 5.5))

for ax, var, title in zip(axes, box_vars, box_titles):
    data_by_cl = [
        df.loc[df['クラスター'] == cl, var].astype(float).values
        for cl in range(1, 5)
    ]
    bp = ax.boxplot(
        data_by_cl,
        patch_artist=True,
        medianprops=dict(color='white', linewidth=2),
        whiskerprops=dict(linewidth=1.2),
        capprops=dict(linewidth=1.2),
        flierprops=dict(marker='o', markersize=4, linestyle='none',
                        markeredgecolor='#555', markerfacecolor='none'),
    )
    for patch, cl in zip(bp['boxes'], range(1, 5)):
        patch.set_facecolor(cluster_palette[cl])
        patch.set_alpha(0.8)

    ax.set_xticklabels([f'C{i}' for i in range(1, 5)], fontsize=9)
    ax.set_title(title, fontsize=9.5, fontweight='bold')
    ax.grid(True, axis='y', alpha=0.3)

fig4.suptitle('クラスター別主要指標の分布（Ward法 4クラスター, 2022年）',
              fontsize=12, fontweight='bold', y=1.01)

legend_patches4 = [
    mpatches.Patch(color=cluster_palette[cl], alpha=0.8, label=f'クラスター{cl}')
    for cl in range(1, 5)
]
fig4.legend(handles=legend_patches4, loc='lower center', ncol=4, fontsize=9,
            bbox_to_anchor=(0.5, -0.08))

fig4.tight_layout()
out4 = os.path.join(FIG_DIR, '2021_U5_2_fig4.png')
fig4.savefig(out4, dpi=150, bbox_inches='tight')
plt.close(fig4)
print(f"  -> {out4}")

# ── 最終確認 ─────────────────────────────────────────────────────────────────
print("\n=== 生成ファイル確認 ===")
for fp in [out1, out2, out3, out4]:
    size = os.path.getsize(fp) / 1024
    print(f"  {os.path.basename(fp)}: {size:.1f} KB")

print("\n=== クラスター別都道府県一覧 ===")
for cl in range(1, 5):
    prefs = df.loc[df['クラスター'] == cl, '都道府県'].tolist()
    avg_emp = df.loc[df['クラスター'] == cl, '女性就業率'].mean()
    avg_age = df.loc[df['クラスター'] == cl, '高齢化率'].mean()
    avg_nur = df.loc[df['クラスター'] == cl, '保育所密度'].mean()
    print(f"  クラスター{cl} (n={len(prefs)}): {', '.join(prefs)}")
    print(f"    女性就業率={avg_emp:.3f}%, 高齢化率={avg_age:.1f}%, 保育所密度={avg_nur:.3f}")

print("\n完了!")