パネルデータを用いた空き家に関する要因分析

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_b = df_b.sort_values(['都道府県', '年度']).reset_index(drop=True)

print(f"パネルデータ形状: {df_b.shape}  ({df_b['都道府県'].nunique()} 都道府県 × {df_b['年度'].nunique()} 年度)")
print(f"年度範囲: {df_b['年度'].min()} 〜 {df_b['年度'].max()}")

df_b['転出率']      = df_b['転出者数（日本人移動者）'] / df_b['総人口'] * 1000   # 人口千対
df_b['転入率']      = df_b['転入者数（日本人移動者）'] / df_b['総人口'] * 1000
df_b['着工率']      = df_b['着工新設住宅戸数'] / df_b['総人口'] * 10000          # 人口1万対
df_b['高齢化率']    = df_b['65歳以上人口'] / df_b['総人口'] * 100                # %
df_b['地価_log']    = np.log(df_b['標準価格（平均価格）（住宅地）'].clip(lower=1))
df_b['求人倍率']    = df_b['月間有効求人数（一般）'] / df_b['月間有効求職者数（一般）'].replace(0, np.nan)
df_b['消費支出_log']= np.log(df_b['消費支出（二人以上の世帯）'].clip(lower=1))
df_b['純移動率']    = (df_b['転入者数（日本人移動者）'] - df_b['転出者数（日本人移動者）']) / df_b['総人口'] * 1000

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

print("\n[図1] 時系列プロット作成中...")
ts = df_b.groupby('年度').agg(
    転出率_mean=('転出率', 'mean'),
    着工率_mean=('着工率', 'mean'),
    転入率_mean=('転入率', 'mean'),
).reset_index()

fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
fig.suptitle('転出率・着工新設住宅の全国平均推移（2012〜2023年）', fontsize=14, fontweight='bold', y=1.01)

# 左：転出率
ax = axes[0]
ax.plot(ts['年度'], ts['転出率_mean'], 'o-', color='#C62828', linewidth=2.2,
        markersize=6, label='転出率（千人対）')
ax.plot(ts['年度'], ts['転入率_mean'], 's--', color='#1565C0', linewidth=2,
        markersize=5, label='転入率（千人対）', alpha=0.8)
ax.fill_between(ts['年度'], ts['転出率_mean'], ts['転入率_mean'],
                where=ts['転出率_mean'] > ts['転入率_mean'],
                alpha=0.15, color='#C62828', label='転出超過')
ax.fill_between(ts['年度'], ts['転出率_mean'], ts['転入率_mean'],
                where=ts['転出率_mean'] <= ts['転入率_mean'],
                alpha=0.15, color='#1565C0')
ax.set_xlabel('年度', fontsize=11)
ax.set_ylabel('率（人口千対）', fontsize=11)
ax.set_title('転出率・転入率', fontsize=12)
ax.legend(fontsize=10)
ax.grid(axis='y', alpha=0.4)
ax.xaxis.set_major_locator(mticker.MultipleLocator(2))

# 右：着工率
ax2 = axes[1]
ax2.plot(ts['年度'], ts['着工率_mean'], 'D-', color='#2E7D32', linewidth=2.2,
         markersize=6, label='着工新設住宅戸数（人口1万対）')
ax2.set_xlabel('年度', fontsize=11)
ax2.set_ylabel('戸数（人口1万対）', fontsize=11)
ax2.set_title('着工新設住宅戸数', fontsize=12)
ax2.legend(fontsize=10)
ax2.grid(axis='y', alpha=0.4)
ax2.xaxis.set_major_locator(mticker.MultipleLocator(2))

# 補足テキスト
axes[0].annotate('コロナ禍\n2020年', xy=(2020, ts.loc[ts['年度']==2020,'転出率_mean'].values[0]),
                 xytext=(2021, ts.loc[ts['年度']==2020,'転出率_mean'].values[0]+0.5),
                 fontsize=9, color='gray',
                 arrowprops=dict(arrowstyle='->', color='gray', lw=1))

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2022_U5_1_fig1_ts.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("  -> 保存: 2022_U5_1_fig1_ts.png")

