都道府県別農業生産と食料自給率の
決定要因分析

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

PREFIX = '2021_H3'

def save_fig(n):
    path = os.path.join(FIG_DIR, f'{PREFIX}_fig{n}.png')
    plt.savefig(path, bbox_inches='tight', dpi=150)
    plt.close()
    print(f'  → {path} 保存完了')

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

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)

# 数値列を変換
num_cols = ['総人口', '65歳以上人口', '15歳未満人口', '15～64歳人口',
            '年平均気温', '消費支出（二人以上の世帯）', '食料費（二人以上の世帯）',
            '一般旅券発行件数']
for c in num_cols:
    if c in df_b.columns:
        df_b[c] = pd.to_numeric(df_b[c], errors='coerce')

# 高齢化率を計算
df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100

print(f'  SSDSE-B: {df_b.shape[0]}行 × {df_b.shape[1]}列')
print(f'  年度: {sorted(df_b["年度"].unique())}')

df_e_raw = pd.read_csv(DATA_E, encoding='cp932', header=None)
col_names = df_e_raw.iloc[2].tolist()
df_e = df_e_raw.iloc[3:].copy()
df_e.columns = col_names
df_e = df_e[df_e['地域コード'].astype(str).str.match(r'^R\d{5}$', na=False)].copy()
df_e = df_e[df_e['地域コード'] != 'R00000'].copy()
df_e = df_e.reset_index(drop=True)

# 数値変換
for c in ['農家数（販売農家）', '農家数（自給的農家）', '耕地面積',
          '総面積（北方地域及び竹島を除く）', '可住地面積',
          '総人口', '65歳以上人口', '県内総生産額（平成27年基準）',
          '従業者数（民営）']:
    if c in df_e.columns:
        df_e[c] = pd.to_numeric(df_e[c], errors='coerce')

print(f'  SSDSE-E: {df_e.shape[0]}行 × {df_e.shape[1]}列')

df2022 = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)

# 食料費比率（食料費/消費支出）
df2022['食料費比率'] = pd.to_numeric(df2022['食料費（二人以上の世帯）'], errors='coerce') / \
                       pd.to_numeric(df2022['消費支出（二人以上の世帯）'], errors='coerce') * 100

# 地域コードの末尾3桁でマッチング
# SSDSE-B: R01000, SSDSE-E: R01000 → 同形式
df_merge = df2022.merge(
    df_e[['地域コード', '農家数（販売農家）', '農家数（自給的農家）', '耕地面積',
          '総面積（北方地域及び竹島を除く）', '可住地面積']],
    on='地域コード', how='inner'
)

print(f'  マージ後: {df_merge.shape[0]}都道府県')

# ── 農業関連指標の計算 ──
# 農家密度（販売農家数 / 可住地面積 [ha → km²=ha/100], 戸/km²）
# 可住地面積はha単位
df_merge['農家密度'] = df_merge['農家数（販売農家）'] / (df_merge['可住地面積'] / 100)

# 耕地率（耕地面積[ha] / 総面積[ha] × 100）
# 総面積はha単位
df_merge['耕地率'] = df_merge['耕地面積'] / df_merge['総面積（北方地域及び竹島を除く）'] * 100

# 人口密度（人/km²）：総面積はha → km²=÷100
df_merge['人口密度'] = df_merge['総人口'] / (df_merge['総面積（北方地域及び竹島を除く）'] / 100)

# 地域列を付与
df_merge['地域'] = df_merge['都道府県'].map(region_map)
df_e['地域'] = df_e['都道府県'].map(region_map)

# 分析用データ確認
print('\n=== 断面データ（2022年度）統計量 ===')
for col in ['農家密度', '耕地率', '高齢化率', '人口密度', '食料費比率', '年平均気温']:
    if col in df_merge.columns:
        s = df_merge[col].describe()
        print(f'  {col}: mean={s["mean"]:.2f}, std={s["std"]:.2f}, '
              f'min={s["min"]:.2f}, max={s["max"]:.2f}')

print('\n=== 図2: 農家密度 vs 人口密度 散布図 ===')

