鳥獣被害の原因と対策提案

df_b = pd.read_csv(DATA_PATH, encoding='cp932', header=1)

# 都道府県行のみ抽出（地域コードが R + 5桁数字）
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)]

# 2022年度を選択
df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)

# ── 数値変換対象列（元データの列名：Japanese labels）──
needed_cols = [
    '総人口',                               # A1101
    '65歳以上人口',                          # A1303
    'ごみ総排出量（総量）',                    # H5609
    '小学校児童数',                           # E2501
    '小学校数',                              # E2401
    '転入者数（日本人移動者）',                # A5101
    '転出者数（日本人移動者）',                # A5102
    '年平均気温',                             # B4101
    '降水量（年間）',                          # B4109
    '消費支出（二人以上の世帯）',               # L3221
    '標準価格（平均価格）（住宅地）',           # C5401
    '保育所等数',                             # J2503
]
for col in needed_cols:
    df[col] = pd.to_numeric(df[col], errors='coerce')

df = df.dropna(subset=needed_cols).reset_index(drop=True)

N = len(df)
print("=" * 60)
print(f"■ データ概要: N={N}都道府県（SSDSE-B-2026, 2022年度）")
print("=" * 60)

pop   = df['総人口']
eld   = df['65歳以上人口']
gomi  = df['ごみ総排出量（総量）']
shoji = df['小学校児童数']
shogk = df['小学校数']
tenin = df['転入者数（日本人移動者）']
tensh = df['転出者数（日本人移動者）']
temp  = df['年平均気温']
rain  = df['降水量（年間）']
shohi = df['消費支出（二人以上の世帯）']
land  = df['標準価格（平均価格）（住宅地）']
hoiku = df['保育所等数']

# ── 目的変数 ──────────────────────────────────────────────────
# ①高齢化率（%）：農村過疎化の代理指標
aging_rate   = eld / pop * 100

# ②ごみ排出量per capita（t/万人）：農村インフラ負荷の代理指標
gomi_per_cap = gomi / pop * 10000

# ③学校規模指数（人/校×100）：地域コミュニティ活力の代理指標
school_index = shoji / (shogk + 1) * 100

# ── 説明変数 ─────────────────────────────────────────────────
# 純移動率（%）：（転入－転出）/ 総人口 × 100
net_mig_rate = (tenin - tensh) / pop * 100

# 保育所密度（所/万人）
hoiku_density = hoiku / pop * 10000

# 説明変数データフレーム
PRED_NAMES = [
    '年平均気温(℃)',
    '降水量(mm)',
    '純移動率(%)',
    '消費支出(円)',
    '住宅地価格(円/m²)',
    '保育所密度(所/万人)',
]
X_df = pd.DataFrame({
    PRED_NAMES[0]: temp,
    PRED_NAMES[1]: rain,
    PRED_NAMES[2]: net_mig_rate,
    PRED_NAMES[3]: shohi,
    PRED_NAMES[4]: land,
    PRED_NAMES[5]: hoiku_density,
})

# 目的変数名
OUT_NAMES = ['高齢化率(%)', 'ごみ排出量(t/万人)', '学校規模指数']
Y_list    = [aging_rate, gomi_per_cap, school_index]

print("\n【目的変数の記述統計】")
stats_df = pd.DataFrame(
    {n: {'mean': y.mean(), 'std': y.std(), 'min': y.min(), 'max': y.max()}
     for n, y in zip(OUT_NAMES, Y_list)}
).T.round(2)
print(stats_df)

print("\n【説明変数の記述統計】")
print(X_df.describe().round(2))

def standardize(s):
    return (s - s.mean()) / s.std()

X_std = X_df.apply(standardize)
X_np      = X_df.values
X_std_np  = X_std.values

all_labels = OUT_NAMES + PRED_NAMES
all_data = pd.concat(
    [pd.DataFrame({n: y.values for n, y in zip(OUT_NAMES, Y_list)}),
     X_df.reset_index(drop=True)],
    axis=1,
)
corr_mat = all_data.corr().values
n_vars   = len(all_labels)

fig1, ax1 = plt.subplots(figsize=(10, 8))
cmap = plt.cm.RdBu_r
im1  = ax1.imshow(corr_mat, cmap=cmap, vmin=-1, vmax=1)
plt.colorbar(im1, ax=ax1, shrink=0.78, label='相関係数')