print("\n" + "="*60)
print("[固定効果モデル 主要指標]")
print(f"  R² (within)  : {fe_res.rsquared:.4f}")
print(f"  R² (between) : {fe_res.rsquared_between:.4f}")
print(f"  観測数       : {fe_res.nobs}")
print(f"  エンティティ : {fe_res.entity_info['total']}")
print()
print("[係数（転出率への影響）]")
for var in EXOG_VARS:
    if var in fe_res.params:
        c = fe_res.params[var]
        p = fe_res.pvalues[var]
        sig = '***' if p < 0.001 else '** ' if p < 0.01 else '*  ' if p < 0.05 else '   '
        print(f"  {var:18s}: {c:+.4f}  {sig}  (p={p:.4f})")
print()
print("[Hausman 検定]")
print(f"  H = {H_stat:.4f},  df = {H_df},  p = {H_pval:.4f}")
if H_pval < 0.05:
    print("  → 固定効果モデルを採用（ランダム効果と有意差あり）")
else:
    print("  → p≥0.05（ランダム効果モデルも統計的に許容）")
print()
print("[生成図一覧]")
for f in ['2022_U5_1_fig1_ts.png', '2022_U5_1_fig2_fe_coef.png',
          '2022_U5_1_fig3_scatter.png', '2022_U5_1_fig4_region.png']:
    path = os.path.join(FIG_DIR, f)
    size = os.path.getsize(path) // 1024 if os.path.exists(path) else -1
    print(f"  {f} ({size} KB)")
print("="*60)
print("\n完了: 全4図を html/figures/ に保存しました。")

print("\n[図2] 固定効果係数プロット作成中...")

coef_df = pd.DataFrame({
    'coef': fe_res.params.drop('const', errors='ignore'),
    'se':   fe_res.std_errors.drop('const', errors='ignore'),
    'pval': fe_res.pvalues.drop('const', errors='ignore'),
})
coef_df['ci_lo'] = coef_df['coef'] - 1.96 * coef_df['se']
coef_df['ci_hi'] = coef_df['coef'] + 1.96 * coef_df['se']
coef_df = coef_df.reindex(EXOG_VARS)

var_labels = {
    '高齢化率':    '高齢化率（%）',
    '地価_log':    '住宅地地価\n（log, 円/㎡）',
    '求人倍率':    '有効求人倍率',
    '消費支出_log':'消費支出\n（log, 円）',
}
labels = [var_labels.get(v, v) for v in coef_df.index]
colors = ['#C62828' if p < 0.05 else '#90A4AE' for p in coef_df['pval']]
sig_markers = ['★' if p < 0.01 else '✓' if p < 0.05 else '' for p in coef_df['pval']]

fig, ax = plt.subplots(figsize=(8, 5))
y_pos = np.arange(len(coef_df))
ax.barh(y_pos, coef_df['coef'], color=colors, height=0.55, alpha=0.85)
ax.errorbar(coef_df['coef'], y_pos,
            xerr=1.96 * coef_df['se'],
            fmt='none', color='#333', elinewidth=1.5, capsize=5)
ax.axvline(0, color='black', linewidth=0.9, linestyle='-')

for i, (coef, sig, pval) in enumerate(zip(coef_df['coef'], sig_markers, coef_df['pval'])):
    xoff = coef_df['ci_hi'].iloc[i] + abs(coef_df['ci_hi'].max()) * 0.03
    pstr = f'p={pval:.3f}' if pval >= 0.001 else 'p<0.001'
    ax.text(xoff, i, f'{sig} {pstr}', va='center', fontsize=9,
            color='#C62828' if pval < 0.05 else '#666')

ax.set_yticks(y_pos)
ax.set_yticklabels(labels, fontsize=11)
ax.set_xlabel('推定係数（従属変数：転出率 人口千対）', fontsize=11)
ax.set_title('固定効果モデル（FE）の推定係数\n（エラーバー：95%信頼区間、標準誤差はクラスタリング）',
             fontsize=12, fontweight='bold')
