住宅費・居住コストの地域格差と決定要因分析

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

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'
os.makedirs(FIG_DIR, exist_ok=True)

# ── データ読み込み ──────────────────────────────────
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)

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

# ── 2022年断面データ ──────────────────────────────
df22 = df_b[df_b['年度'] == 2022].copy()
df22['地域'] = df22['都道府県'].map(region_map)

# 派生変数
df22['転入率'] = df22['転入者数（日本人移動者）'] / df22['総人口'] * 1000   # 人口千人あたり転入者数
df22['高齢化率'] = df22['65歳以上人口'] / df22['総人口'] * 100              # %
df22['住宅地価格'] = df22['標準価格（平均価格）（住宅地）']                   # 円/m²（都市集積の代理）
df22['住宅地価格_log'] = np.log(df22['住宅地価格'])
df22 = df22.sort_values('住居費（二人以上の世帯）', ascending=False).reset_index(drop=True)

Y_COL  = '住居費（二人以上の世帯）'       # 目的変数
X_COLS = ['住宅地価格_log', '転入率', '消費支出（二人以上の世帯）', '高齢化率', '合計特殊出生率']
X_LABELS = ['住宅地価格\n(log, 都市集積)', '転入率\n(需要代理)', '消費支出\n(所得水準)', '高齢化率\n(%)', '合計特殊\n出生率']

# ── OLS重回帰 ──────────────────────────────────────
Y = df22[Y_COL]
X_raw = df22[X_COLS].copy()
X_mat = sm.add_constant(X_raw)
model = sm.OLS(Y, X_mat).fit()

# 標準化偏回帰係数
Y_std = (Y - Y.mean()) / Y.std()
X_std_df = (X_raw - X_raw.mean()) / X_raw.std()
X_std_mat = sm.add_constant(X_std_df)
model_std = sm.OLS(Y_std, X_std_mat).fit()
std_coefs  = model_std.params[1:].values
std_pvals  = model_std.pvalues[1:].values

print("=" * 60)
print("OLS 重回帰分析（住居費の決定要因）")
print("=" * 60)
print(model.summary())
print()
print("標準化偏回帰係数:")
for lbl, c, p in zip(X_LABELS, std_coefs, std_pvals):
    sig = "**" if p < 0.01 else ("*" if p < 0.05 else "n.s.")
    print(f"  {lbl.replace(chr(10),'')}: β={c:.4f}  p={p:.4f}  {sig}")
print()
print(f"R² = {model.rsquared:.4f},  調整済みR² = {model.rsquared_adj:.4f}")
print(f"F統計量 = {model.fvalue:.3f},  p(F) = {model.f_pvalue:.4f}")

# ── 記述統計（出力用） ──────────────────────────────
max_val = df22[Y_COL].max()
min_val = df22[Y_COL].min()
mean_val = df22[Y_COL].mean()
max_pref = df22.loc[df22[Y_COL].idxmax(), '都道府県']
min_pref = df22.loc[df22[Y_COL].idxmin(), '都道府県']
print()
print(f"住居費 最大: {max_pref} {max_val:,.0f}円  最小: {min_pref} {min_val:,.0f}円  全国平均: {mean_val:,.0f}円")
print(f"最大/最小 比率: {max_val/min_val:.1f}倍")

# 転入率と住居費の相関
r_nyusu, p_nyusu = stats.pearsonr(df22['転入率'], df22[Y_COL])
r_price, p_price = stats.pearsonr(df22['住宅地価格'], df22[Y_COL])
r_aging, p_aging = stats.pearsonr(df22['高齢化率'], df22[Y_COL])
print(f"\n相関係数:")
print(f"  転入率 vs 住居費: r={r_nyusu:.3f}, p={p_nyusu:.4f}")
print(f"  住宅地価格 vs 住居費: r={r_price:.3f}, p={p_price:.4f}")
print(f"  高齢化率 vs 住居費: r={r_aging:.3f}, p={p_aging:.4f}")

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

df_sorted = df22.sort_values(Y_COL, ascending=True)
bar_colors = [region_colors[r] for r in df_sorted['地域']]

bars = ax1.barh(df_sorted['都道府県'], df_sorted[Y_COL] / 1000, color=bar_colors,
                edgecolor='none', height=0.75)

# 全国平均線
mean_k = mean_val / 1000
ax1.axvline(mean_k, color='#333', linewidth=1.5, linestyle='--', label=f'全国平均 {mean_val:,.0f}円')
ax1.set_xlabel('住居費（二人以上の世帯）[千円/月]', fontsize=12)
ax1.set_title('図1: 都道府県別 住居費ランキング（2022年度）', fontsize=14, fontweight='bold', pad=12)
ax1.tick_params(axis='y', labelsize=8.5)

# 凡例（地域）
from matplotlib.patches import Patch
legend_handles = [Patch(facecolor=c, label=r) for r, c in region_colors.items()]
ax1.legend(handles=legend_handles, loc='lower right', fontsize=9, framealpha=0.9)
ax1.text(mean_k + 0.3, 2, '全国平均', fontsize=9, color='#333')

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2020_U5_3_fig1.png'), dpi=150, bbox_inches='tight')
plt.close(fig1)
print("\n図1 保存完了: 2020_U5_3_fig1.png")

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