df_plot = df_merge[['都道府県', '地域', '農家密度', '人口密度', '耕地率']].dropna()

fig, ax = plt.subplots(figsize=(10, 7))

# 回帰線
x_all = df_plot['人口密度'].values
y_all = df_plot['農家密度'].values
slope, intercept, r_val, p_val, se = stats.linregress(x_all, y_all)
x_line = np.linspace(x_all.min(), x_all.max(), 200)
ax.plot(x_line, slope * x_line + intercept,
        color='#333333', linestyle='--', linewidth=1.5,
        label=f'回帰直線 (r={r_val:.3f}, p={p_val:.3f})')

# 散布図
for region in REGIONS:
    mask = df_plot['地域'] == region
    if mask.sum() == 0:
        continue
    ax.scatter(df_plot.loc[mask, '人口密度'],
               df_plot.loc[mask, '農家密度'],
               color=region_colors[region], s=70,
               alpha=0.85, edgecolors='white', linewidth=0.5,
               label=region)

# 都道府県ラベル（全47）
for _, row in df_plot.iterrows():
    ax.annotate(row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', ''),
                xy=(row['人口密度'], row['農家密度']),
                fontsize=7, color='#333333',
                xytext=(3, 3), textcoords='offset points')

ax.set_xlabel('人口密度（人/km²）', fontsize=12)
ax.set_ylabel('農家密度（販売農家数/可住地面積, 戸/km²）', fontsize=12)
ax.set_title('農家密度 × 人口密度（47都道府県）\n'
             '— 都市化が進むほど農家密度は低下 —', fontsize=13)
ax.legend(loc='upper right', fontsize=9, ncol=2)
ax.grid(True, alpha=0.3)

# 相関係数注釈
ax.text(0.03, 0.97,
        f'r = {r_val:.3f}  (p = {p_val:.4f})\nデータ: SSDSE-B 2022 + SSDSE-E 2019/2024',
        transform=ax.transAxes, fontsize=9,
        verticalalignment='top',
        bbox=dict(boxstyle='round,pad=0.4', facecolor='#EFF3FF', alpha=0.8))

fig.tight_layout()
save_fig(2)

print(f'  農家密度 × 人口密度: r={r_val:.3f}, p={p_val:.4f}')

print('\n' + '='*60)
print('=== HTML用 統計値まとめ ===')
print('='*60)

# 高齢化率（2022年、地域別）
print('\n【地域別高齢化率（2022年）】')
df_b22 = df_b[df_b['年度'] == 2022].copy()
df_b22['地域'] = df_b22['都道府県'].map(region_map)
for region in REGIONS:
    mask = df_b22['地域'] == region
    mean_val = df_b22.loc[mask, '高齢化率'].mean()
    print(f'  {region}: {mean_val:.1f}%')

print(f'\n  全国平均高齢化率（2022）: {df_b22["高齢化率"].mean():.1f}%')
print(f'  高齢化率最高: {df_b22.loc[df_b22["高齢化率"].idxmax(), "都道府県"]} '
      f'({df_b22["高齢化率"].max():.1f}%)')
print(f'  高齢化率最低: {df_b22.loc[df_b22["高齢化率"].idxmin(), "都道府県"]} '
      f'({df_b22["高齢化率"].min():.1f}%)')

print('\n【OLS回帰分析 主要結果】')
print(f'  N={len(df_ols)}, R²={model_raw.rsquared:.3f}, adj.R²={model_raw.rsquared_adj:.3f}')
print(f'  F={model_raw.fvalue:.2f}, p={model_raw.f_pvalue:.4f}')

print('\n【農家密度 Top5・Bottom5】')
print('  Top5:')
for i in range(min(5, len(df_rank))):
    r = df_rank.iloc[i]
    print(f'    {i+1}. {r["都道府県"]}: {r["農家密度"]:.3f} 戸/km²')
print('  Bottom5:')
for i in range(len(df_rank)-5, len(df_rank)):
    r = df_rank.iloc[i]
    print(f'    {i+1}. {r["都道府県"]}: {r["農家密度"]:.3f} 戸/km²')