ax.grid(axis='x', alpha=0.4, linestyle='--')

# 凡例
from matplotlib.patches import Patch
legend_els = [
    Patch(facecolor='#C62828', label='p < 0.05（統計的有意）'),
    Patch(facecolor='#90A4AE', label='p ≥ 0.05（非有意）'),
]
ax.legend(handles=legend_els, fontsize=10, loc='lower right')

# Hausman 結果をテキスト注釈
h_text = (f'Hausman検定: H={H_stat:.2f}, p={H_pval:.3f}\n'
          f'→ 固定効果モデル{"適切" if H_pval < 0.05 else "（RE許容）"}')
ax.text(0.02, 0.98, h_text, transform=ax.transAxes,
        fontsize=9, va='top', color='#333',
        bbox=dict(boxstyle='round,pad=0.4', facecolor='#EFF3FF', alpha=0.8))

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2022_U5_1_fig2_fe_coef.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("  -> 保存: 2022_U5_1_fig2_fe_coef.png")

print("\n[図3] 散布図（転出率 vs 住宅地地価）作成中...")

# 2012〜2023 年の平均で集計
df_avg = df_b.groupby('都道府県').agg(
    転出率_avg=('転出率', 'mean'),
    地価_avg=('標準価格（平均価格）（住宅地）', 'mean'),
    地域=('地域', 'first'),
).reset_index()

region_colors = {
    '北海道・東北': '#1565C0',
    '関東':        '#C62828',
    '中部':        '#2E7D32',
    '近畿':        '#E65100',
    '中国・四国':  '#6A1B9A',
    '九州・沖縄':  '#00838F',
}

fig, ax = plt.subplots(figsize=(10, 7))
for region, grp in df_avg.groupby('地域'):
    ax.scatter(grp['地価_avg'], grp['転出率_avg'],
               color=region_colors.get(region, 'gray'),
               s=70, alpha=0.85, label=region, zorder=3)

# 都道府県ラベル（一部省略せず全て）
for _, row in df_avg.iterrows():
    pref_short = row['都道府県'].replace('県','').replace('府','').replace('都','').replace('道','')
    ax.annotate(pref_short, (row['地価_avg'], row['転出率_avg']),
                textcoords='offset points', xytext=(4, 2),
                fontsize=7, alpha=0.85, color='#333')

# OLS 回帰線
x_all = df_avg['地価_avg'].values
y_all = df_avg['転出率_avg'].values
slope, intercept, r_val, p_val, _ = stats.linregress(x_all, y_all)
xline = np.linspace(x_all.min(), x_all.max(), 100)
ax.plot(xline, slope * xline + intercept, '--', color='gray', linewidth=1.5,
        label=f'回帰線 (r={r_val:.2f}, p={p_val:.3f})')

ax.set_xlabel('住宅地平均地価（円/㎡）\n2012〜2023年平均', fontsize=12)
ax.set_ylabel('転出率（人口千対）\n2012〜2023年平均', fontsize=12)
ax.set_title('転出率 vs 住宅地地価（都道府県別）', fontsize=13, fontweight='bold')
ax.legend(fontsize=10, ncol=2, loc='upper right')
ax.grid(alpha=0.3, linestyle='--')
ax.xaxis.set_major_formatter(mticker.FuncFormatter(lambda x, _: f'{int(x):,}'))

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2022_U5_1_fig3_scatter.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("  -> 保存: 2022_U5_1_fig3_scatter.png")

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

from linearmodels.panel import PanelOLS, RandomEffects
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

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

print("\n[パネル分析] 固定効果モデル推定中...")

EXOG_VARS = ['高齢化率', '地価_log', '求人倍率', '消費支出_log']
ENDOG_VAR  = '転出率'

df_panel = df_b[['都道府県', '年度', ENDOG_VAR] + EXOG_VARS].dropna().copy()
print(f"  分析対象: {len(df_panel)} 観測 ({df_panel['都道府県'].nunique()} 都道府県)")

# MultiIndex 設定
df_panel = df_panel.set_index(['都道府県', '年度'])

