都市農村間の所得格差
消費支出と雇用環境の地域差分析

df_raw = pd.read_csv(DATA_B, encoding='shift-jis', header=0)
# 行0がラベル行、行1以降が実データ
df_all = df_raw.iloc[1:].copy()
df_all.columns = df_raw.columns

# 数値列に変換
NUM_COLS = ['A1101', 'A1303', 'A5101', 'A5102',
            'C5401', 'E6302', 'F3102', 'F3103', 'L3221']
for c in NUM_COLS:
    df_all[c] = pd.to_numeric(df_all[c], errors='coerce')

df_all['年度']    = df_all['SSDSE-B-2026'].astype(int)
df_all['都道府県'] = df_all['Prefecture'].astype(str)

print("=" * 65)
print("■ SSDSE-B-2026 読み込み完了")
print(f"  年度範囲: {df_all['年度'].min()}〜{df_all['年度'].max()}")
print(f"  都道府県数（各年）: {df_all.groupby('年度')['都道府県'].count().iloc[0]}件")
print("=" * 65)

df_all['消費支出_万円']   = df_all['L3221'] / 10000
df_all['求人倍率']        = df_all['F3103'] / df_all['F3102']
df_all['住宅地価_log']    = np.log(df_all['C5401'])
df_all['高齢化率']        = df_all['A1303'] / df_all['A1101'] * 100
df_all['大学学生数_千人'] = df_all['E6302'] / 1000
df_all['純移動率']        = (df_all['A5101'] - df_all['A5102']) / df_all['A1101'] * 1000

# 2022年データ抽出
df_2022 = df_all[df_all['年度'] == 2022].copy().reset_index(drop=True)

# 都市部 / 地方 分類（住宅地価の全国中央値を閾値）
land_median = df_2022['C5401'].median()
df_2022['地域区分'] = df_2022['C5401'].apply(
    lambda x: '都市部' if x >= land_median else '地方'
)
df_all['地域区分'] = df_all.apply(
    lambda row: '都市部' if row['C5401'] >= land_median else '地方',
    axis=1
)

print(f"\n■ 2022年データ: {len(df_2022)}都道府県")
print(f"  住宅地価 中央値（分類閾値）: {land_median:,.0f} 円/m²")
print(f"  都市部: {(df_2022['地域区分']=='都市部').sum()}都道府県")
print(f"  地方 : {(df_2022['地域区分']=='地方').sum()}都道府県")
print(f"  消費支出（平均）: {df_2022['消費支出_万円'].mean():.2f}万円/月")
print(f"  求人倍率 （平均）: {df_2022['求人倍率'].mean():.3f}")
print(f"  高齢化率 （平均）: {df_2022['高齢化率'].mean():.1f}%")

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

print("\n" + "=" * 65)
print("■ 図1: 消費支出の時系列推移（地域別平均）")
print("=" * 65)

df_region_ts = (
    df_all.groupby(['年度', '地方'])['消費支出_万円']
    .mean()
    .reset_index()
)

fig1, ax1 = plt.subplots(figsize=(11, 6))

for region in REGION_ORDER:
    sub = df_region_ts[df_region_ts['地方'] == region].sort_values('年度')
    if len(sub) == 0:
        continue
    ax1.plot(sub['年度'], sub['消費支出_万円'],
             marker='o', linewidth=2.2, markersize=5.5,
             color=REGION_COLORS[region], label=region)

# 2020年にCOVID-19 影響を示す縦線
ax1.axvline(2020, color='#FF8F00', linestyle=':', linewidth=1.5, alpha=0.8)
ax1.text(2020.08, ax1.get_ylim()[0] + 0.5, 'COVID-19\n(2020)',
         ha='left', va='bottom', fontsize=8.5, color='#FF8F00')

# 2022年に分析基準年の縦線
ax1.axvline(2022, color='#C62828', linestyle='--', linewidth=1.5, alpha=0.8)
ax1.text(2022.08, ax1.get_ylim()[1] - 0.8, '2022年\n（分析基準年）',
         ha='left', va='top', fontsize=8.5, color='#C62828', fontweight='bold')

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('消費支出（万円/月）', fontsize=12)
ax1.set_title('消費支出の時系列推移（地方別平均, 2012〜2023年）\n〜関東が最高、九州・東北が低水準で推移〜',
              fontsize=13, fontweight='bold')
