""" 教育用再現コード: 2021年度 統計数理賞（高校生の部） ================================================================= 論文タイトル：都道府県別農業生産と食料自給率の決定要因分析 受賞：統計数理賞（高校生の部） 【分析概要】 データ：SSDSE-B-2026（47都道府県・2012〜2022年度） SSDSE-E-2026（47都道府県・農家数・耕地面積等） 食料自給率と農業生産の地域差に着目し、重回帰分析・相関分析・時系列分析で 決定要因を明らかにする。SSDSE-Bには直接の農業就業者数データが含まれないため、 SSDSE-Eの農家数（販売農家）・耕地面積と、SSDSE-Bの高齢化率・人口密度等を 組み合わせて代理指標を設計する。 【変数設計（代理指標）】 - 農家密度 (販売農家/可住地面積, 戸/km²) ：農業生産の地域集積を表す主指標 - 耕地率 (耕地面積/総面積, %) ：農地の面積的広がり - 高齢化率 (65歳以上人口/総人口, %) ：農業従事者の高齢化・後継者問題の代理 - 人口密度 (総人口/総面積, 人/km²) ：都市化度・農業地域かどうかの代理 - 食料費比率(食料費/消費支出, %) ：食生活の地域特性、農村的生活の代理 - 気温 (年平均気温, ℃) ：農業生産適地の気候条件 【時系列分析（図1）】 SSDSE-B 2012〜2022年の高齢化率推移を地域別に可視化 （農業地域の高齢化加速を示す代理時系列として活用） 【学習ポイント 4項目】 1. 代理変数の設計：直接データがない場合の変数選択 2. 重回帰分析（OLS）：statsmodelsによる実装と診断 3. 標準化偏回帰係数：偏回帰係数の大きさ比較 4. 相関分析と散布図：2変数関係の可視化 【使用データ】 SSDSE-B-2026.csv（SSDSE: 社会・人口統計体系データセット） SSDSE-E-2026.csv（SSDSE: 都道府県比較データセット） ※ 合成データ・乱数生成（np.random.seed等）は一切使用しない。 ================================================================= """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・SSDSE-B-2026.csv # → data/raw/SSDSE-B-2026.csv に配置 # ・SSDSE-E-2026.csv # → data/raw/SSDSE-E-2026.csv に配置 # # ダウンロード先: # https://www.nstac.go.jp/use/literacy/ssdse/ # （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード） # （SSDSE-E（社会・人口統計体系 都道府県の指標2） の CSV をダウンロード） # # フォルダ配置（プロジェクトルートからの相対パス）: # code/ ← このスクリプトの場所 # data/raw/ ← CSV ファイルをここに配置 # html/figures/ ← 図の出力先（自動生成） # # 実行方法（ファイルを一切編集せず実行可能）: # python3 code/2021_H3_suri.py # ============================================================ 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) # ────────────────────────────────────────────────────────────── # 共通設定 # ────────────────────────────────────────────────────────────── 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()) # ────────────────────────────────────────────────────────────── # データ読み込み：SSDSE-B（2012〜2022年、47都道府県） # ────────────────────────────────────────────────────────────── 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())}') # ────────────────────────────────────────────────────────────── # データ読み込み：SSDSE-E（断面データ、農業統計） # ────────────────────────────────────────────────────────────── 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]}列') # ────────────────────────────────────────────────────────────── # 断面データ作成（2022年 SSDSE-B + SSDSE-E） # ────────────────────────────────────────────────────────────── 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}') # ────────────────────────────────────────────────────────────── # 図1: 時系列グラフ — 地域別高齢化率の推移（2012〜2022年） # ────────────────────────────────────────────────────────────── 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) # ────────────────────────────────────────────────────────────── # 図2: 散布図 — 農家密度 vs 人口密度（47都道府県、地域色分け） # ────────────────────────────────────────────────────────────── 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}') # ────────────────────────────────────────────────────────────── # OLS重回帰分析 # 目的変数 Y : 農家密度（農業生産性の代理） # 説明変数 : 耕地率、高齢化率、人口密度（log）、食料費比率、年平均気温 # ────────────────────────────────────────────────────────────── 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}') # ────────────────────────────────────────────────────────────── # 図3: 標準化偏回帰係数プロット # ────────────────────────────────────────────────────────────── 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) # ────────────────────────────────────────────────────────────── # 図4: 棒グラフ — 都道府県別農家密度ランキング（地域色分け） # ────────────────────────────────────────────────────────────── 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})') # ────────────────────────────────────────────────────────────── # 統計量まとめ出力（HTMLへ反映用） # ────────────────────────────────────────────────────────────── 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=== 全図保存完了 ===')