# ── 固定効果モデル（FE）
fe_model = PanelOLS(
    df_panel[ENDOG_VAR],
    sm.add_constant(df_panel[EXOG_VARS]),
    entity_effects=True,
    time_effects=False,
)
fe_res = fe_model.fit(cov_type='clustered', cluster_entity=True)
print("\n[固定効果モデル推定結果]")
print(fe_res.summary.tables[1])

# ── ランダム効果モデル（RE）— Hausman 検定用
re_model = RandomEffects(
    df_panel[ENDOG_VAR],
    sm.add_constant(df_panel[EXOG_VARS]),
)
re_res = re_model.fit(cov_type='robust')

print("\n[Hausman検定] FE vs RE 係数差の検定...")
fe_params = fe_res.params.drop('const', errors='ignore')
re_params = re_res.params.drop('const', errors='ignore')
common_vars = fe_params.index.intersection(re_params.index)

b_diff = fe_params[common_vars].values - re_params[common_vars].values

# 共分散行列差
fe_cov = fe_res.cov.loc[common_vars, common_vars].values
re_cov = re_res.cov.loc[common_vars, common_vars].values
cov_diff = fe_cov - re_cov

try:
    cov_inv = np.linalg.pinv(cov_diff)
    H_stat = float(b_diff @ cov_inv @ b_diff)
    H_df   = len(common_vars)
    H_pval = 1 - stats.chi2.cdf(H_stat, H_df)
    hausman_ok = True
except Exception:
    H_stat, H_df, H_pval = np.nan, len(common_vars), np.nan
    hausman_ok = False

print(f"  Hausman 統計量 H = {H_stat:.4f}")
print(f"  自由度 = {H_df}")
print(f"  p値 = {H_pval:.4f}")
if hausman_ok:
    if H_pval < 0.05:
        print("  → p<0.05: 固定効果モデルが適切（個体固有の偏りを制御）")
    else:
        print("  → p≥0.05: ランダム効果モデルも許容")

print("\n[図4] 地域別転出率比較作成中...")

# 2022年（最新の完全年）で比較（全期間の平均も補助で）
df_region_data = []
for region in region_order:
    vals = df_b[df_b['地域'] == region]['転出率'].dropna().values
    df_region_data.append(vals)

fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))
fig.suptitle('地域別 転出率の分布（2012〜2023年）', fontsize=14, fontweight='bold')

# 左：箱ひげ図
ax = axes[0]
bp = ax.boxplot(df_region_data, patch_artist=True, notch=False,
                medianprops=dict(color='white', linewidth=2.5),
                whiskerprops=dict(linewidth=1.5),
                flierprops=dict(marker='o', markersize=4, alpha=0.5))
for patch, region in zip(bp['boxes'], region_order):
    patch.set_facecolor(region_colors[region])
    patch.set_alpha(0.8)
ax.set_xticks(range(1, len(region_order)+1))
ax.set_xticklabels([r.replace('・', '\n') for r in region_order], fontsize=9)
ax.set_ylabel('転出率（人口千対）', fontsize=11)
ax.set_title('地域別転出率（箱ひげ図）', fontsize=12)
ax.grid(axis='y', alpha=0.4)

# 右：年別推移（地域別折れ線）
ax2 = axes[1]
for region in region_order:
    ts_r = df_b[df_b['地域'] == region].groupby('年度')['転出率'].mean()
    ax2.plot(ts_r.index, ts_r.values, 'o-', color=region_colors[region],
             linewidth=1.8, markersize=4, label=region, alpha=0.9)

ax2.set_xlabel('年度', fontsize=11)
ax2.set_ylabel('転出率（人口千対）', fontsize=11)
ax2.set_title('地域別 転出率の年次推移', fontsize=12)
ax2.legend(fontsize=9, loc='upper right', ncol=2)
ax2.grid(alpha=0.4)
ax2.xaxis.set_major_locator(mticker.MultipleLocator(2))

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2022_U5_1_fig4_region.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("  -> 保存: 2022_U5_1_fig4_region.png")