ax1.set_xticks(sorted(df_all['年度'].unique()))
ax1.xaxis.set_tick_params(rotation=45)
ax1.legend(loc='lower right', fontsize=9.5, framealpha=0.88)
ax1.grid(True, alpha=0.3)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2022_H5_14_fig1_timeseries.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig1)
print("図1保存: 2022_H5_14_fig1_timeseries.png")

print("\n■ 図2: 住宅地価 vs 消費支出 散布図（2022年）")

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

color_map = {'都市部': '#C62828', '地方': '#1565C0'}
marker_map = {'都市部': 'o', '地方': '^'}

for grp in ['都市部', '地方']:
    sub2 = df_2022[df_2022['地域区分'] == grp]
    ax2.scatter(sub2['住宅地価_log'], sub2['消費支出_万円'],
                c=color_map[grp], marker=marker_map[grp],
                s=80, alpha=0.80, edgecolors='#333', linewidth=0.4,
                label=grp, zorder=3)

# 回帰直線（全データ）
valid2 = df_2022[['住宅地価_log', '消費支出_万円']].dropna()
x2v = valid2['住宅地価_log'].values
y2v = valid2['消費支出_万円'].values
z2  = np.polyfit(x2v, y2v, 1)
xs2 = np.linspace(x2v.min(), x2v.max(), 100)
ax2.plot(xs2, np.poly1d(z2)(xs2), 'k--', linewidth=1.8, alpha=0.65,
         label=f'回帰直線 (r={r_land:.3f}, p={p_land:.4f})')

# 中央値の閾値ライン
log_med = np.log(land_median)
ax2.axvline(log_med, color='#FF8F00', linestyle=':', linewidth=1.5, alpha=0.8)
ax2.text(log_med + 0.02, ax2.get_ylim()[0] + 0.3,
         f'中央値\n({land_median:,.0f}円/m²)',
         ha='left', va='bottom', fontsize=8, color='#FF8F00')

# 都道府県ラベル（注目都市）
LABEL_PREFS = {
    '東京都', '北海道', '沖縄県', '秋田県', '青森県',
    '神奈川県', '愛知県', '大阪府', '京都府', '福岡県',
    '高知県', '島根県', '鳥取県', '宮崎県',
}
for _, row2 in df_2022.iterrows():
    if row2['都道府県'] in LABEL_PREFS:
        short = row2['都道府県'].replace('県', '').replace('都', '').replace('道', '').replace('府', '')
        ax2.annotate(short, (row2['住宅地価_log'], row2['消費支出_万円']),
                     xytext=(4, 3), textcoords='offset points',
                     fontsize=7.5, color='#1A237E', fontweight='bold')

ax2.set_xlabel('住宅地価（log 変換, 円/m²）', fontsize=12)
ax2.set_ylabel('消費支出（万円/月）', fontsize=12)
ax2.set_title('住宅地価（都市度の代理）と消費支出の関係（2022年, 都道府県別）\n'
              '〜都市部（赤○）は高地価・高消費支出, 地方（青△）は低水準〜',
              fontsize=13, fontweight='bold')
ax2.legend(fontsize=9.5)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2022_H5_14_fig2_scatter.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2022_H5_14_fig2_scatter.png")

print("\n■ 図3: OLS回帰係数プロット")

coefs3 = np.asarray(ols_model.params[1:])
ses3   = np.asarray(ols_model.bse[1:])
pvals3 = np.asarray(ols_model.pvalues[1:])

COEF_COLORS = []
for p in pvals3:
    if p < 0.01:
        COEF_COLORS.append('#C62828')
    elif p < 0.05:
        COEF_COLORS.append('#FF8F00')
    elif p < 0.10:
        COEF_COLORS.append('#FDD835')
    else:
        COEF_COLORS.append('#9E9E9E')

fig3, ax3 = plt.subplots(figsize=(9, 5.5))
y_pos3 = np.arange(len(REG_LABELS))

ax3.barh(y_pos3, coefs3, color=COEF_COLORS, alpha=0.82,
         edgecolor='white', height=0.55)
ax3.errorbar(coefs3, y_pos3, xerr=1.96 * ses3,
             fmt='none', color='#222', capsize=5, linewidth=1.5)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(y_pos3)
ax3.set_yticklabels(REG_LABELS, fontsize=11)
ax3.set_xlabel('標準化回帰係数（±1.96SE）', fontsize=11)
ax3.set_title(
    f'消費支出の決定要因 — OLS重回帰係数（2022年, N={len(df_reg)}都道府県）\n'
    f'R²={ols_model.rsquared:.3f}（adj. R²={ols_model.rsquared_adj:.3f}）',
    fontsize=12, fontweight='bold'
)
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)