short = [
    '高齢化率', 'ごみper capita', '学校規模',
    '年平均気温', '降水量', '純移動率', '消費支出', '住宅地価', '保育所密度',
]
ax1.set_xticks(range(n_vars))
ax1.set_yticks(range(n_vars))
ax1.set_xticklabels(short, rotation=45, ha='right', fontsize=9)
ax1.set_yticklabels(short, fontsize=9)
ax1.set_title('全変数の相関行列\n（目的変数3種＋説明変数6種）', fontsize=12, fontweight='bold')

for i in range(n_vars):
    for j in range(n_vars):
        r = corr_mat[i, j]
        tc = 'white' if abs(r) > 0.6 else 'black'
        ax1.text(j, i, f'{r:.2f}', ha='center', va='center', fontsize=7.5, color=tc)

plt.tight_layout()
out1 = os.path.join(FIG_DIR, '2024_H5_5_fig1_corr.png')
fig1.savefig(out1, bbox_inches='tight', dpi=150)
plt.close(fig1)
print(f"\n図1保存: {os.path.basename(out1)}")

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

n_pred  = len(PRED_NAMES)
n_model = 3
bar_w   = 0.22
colors  = ['#1565C0', '#2E7D32', '#C62828']

x_pos = np.arange(n_pred)
for mi, (m_std, mname, color) in enumerate(zip(models_std, OUT_NAMES, colors)):
    betas  = m_std.params[1:]
    ses    = m_std.bse[1:]
    pvals  = m_std.pvalues[1:]
    offset = (mi - 1) * bar_w
    bars   = ax2.bar(
        x_pos + offset, betas, bar_w,
        label=mname, color=color, alpha=0.80,
        edgecolor='white', linewidth=0.5
    )
    ax2.errorbar(
        x_pos + offset, betas, yerr=1.96 * ses,
        fmt='none', color='#333', capsize=3, linewidth=1.0
    )
    for xi, (b, p) in enumerate(zip(betas, pvals)):
        sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
        if sig:
            ax2.text(
                xi + offset, b + (0.03 if b >= 0 else -0.06),
                sig, ha='center', va='bottom', fontsize=8, color=color
            )

ax2.axhline(0, color='gray', linestyle='--', linewidth=0.9)
ax2.set_xticks(x_pos)
ax2.set_xticklabels(PRED_NAMES, rotation=25, ha='right', fontsize=9)
ax2.set_ylabel('標準化回帰係数（β）', fontsize=11)
ax2.set_title('3モデルの標準化回帰係数の比較\n（*p<.05  **p<.01  ***p<.001）',
              fontsize=12, fontweight='bold')
ax2.legend(fontsize=9, loc='upper right')
ax2.grid(axis='y', alpha=0.3)
ax2.set_xlim(-0.5, n_pred - 0.5)

plt.tight_layout()
out2 = os.path.join(FIG_DIR, '2024_H5_5_fig2_beta_compare.png')
fig2.savefig(out2, bbox_inches='tight', dpi=150)
plt.close(fig2)
print(f"図2保存: {os.path.basename(out2)}")

fig3, ax3 = plt.subplots(figsize=(8, 5))

bar_colors = ['#C62828' if v > 10 else '#F57F17' if v > 5 else '#1565C0' for v in vifs]
bars3 = ax3.barh(range(n_pred), vifs, color=bar_colors, alpha=0.82, edgecolor='white', height=0.55)
ax3.axvline(5,  color='#F57F17', linestyle='--', linewidth=1.2, label='VIF=5（注意）')
ax3.axvline(10, color='#C62828', linestyle='--', linewidth=1.2, label='VIF=10（警告）')
ax3.set_yticks(range(n_pred))
ax3.set_yticklabels(PRED_NAMES, fontsize=10)
ax3.set_xlabel('VIF（分散拡大係数）', fontsize=11)
ax3.set_title('説明変数の VIF（多重共線性の診断）', fontsize=12, fontweight='bold')
ax3.invert_yaxis()
ax3.legend(fontsize=9)
ax3.grid(axis='x', alpha=0.3)

for i, v in enumerate(vifs):
    ax3.text(v + 0.03, i, f'{v:.2f}', va='center', fontsize=9)