print('\n【農家密度 × 人口密度 相関】')
print(f'  r = {r_val:.3f}, p = {p_val:.4f}')

print('\n=== 全図保存完了 ===')

import os
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
import warnings
warnings.filterwarnings('ignore')

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

print('\n=== 図1: 地域別高齢化率の推移 ===')

years = sorted(df_b['年度'].unique())

fig, ax = plt.subplots(figsize=(10, 6))

for region in REGIONS:
    prefs = [p for p, r in region_map.items() if r == region]
    region_data = df_b[df_b['都道府県'].isin(prefs)].groupby('年度')['高齢化率'].mean()
    region_data = region_data.reindex(years)
    ax.plot(years, region_data.values,
            color=region_colors[region], linewidth=2.2,
            marker='o', markersize=4, label=region)

ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('高齢化率（65歳以上人口比率, %）', fontsize=12)
ax.set_title('地域別 高齢化率の推移（2012〜2022年）\n'
             '— 農業地域における後継者問題の背景 —', fontsize=13)
ax.legend(loc='upper left', fontsize=10)
ax.set_xticks(years)
ax.set_xticklabels([str(y) for y in years], rotation=45)
ax.grid(True, alpha=0.35)

# 農業高齢化の注釈
ax.annotate('農村部（北海道・東北、中国・四国）\nで高齢化が加速',
            xy=(2022, 33), xytext=(2018, 36),
            fontsize=9, color='#4e9af1',
            arrowprops=dict(arrowstyle='->', color='#4e9af1', lw=1.5))

fig.tight_layout()
save_fig(1)

print('\n=== OLS 重回帰分析 ===')

df_ols = df_merge[['都道府県', '農家密度', '耕地率', '高齢化率',
                   '人口密度', '食料費比率', '年平均気温']].dropna().copy()

df_ols['ln_人口密度'] = np.log(df_ols['人口密度'] + 1)

y_ols  = df_ols['農家密度']
X_vars = ['耕地率', '高齢化率', 'ln_人口密度', '食料費比率', '年平均気温']
X_raw  = df_ols[X_vars]

# 標準化（標準化偏回帰係数のため）
X_std_vals = (X_raw - X_raw.mean()) / X_raw.std()
y_std_vals = (y_ols - y_ols.mean()) / y_ols.std()

X_std_c = sm.add_constant(X_std_vals)
model_std = sm.OLS(y_std_vals, X_std_c).fit()

X_c = sm.add_constant(X_raw)
model_raw = sm.OLS(y_ols, X_c).fit()

print(f'  N={len(df_ols)}, R²={model_raw.rsquared:.3f}, adj.R²={model_raw.rsquared_adj:.3f}')
print(f'  F統計量={model_raw.fvalue:.2f}, p={model_raw.f_pvalue:.4f}')
print()

var_labels = {
    '耕地率':      '耕地率（%）',
    '高齢化率':    '高齢化率（%）',
    'ln_人口密度': 'ln（人口密度）',
    '食料費比率':  '食料費比率（%）',
    '年平均気温':  '年平均気温（℃）'
}

print('  変数           | 標準化係数  |  p値')
print('  ' + '-'*45)
for var in X_vars:
    coef = model_std.params[var]
    pval = model_std.pvalues[var]
    sig  = '**' if pval < 0.01 else ('*' if pval < 0.05 else '')
    label = var_labels.get(var, var)
    print(f'  {label:<20s} | {coef:+.4f}  |  {pval:.4f} {sig}')

print('\n=== 図3: 標準化偏回帰係数プロット ===')

coefs = model_std.params[X_vars].values
pvals = model_std.pvalues[X_vars].values
ci_lo = model_std.conf_int().loc[X_vars, 0].values
ci_hi = model_std.conf_int().loc[X_vars, 1].values
labels = [var_labels.get(v, v) for v in X_vars]