for region, grp in df22.groupby('地域'):
    ax2.scatter(grp['住宅地価格'], grp[Y_COL] / 1000,
                color=region_colors[region], label=region,
                s=55, alpha=0.85, zorder=3)

# 都道府県ラベル
for _, row in df22.iterrows():
    ax2.annotate(row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', ''),
                 (row['住宅地価格'], row[Y_COL] / 1000),
                 fontsize=6.5, xytext=(3, 2), textcoords='offset points', color='#444')

# 回帰線（対数スケール → 線形）
x_range = np.linspace(df22['住宅地価格'].min(), df22['住宅地価格'].max(), 200)
# Fit simple regression log(price) → cost
from numpy.polynomial import polynomial as P
log_x = np.log(df22['住宅地価格'])
c_fit = np.polyfit(log_x, df22[Y_COL] / 1000, 1)
ax2.plot(x_range, np.polyval(c_fit, np.log(x_range)),
         color='#333', linewidth=1.8, linestyle='--', label=f'回帰線 r={r_price:.3f}')

ax2.set_xscale('log')
ax2.set_xlabel('標準価格（住宅地）[円/m²]（対数スケール）', fontsize=11)
ax2.set_ylabel('住居費（二人以上の世帯）[千円/月]', fontsize=11)
ax2.set_title('図2: 住宅地価格（都市集積）と住居費の関係（2022年度）', fontsize=13, fontweight='bold', pad=10)
ax2.legend(fontsize=9, framealpha=0.9)
ax2.grid(True, which='major', alpha=0.3)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2020_U5_3_fig2.png'), dpi=150, bbox_inches='tight')
plt.close(fig2)
print("図2 保存完了: 2020_U5_3_fig2.png")

ts_df = df_b[df_b['年度'].between(2012, 2023)].copy()
ts_df['地域'] = ts_df['都道府県'].map(region_map)

# 地域別・年度別 平均住居費
ts_region = ts_df.groupby(['年度', '地域'])[Y_COL].mean().reset_index()

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

region_order = ['関東', '近畿', '中部', '九州・沖縄', '北海道・東北', '中国・四国']
for region in region_order:
    sub = ts_region[ts_region['地域'] == region].sort_values('年度')
    ax3.plot(sub['年度'], sub[Y_COL] / 1000,
             color=region_colors[region], linewidth=2.2,
             marker='o', markersize=5, label=region)

ax3.set_xticks(range(2012, 2024))
ax3.set_xticklabels([str(y) for y in range(2012, 2024)], rotation=45, fontsize=9)
ax3.set_xlabel('年度', fontsize=11)
ax3.set_ylabel('住居費（二人以上の世帯）[千円/月]', fontsize=11)
ax3.set_title('図3: 地域別 住居費の推移（2012〜2023年度）', fontsize=13, fontweight='bold', pad=10)
ax3.legend(loc='upper left', fontsize=10, framealpha=0.9)
ax3.grid(True, alpha=0.3)
ax3.set_xlim(2011.5, 2023.5)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2020_U5_3_fig3.png'), dpi=150, bbox_inches='tight')
plt.close(fig3)
print("図3 保存完了: 2020_U5_3_fig3.png")

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

colors_coef = []
for c, p in zip(std_coefs, std_pvals):
    if p < 0.05:
        colors_coef.append('#e05c5c' if c < 0 else '#2e86ab')
    else:
        colors_coef.append('#aaa')

bars4 = ax4.barh(X_LABELS, std_coefs, color=colors_coef,
                 edgecolor='none', height=0.55)
ax4.axvline(0, color='#333', linewidth=1.2)

# p値の注釈
for i, (c, p) in enumerate(zip(std_coefs, std_pvals)):
    sig_mark = '**' if p < 0.01 else ('*' if p < 0.05 else '')
    offset = 0.01 if c >= 0 else -0.01
    ha = 'left' if c >= 0 else 'right'
    ax4.text(c + offset, i, f' β={c:.3f}{sig_mark}', va='center', ha=ha, fontsize=10, color='#222')

ax4.set_xlabel('標準化偏回帰係数（β）', fontsize=11)
ax4.set_title('図4: 住居費に対する標準化偏回帰係数（OLS，2022年度）', fontsize=12, fontweight='bold', pad=10)
ax4.set_xlim(-0.5, 0.65)
ax4.grid(True, axis='x', alpha=0.3)

# 凡例
from matplotlib.patches import Patch
legend_h = [
    Patch(facecolor='#2e86ab', label='正の効果（p<0.05）'),
    Patch(facecolor='#e05c5c', label='負の効果（p<0.05）'),
    Patch(facecolor='#aaa',    label='非有意（p>=0.05）'),
]
ax4.legend(handles=legend_h, loc='lower right', fontsize=9, framealpha=0.9)

# 統計量テキスト
ax4.text(0.02, 0.97, f'R2={model.rsquared:.3f}  調整済みR2={model.rsquared_adj:.3f}\nN=47都道府県  *p<0.05  **p<0.01',
         transform=ax4.transAxes, fontsize=9, va='top', color='#555',
         bbox=dict(boxstyle='round,pad=0.4', fc='white', ec='#ccc', alpha=0.9))

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2020_U5_3_fig4.png'), dpi=150, bbox_inches='tight')
plt.close(fig4)
print("図4 保存完了: 2020_U5_3_fig4.png")
print("\n全図の生成が完了しました。")