# p値ラベル
for i, (c, p) in enumerate(zip(coefs3, pvals3)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    offset = 0.01
    ha = 'left' if c >= 0 else 'right'
    ax3.text(c + (offset if c >= 0 else -offset), i,
             f' {c:.3f} {sig}',
             va='center', ha=ha, fontsize=8.5)

from matplotlib.patches import Patch
ax3.legend(handles=[
    Patch(color='#C62828', alpha=0.85, label='p<0.01'),
    Patch(color='#FF8F00', alpha=0.85, label='p<0.05'),
    Patch(color='#FDD835', alpha=0.85, label='p<0.10'),
    Patch(color='#9E9E9E', alpha=0.85, label='n.s.'),
], fontsize=9, loc='lower right')
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2022_H5_14_fig3_coef.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2022_H5_14_fig3_coef.png")

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

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

print("\n" + "=" * 65)
print("■ OLS重回帰分析（2022年）")
print("=" * 65)

REG_VARS   = ['求人倍率', '住宅地価_log', '高齢化率', '大学学生数_千人', '純移動率']
REG_LABELS = ['求人倍率', '住宅地価（log）', '高齢化率（%）', '大学学生数（千人）', '純移動率（‰）']

df_reg = df_2022[['消費支出_万円'] + REG_VARS].dropna()
y_ols  = df_reg['消費支出_万円'].values
X_ols  = df_reg[REG_VARS].values

# 標準化（係数の大きさを比較するため）
X_std = (X_ols - X_ols.mean(axis=0)) / X_ols.std(axis=0)
X_with_const = sm.add_constant(X_std)
ols_model = sm.OLS(y_ols, X_with_const).fit()

print(ols_model.summary2())
print(f"\n  R²        = {ols_model.rsquared:.4f}")
print(f"  自由度修正済R² = {ols_model.rsquared_adj:.4f}")

# 個別相関分析
r_land, p_land = stats.pearsonr(df_2022['住宅地価_log'].dropna(),
                                 df_2022['消費支出_万円'].dropna())
r_job,  p_job  = stats.pearsonr(df_2022['求人倍率'].dropna(),
                                 df_2022['消費支出_万円'].dropna())
print(f"\n■ 相関分析（住宅地価_log vs 消費支出）: r={r_land:.4f}, p={p_land:.4f}")
print(f"■ 相関分析（求人倍率 vs 消費支出）    : r={r_job:.4f},  p={p_job:.4f}")

# Mann-Whitney U検定（都市部 vs 地方 消費支出）
urban = df_2022[df_2022['地域区分']=='都市部']['消費支出_万円'].dropna()
rural = df_2022[df_2022['地域区分']=='地方' ]['消費支出_万円'].dropna()
u_stat, u_pval = stats.mannwhitneyu(urban, rural, alternative='two-sided')
print(f"\n■ Mann-Whitney U検定（都市部 vs 地方 消費支出）:")
print(f"  U={u_stat:.1f}, p={u_pval:.4f}")
print(f"  都市部: n={len(urban)}, 中央値={urban.median():.2f}万円")
print(f"  地方 : n={len(rural)}, 中央値={rural.median():.2f}万円")

# ジニ係数の計算
def gini(arr):
    """ジニ係数を計算する（実データ用）"""
    a = np.sort(arr)
    n = len(a)
    return (2 * np.sum(np.arange(1, n + 1) * a) - (n + 1) * np.sum(a)) / (n * np.sum(a))

g_all   = gini(df_2022['消費支出_万円'].dropna().values)
g_urban = gini(urban.values)
g_rural = gini(rural.values)
print(f"\n■ ジニ係数（消費支出, 2022年）:")
print(f"  全国: {g_all:.4f}")
print(f"  都市部のみ: {g_urban:.4f}")
print(f"  地方のみ  : {g_rural:.4f}")

print("\n■ 図4: 都市部 vs 地方 消費支出比較 箱ひげ図")

fig4, axes4 = plt.subplots(1, 2, figsize=(12, 6))

# 左パネル: 箱ひげ図
ax4L = axes4[0]
data_urban = urban.values
data_rural = rural.values
bp = ax4L.boxplot([data_urban, data_rural],
                  labels=['都市部', '地方'],
                  patch_artist=True,
                  widths=0.5,
                  medianprops=dict(color='white', linewidth=2.5),
                  boxprops=dict(linewidth=1.5),
                  whiskerprops=dict(linewidth=1.5),
                  capprops=dict(linewidth=1.5),
                  flierprops=dict(marker='o', markersize=5, alpha=0.6))

bp['boxes'][0].set_facecolor('#C62828')
bp['boxes'][0].set_alpha(0.78)
bp['boxes'][1].set_facecolor('#1565C0')
bp['boxes'][1].set_alpha(0.78)

# 個別データ点をジッタで表示
rng = np.random.default_rng(seed=42)   # 表示用のジッタのみ（データ生成ではない）
for i, (data, color) in enumerate(zip([data_urban, data_rural], ['#C62828', '#1565C0'])):
    jitter = rng.uniform(-0.12, 0.12, size=len(data))
    ax4L.scatter(np.ones(len(data)) * (i + 1) + jitter, data,
                 color=color, alpha=0.5, s=30, zorder=3)

# 平均値をダイヤモンドでプロット
for i, data in enumerate([data_urban, data_rural]):
    ax4L.scatter(i + 1, data.mean(), marker='D', color='white', s=60, zorder=5,
                 edgecolors='#333', linewidths=1.5)

# Mann-Whitney U検定の結果を注記
sig_txt = '***' if u_pval < 0.001 else '**' if u_pval < 0.01 else '*' if u_pval < 0.05 else 'n.s.'
ax4L.text(1.5, max(data_urban.max(), data_rural.max()) * 1.01,
          f'Mann-Whitney U検定\nU={u_stat:.0f}, p={u_pval:.4f} {sig_txt}',
          ha='center', va='bottom', fontsize=9, color='#333',
          bbox=dict(boxstyle='round,pad=0.3', facecolor='#FFF9C4', edgecolor='#F9A825'))

ax4L.set_ylabel('消費支出（万円/月）', fontsize=11)
ax4L.set_title('都市部 vs 地方\n消費支出の分布（2022年）', fontsize=12, fontweight='bold')
ax4L.grid(axis='y', alpha=0.3)

# 右パネル: 年度別 格差推移（都市部平均 - 地方平均）
ax4R = axes4[1]
gap_ts = (
    df_all.groupby(['年度', '地域区分'])['消費支出_万円']
    .mean()
    .unstack('地域区分')
)
if '都市部' in gap_ts.columns and '地方' in gap_ts.columns:
    gap_ts['格差'] = gap_ts['都市部'] - gap_ts['地方']
    years = gap_ts.index.tolist()

    ax4R.fill_between(years, gap_ts['格差'], alpha=0.3, color='#C62828')
    ax4R.plot(years, gap_ts['格差'], marker='o', linewidth=2.2, markersize=5.5,
              color='#C62828', label='都市部 - 地方（格差）')
    ax4R.axhline(0, color='gray', linestyle='--', linewidth=1.0)

    ax4R.set_xlabel('年度', fontsize=11)
    ax4R.set_ylabel('消費支出の格差（万円/月）', fontsize=11)
    ax4R.set_title('都市農村間 消費支出格差の推移\n（都市部平均 − 地方平均）',
                   fontsize=12, fontweight='bold')
    ax4R.set_xticks(years)
    ax4R.xaxis.set_tick_params(rotation=45)
    ax4R.legend(fontsize=9.5)
    ax4R.grid(True, alpha=0.3)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2022_H5_14_fig4_boxplot.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig4)
print("図4保存: 2022_H5_14_fig4_boxplot.png")

print("\n" + "=" * 65)
print("■ 全図生成完了（4枚）")
print(f"  fig1_timeseries.png : 消費支出の時系列推移（地方別平均）")
print(f"  fig2_scatter.png    : 住宅地価 vs 消費支出 散布図（2022年）")
print(f"  fig3_coef.png       : OLS回帰係数プロット")
print(f"  fig4_boxplot.png    : 都市部 vs 地方 消費支出比較")
print("=" * 65)
print(f"\n  OLS結果サマリ（2022年, N={len(df_reg)}）:")
print(f"    R²={ols_model.rsquared:.3f}, adj.R²={ols_model.rsquared_adj:.3f}")
for label, coef, pval in zip(REG_LABELS, coefs3, pvals3):
    sig = '***' if pval < 0.001 else '**' if pval < 0.01 else '*' if pval < 0.05 else 'n.s.'
    print(f"    {label}: β={coef:.3f} (p={pval:.4f}) {sig}")
print(f"\n  ジニ係数（消費支出, 2022年）: 全国={g_all:.4f}, 都市部={g_urban:.4f}, 地方={g_rural:.4f}")