# 係数の大きさでソート
order = np.argsort(coefs)
coefs_s  = coefs[order]
pvals_s  = pvals[order]
labels_s = [labels[i] for i in order]
ci_lo_s  = ci_lo[order]
ci_hi_s  = ci_hi[order]
xerr_lo  = coefs_s - ci_lo_s
xerr_hi  = ci_hi_s - coefs_s

colors_bar = ['#e05c5c' if p < 0.05 else '#aaaaaa' for p in pvals_s]

fig, ax = plt.subplots(figsize=(9, 5))
bars = ax.barh(range(len(labels_s)), coefs_s, color=colors_bar,
               edgecolor='white', height=0.55)
ax.errorbar(coefs_s, range(len(labels_s)),
            xerr=[xerr_lo, xerr_hi],
            fmt='none', color='#444444', capsize=4, linewidth=1.5)
ax.axvline(0, color='#333333', linewidth=1.2)

for i, (c, p) in enumerate(zip(coefs_s, pvals_s)):
    sig_str = '**' if p < 0.01 else ('*' if p < 0.05 else 'n.s.')
    ax.text(c + (0.01 if c >= 0 else -0.01), i,
            f'{c:+.3f} ({sig_str})',
            va='center', ha='left' if c >= 0 else 'right',
            fontsize=9, color='#222222')

ax.set_yticks(range(len(labels_s)))
ax.set_yticklabels(labels_s, fontsize=11)
ax.set_xlabel('標準化偏回帰係数（β）', fontsize=12)
ax.set_title('農家密度の決定要因：標準化偏回帰係数\n'
             f'R²={model_raw.rsquared:.3f}, adj.R²={model_raw.rsquared_adj:.3f}, '
             f'N={len(df_ols)}',
             fontsize=13)

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

fig.tight_layout()
save_fig(3)

print('\n=== 図4: 農家密度ランキング ===')

df_rank = df_merge[['都道府県', '地域', '農家密度', '耕地率']].dropna().copy()
df_rank = df_rank.sort_values('農家密度', ascending=False).reset_index(drop=True)

fig, ax = plt.subplots(figsize=(14, 7))

bar_colors = [region_colors.get(r, '#888888') for r in df_rank['地域']]
ax.bar(range(len(df_rank)), df_rank['農家密度'], color=bar_colors,
       edgecolor='white', linewidth=0.5)

ax.set_xticks(range(len(df_rank)))
ax.set_xticklabels(
    [p.replace('県', '').replace('府', '').replace('都', '').replace('道', '')
     for p in df_rank['都道府県']],
    rotation=75, ha='right', fontsize=8.5
)
ax.set_ylabel('農家密度（販売農家数/可住地面積, 戸/km²）', fontsize=11)
ax.set_title('都道府県別 農家密度ランキング\n'
             '— 農業生産集積の地域差（SSDSE-E 2019農業センサス + 面積データ）—',
             fontsize=13)

# 地域凡例
patches = [mpatches.Patch(color=region_colors[r], label=r) for r in REGIONS]
ax.legend(handles=patches, loc='upper right', fontsize=9, ncol=2)
ax.grid(True, axis='y', alpha=0.3)

# 上位・下位のラベル
for i in range(min(5, len(df_rank))):
    ax.text(i, df_rank['農家密度'].iloc[i] + 0.002,
            f"{df_rank['農家密度'].iloc[i]:.3f}",
            ha='center', va='bottom', fontsize=7.5, color='#333333')
for i in range(len(df_rank)-5, len(df_rank)):
    ax.text(i, df_rank['農家密度'].iloc[i] + 0.002,
            f"{df_rank['農家密度'].iloc[i]:.3f}",
            ha='center', va='bottom', fontsize=7.5, color='#333333')

fig.tight_layout()
save_fig(4)

print(f'\n  全国農家密度 平均: {df_rank["農家密度"].mean():.3f} 戸/km²')
print(f'  最高: {df_rank.iloc[0]["都道府県"]} ({df_rank.iloc[0]["農家密度"]:.3f})')
print(f'  最低: {df_rank.iloc[-1]["都道府県"]} ({df_rank.iloc[-1]["農家密度"]:.3f})')