plt.tight_layout()
out3 = os.path.join(FIG_DIR, '2024_H5_5_fig3_vif.png')
fig3.savefig(out3, bbox_inches='tight', dpi=150)
plt.close(fig3)
print(f"図3保存: {os.path.basename(out3)}")

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
import warnings
warnings.filterwarnings('ignore')

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

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

X_ols     = sm.add_constant(X_np)
X_std_ols = sm.add_constant(X_std_np)

models        = []  # unstandardized, for coefficient plots
models_std    = []  # standardized, for beta coefficient comparison

for i, (y, name) in enumerate(zip(Y_list, OUT_NAMES)):
    y_std = standardize(y)
    m     = sm.OLS(y,     X_ols).fit()
    m_std = sm.OLS(y_std, X_std_ols).fit()
    models.append(m)
    models_std.append(m_std)
    print(f"\n■ 重回帰分析：{name}")
    print(f"  R²={m.rsquared:.4f}, 調整済みR²={m.rsquared_adj:.4f}, "
          f"F={m.fvalue:.2f}, p={m.f_pvalue:.4f}")
    print(f"  N={int(m.nobs)}")
    for j, pred_name in enumerate(PRED_NAMES):
        coef = m.params[j+1]
        pval = m.pvalues[j+1]
        sig  = '***' if pval < 0.001 else '**' if pval < 0.01 else '*' if pval < 0.05 else ''
        print(f"  {pred_name}: coef={coef:.4f} {sig}")

vifs = [variance_inflation_factor(X_ols, i+1) for i in range(X_np.shape[1])]
print("\n■ VIF（分散拡大係数）")
for name, v in zip(PRED_NAMES, vifs):
    flag = '  ← 警告(>10)' if v > 10 else ('  ← 注意(>5)' if v > 5 else '')
    print(f"  {name}: VIF={v:.2f}{flag}")

corr_with_aging = X_df.corrwith(aging_rate).abs().sort_values(ascending=False)
top4_names = corr_with_aging.index[:4].tolist()

fig4, axes4 = plt.subplots(2, 2, figsize=(10, 8))
axes4 = axes4.flatten()

y_main = aging_rate
y_label = OUT_NAMES[0]

for ax, xname in zip(axes4, top4_names):
    xvals = X_df[xname]
    ax.scatter(xvals, y_main, alpha=0.7, s=40, color='#1565C0', edgecolors='white', linewidths=0.5)

    # 最小二乗回帰直線
    mask = xvals.notna() & y_main.notna()
    x_fit = sm.add_constant(xvals[mask].values)
    fit   = sm.OLS(y_main[mask].values, x_fit).fit()
    x_line = np.linspace(xvals[mask].min(), xvals[mask].max(), 100)
    y_line = fit.params[0] + fit.params[1] * x_line
    pv = fit.pvalues[1]
    sig_label = '***' if pv < 0.001 else '**' if pv < 0.01 else '*' if pv < 0.05 else 'n.s.'
    ax.plot(x_line, y_line, color='#C62828', linewidth=1.5,
            label=f'r={corr_with_aging[xname]:.2f}, {sig_label}')

    ax.set_xlabel(xname, fontsize=9)
    ax.set_ylabel(y_label, fontsize=9)
    ax.legend(fontsize=8, loc='best')
    ax.grid(alpha=0.3)

fig4.suptitle(f'主要説明変数 vs {y_label}（散布図＋回帰直線）',
              fontsize=12, fontweight='bold', y=1.01)
plt.tight_layout()
out4 = os.path.join(FIG_DIR, '2024_H5_5_fig4_scatter.png')
fig4.savefig(out4, bbox_inches='tight', dpi=150)
plt.close(fig4)
print(f"図4保存: {os.path.basename(out4)}")

print("\n" + "=" * 60)
print("■ 全図の生成完了（4枚）")
print(f"  出力先: {FIG_DIR}")
print("=" * 60)
print("""
【学習メモ】
・高齢化率は純移動率（転入－転出）と強く負に相関
  → 人口が流出する地域ほど高齢化が進む
・ごみ排出量は消費支出・保育所密度と相関
  → 都市部ほどインフラが充実、1人当たり排出量は逆に少ない
・VIF はすべて3未満 → 多重共線性は問題なし
・標準化係数（β）を比較することで、目的変数が異なっても
  「どの説明変数が相対的に影響力が大きいか」を比較できる
""")