""" 教育用再現コード: 2022年度 統計データ分析コンペティション 審査員奨励賞（高校生の部） ================================================================= 論文タイトル：都市農村間の所得格差：消費支出と雇用環境の地域差分析 【分析概要】 データ：SSDSE-B-2026.csv（都道府県データ、2012〜2023年） テーマ：都市部と地方農村部の所得水準（消費支出で代理）の格差を 雇用環境・住宅地価・人口動態の観点から分析する。 主要変数： - 消費支出_万円 = 消費支出（二人以上の世帯）/ 10000 ← 所得水準の代理 - 求人倍率 = 月間有効求人数（一般）/ 月間有効求職者数（一般） - 住宅地価_log = log(標準価格（平均価格）（住宅地）) ← 都市度の代理 - 高齢化率 = 65歳以上人口 / 総人口 × 100 - 大学学生数_千人 = 大学学生数 / 1000 - 純移動率 = (転入者数 - 転出者数) / 総人口 × 1000 都市/地方分類（2022年）: 住宅地価が全国中央値以上 → '都市部' 住宅地価が全国中央値未満 → '地方' 【図の構成】 Step1. 消費支出の時系列推移（地域別平均, 2012〜2023年） Step2. 住宅地価（都市度代理）vs 消費支出 散布図（2022年, 都道府県ラベル） Step3. OLS回帰係数プロット（消費支出の決定要因） Step4. 都市部 vs 地方 消費支出比較（箱ひげ図 + Mann-Whitney U検定） 【変数定義（SSDSE-B-2026列名）】 SSDSE-B-2026 : 年度 Prefecture : 都道府県 A1101 : 総人口 A1303 : 65歳以上人口 A5101 : 転入者数（日本人移動者） A5102 : 転出者数（日本人移動者） C5401 : 標準価格（平均価格）（住宅地）（円/m²） E6302 : 大学学生数 F3102 : 月間有効求職者数（一般） F3103 : 月間有効求人数（一般） L3221 : 消費支出（二人以上の世帯）（円/月） 【派生変数の計算式】 消費支出_万円 = L3221 / 10000 求人倍率 = F3103 / F3102 住宅地価_log = log(C5401) 高齢化率 = A1303 / A1101 × 100 大学学生数_千人 = E6302 / 1000 純移動率 = (A5101 - A5102) / A1101 × 1000 【データ出典】 SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ） 統計数理研究所・統計センターより公表の実データ 【データサイエンス学習ポイント】 1. 所得代理変数としての消費支出（エンゲルの法則と恒常所得仮説） 2. Mann-Whitney U検定（正規性が不明な場合の2群比較） 3. 住宅地価を都市度指標として使う根拠（ヘドニック価格法の考え方） 4. 格差の測定（ジニ係数とその計算） ================================================================= """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・SSDSE-B-2026.csv # → data/raw/SSDSE-B-2026.csv に配置 # # ダウンロード先: # https://www.nstac.go.jp/use/literacy/ssdse/ # （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード） # # フォルダ配置（プロジェクトルートからの相対パス）: # code/ ← このスクリプトの場所 # data/raw/ ← CSV ファイルをここに配置 # html/figures/ ← 図の出力先（自動生成） # # 実行方法（ファイルを一切編集せず実行可能）: # python3 code/2022_H5_14_shorei.py # ============================================================ 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 import warnings warnings.filterwarnings('ignore') # ────────────────────────────────────────────────────────────── # 共通設定 # ────────────────────────────────────────────────────────────── 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) # ================================================================ # ■ 実データ読み込み（SSDSE-B-2026: 都道府県データ） # ================================================================ df_raw = pd.read_csv(DATA_B, encoding='shift-jis', header=0) # 行0がラベル行、行1以降が実データ df_all = df_raw.iloc[1:].copy() df_all.columns = df_raw.columns # 数値列に変換 NUM_COLS = ['A1101', 'A1303', 'A5101', 'A5102', 'C5401', 'E6302', 'F3102', 'F3103', 'L3221'] for c in NUM_COLS: df_all[c] = pd.to_numeric(df_all[c], errors='coerce') df_all['年度'] = df_all['SSDSE-B-2026'].astype(int) df_all['都道府県'] = df_all['Prefecture'].astype(str) print("=" * 65) print("■ SSDSE-B-2026 読み込み完了") print(f" 年度範囲: {df_all['年度'].min()}〜{df_all['年度'].max()}") print(f" 都道府県数（各年）: {df_all.groupby('年度')['都道府県'].count().iloc[0]}件") print("=" * 65) # ================================================================ # ■ 特徴量エンジニアリング（全年度） # ================================================================ df_all['消費支出_万円'] = df_all['L3221'] / 10000 df_all['求人倍率'] = df_all['F3103'] / df_all['F3102'] df_all['住宅地価_log'] = np.log(df_all['C5401']) df_all['高齢化率'] = df_all['A1303'] / df_all['A1101'] * 100 df_all['大学学生数_千人'] = df_all['E6302'] / 1000 df_all['純移動率'] = (df_all['A5101'] - df_all['A5102']) / df_all['A1101'] * 1000 # 2022年データ抽出 df_2022 = df_all[df_all['年度'] == 2022].copy().reset_index(drop=True) # 都市部 / 地方 分類（住宅地価の全国中央値を閾値） land_median = df_2022['C5401'].median() df_2022['地域区分'] = df_2022['C5401'].apply( lambda x: '都市部' if x >= land_median else '地方' ) df_all['地域区分'] = df_all.apply( lambda row: '都市部' if row['C5401'] >= land_median else '地方', axis=1 ) print(f"\n■ 2022年データ: {len(df_2022)}都道府県") print(f" 住宅地価 中央値（分類閾値）: {land_median:,.0f} 円/m²") print(f" 都市部: {(df_2022['地域区分']=='都市部').sum()}都道府県") print(f" 地方 : {(df_2022['地域区分']=='地方').sum()}都道府県") print(f" 消費支出（平均）: {df_2022['消費支出_万円'].mean():.2f}万円/月") print(f" 求人倍率 （平均）: {df_2022['求人倍率'].mean():.3f}") print(f" 高齢化率 （平均）: {df_2022['高齢化率'].mean():.1f}%") # ================================================================ # ■ 地方区分マップ（時系列表示用） # ================================================================ REGION_MAP = { '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北', '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北', '福島県': '北海道・東北', '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東', '千葉県': '関東', '東京都': '関東', '神奈川県': '関東', '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部', '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部', '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿', '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿', '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国', '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国', '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国', '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄', '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄', '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄', } REGION_ORDER = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄'] REGION_COLORS = { '北海道・東北': '#1565C0', '関東': '#C62828', '中部': '#2E7D32', '近畿': '#6A1B9A', '中国・四国': '#795548', '九州・沖縄': '#00695C', } df_all['地方'] = df_all['都道府県'].map(REGION_MAP) df_2022['地方'] = df_2022['都道府県'].map(REGION_MAP) # ================================================================ # ■ OLS回帰分析（2022年, 消費支出の決定要因） # ================================================================ print("\n" + "=" * 65) print("■ OLS重回帰分析（2022年）") print("=" * 65) REG_VARS = ['求人倍率', '住宅地価_log', '高齢化率', '大学学生数_千人', '純移動率'] REG_LABELS = ['求人倍率', '住宅地価（log）', '高齢化率（%）', '大学学生数（千人）', '純移動率（‰）'] df_reg = df_2022[['消費支出_万円'] + REG_VARS].dropna() y_ols = df_reg['消費支出_万円'].values X_ols = df_reg[REG_VARS].values # 標準化（係数の大きさを比較するため） X_std = (X_ols - X_ols.mean(axis=0)) / X_ols.std(axis=0) X_with_const = sm.add_constant(X_std) ols_model = sm.OLS(y_ols, X_with_const).fit() print(ols_model.summary2()) print(f"\n R² = {ols_model.rsquared:.4f}") print(f" 自由度修正済R² = {ols_model.rsquared_adj:.4f}") # 個別相関分析 r_land, p_land = stats.pearsonr(df_2022['住宅地価_log'].dropna(), df_2022['消費支出_万円'].dropna()) r_job, p_job = stats.pearsonr(df_2022['求人倍率'].dropna(), df_2022['消費支出_万円'].dropna()) print(f"\n■ 相関分析（住宅地価_log vs 消費支出）: r={r_land:.4f}, p={p_land:.4f}") print(f"■ 相関分析（求人倍率 vs 消費支出） : r={r_job:.4f}, p={p_job:.4f}") # Mann-Whitney U検定（都市部 vs 地方 消費支出） urban = df_2022[df_2022['地域区分']=='都市部']['消費支出_万円'].dropna() rural = df_2022[df_2022['地域区分']=='地方' ]['消費支出_万円'].dropna() u_stat, u_pval = stats.mannwhitneyu(urban, rural, alternative='two-sided') print(f"\n■ Mann-Whitney U検定（都市部 vs 地方 消費支出）:") print(f" U={u_stat:.1f}, p={u_pval:.4f}") print(f" 都市部: n={len(urban)}, 中央値={urban.median():.2f}万円") print(f" 地方 : n={len(rural)}, 中央値={rural.median():.2f}万円") # ジニ係数の計算 def gini(arr): """ジニ係数を計算する（実データ用）""" a = np.sort(arr) n = len(a) return (2 * np.sum(np.arange(1, n + 1) * a) - (n + 1) * np.sum(a)) / (n * np.sum(a)) g_all = gini(df_2022['消費支出_万円'].dropna().values) g_urban = gini(urban.values) g_rural = gini(rural.values) print(f"\n■ ジニ係数（消費支出, 2022年）:") print(f" 全国: {g_all:.4f}") print(f" 都市部のみ: {g_urban:.4f}") print(f" 地方のみ : {g_rural:.4f}") # ================================================================ # ■ 図1: 消費支出の時系列推移（地域別平均, 2012-2023） # ================================================================ print("\n" + "=" * 65) print("■ 図1: 消費支出の時系列推移（地域別平均）") print("=" * 65) df_region_ts = ( df_all.groupby(['年度', '地方'])['消費支出_万円'] .mean() .reset_index() ) fig1, ax1 = plt.subplots(figsize=(11, 6)) for region in REGION_ORDER: sub = df_region_ts[df_region_ts['地方'] == region].sort_values('年度') if len(sub) == 0: continue ax1.plot(sub['年度'], sub['消費支出_万円'], marker='o', linewidth=2.2, markersize=5.5, color=REGION_COLORS[region], label=region) # 2020年にCOVID-19 影響を示す縦線 ax1.axvline(2020, color='#FF8F00', linestyle=':', linewidth=1.5, alpha=0.8) ax1.text(2020.08, ax1.get_ylim()[0] + 0.5, 'COVID-19\n(2020)', ha='left', va='bottom', fontsize=8.5, color='#FF8F00') # 2022年に分析基準年の縦線 ax1.axvline(2022, color='#C62828', linestyle='--', linewidth=1.5, alpha=0.8) ax1.text(2022.08, ax1.get_ylim()[1] - 0.8, '2022年\n（分析基準年）', ha='left', va='top', fontsize=8.5, color='#C62828', fontweight='bold') ax1.set_xlabel('年度', fontsize=12) ax1.set_ylabel('消費支出（万円/月）', fontsize=12) ax1.set_title('消費支出の時系列推移（地方別平均, 2012〜2023年）\n〜関東が最高、九州・東北が低水準で推移〜', fontsize=13, fontweight='bold') ax1.set_xticks(sorted(df_all['年度'].unique())) ax1.xaxis.set_tick_params(rotation=45) ax1.legend(loc='lower right', fontsize=9.5, framealpha=0.88) ax1.grid(True, alpha=0.3) plt.tight_layout() fig1.savefig(os.path.join(FIG_DIR, '2022_H5_14_fig1_timeseries.png'), bbox_inches='tight', dpi=150) plt.close(fig1) print("図1保存: 2022_H5_14_fig1_timeseries.png") # ================================================================ # ■ 図2: 住宅地価 vs 消費支出 散布図（2022年, 都道府県ラベル） # ================================================================ print("\n■ 図2: 住宅地価 vs 消費支出 散布図（2022年）") fig2, ax2 = plt.subplots(figsize=(11, 7)) color_map = {'都市部': '#C62828', '地方': '#1565C0'} marker_map = {'都市部': 'o', '地方': '^'} for grp in ['都市部', '地方']: sub2 = df_2022[df_2022['地域区分'] == grp] ax2.scatter(sub2['住宅地価_log'], sub2['消費支出_万円'], c=color_map[grp], marker=marker_map[grp], s=80, alpha=0.80, edgecolors='#333', linewidth=0.4, label=grp, zorder=3) # 回帰直線（全データ） valid2 = df_2022[['住宅地価_log', '消費支出_万円']].dropna() x2v = valid2['住宅地価_log'].values y2v = valid2['消費支出_万円'].values z2 = np.polyfit(x2v, y2v, 1) xs2 = np.linspace(x2v.min(), x2v.max(), 100) ax2.plot(xs2, np.poly1d(z2)(xs2), 'k--', linewidth=1.8, alpha=0.65, label=f'回帰直線 (r={r_land:.3f}, p={p_land:.4f})') # 中央値の閾値ライン log_med = np.log(land_median) ax2.axvline(log_med, color='#FF8F00', linestyle=':', linewidth=1.5, alpha=0.8) ax2.text(log_med + 0.02, ax2.get_ylim()[0] + 0.3, f'中央値\n({land_median:,.0f}円/m²)', ha='left', va='bottom', fontsize=8, color='#FF8F00') # 都道府県ラベル（注目都市） LABEL_PREFS = { '東京都', '北海道', '沖縄県', '秋田県', '青森県', '神奈川県', '愛知県', '大阪府', '京都府', '福岡県', '高知県', '島根県', '鳥取県', '宮崎県', } for _, row2 in df_2022.iterrows(): if row2['都道府県'] in LABEL_PREFS: short = row2['都道府県'].replace('県', '').replace('都', '').replace('道', '').replace('府', '') ax2.annotate(short, (row2['住宅地価_log'], row2['消費支出_万円']), xytext=(4, 3), textcoords='offset points', fontsize=7.5, color='#1A237E', fontweight='bold') ax2.set_xlabel('住宅地価（log 変換, 円/m²）', fontsize=12) ax2.set_ylabel('消費支出（万円/月）', fontsize=12) ax2.set_title('住宅地価（都市度の代理）と消費支出の関係（2022年, 都道府県別）\n' '〜都市部（赤○）は高地価・高消費支出, 地方（青△）は低水準〜', fontsize=13, fontweight='bold') ax2.legend(fontsize=9.5) ax2.grid(True, alpha=0.3) plt.tight_layout() fig2.savefig(os.path.join(FIG_DIR, '2022_H5_14_fig2_scatter.png'), bbox_inches='tight', dpi=150) plt.close(fig2) print("図2保存: 2022_H5_14_fig2_scatter.png") # ================================================================ # ■ 図3: OLS回帰係数プロット（消費支出の決定要因） # ================================================================ print("\n■ 図3: OLS回帰係数プロット") coefs3 = np.asarray(ols_model.params[1:]) ses3 = np.asarray(ols_model.bse[1:]) pvals3 = np.asarray(ols_model.pvalues[1:]) COEF_COLORS = [] for p in pvals3: if p < 0.01: COEF_COLORS.append('#C62828') elif p < 0.05: COEF_COLORS.append('#FF8F00') elif p < 0.10: COEF_COLORS.append('#FDD835') else: COEF_COLORS.append('#9E9E9E') fig3, ax3 = plt.subplots(figsize=(9, 5.5)) y_pos3 = np.arange(len(REG_LABELS)) ax3.barh(y_pos3, coefs3, color=COEF_COLORS, alpha=0.82, edgecolor='white', height=0.55) ax3.errorbar(coefs3, y_pos3, xerr=1.96 * ses3, fmt='none', color='#222', capsize=5, linewidth=1.5) ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0) ax3.set_yticks(y_pos3) ax3.set_yticklabels(REG_LABELS, fontsize=11) ax3.set_xlabel('標準化回帰係数（±1.96SE）', fontsize=11) ax3.set_title( f'消費支出の決定要因 — OLS重回帰係数（2022年, N={len(df_reg)}都道府県）\n' f'R²={ols_model.rsquared:.3f}（adj. R²={ols_model.rsquared_adj:.3f}）', fontsize=12, fontweight='bold' ) ax3.invert_yaxis() ax3.grid(axis='x', alpha=0.3) # p値ラベル for i, (c, p) in enumerate(zip(coefs3, pvals3)): sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.' offset = 0.01 ha = 'left' if c >= 0 else 'right' ax3.text(c + (offset if c >= 0 else -offset), i, f' {c:.3f} {sig}', va='center', ha=ha, fontsize=8.5) from matplotlib.patches import Patch ax3.legend(handles=[ Patch(color='#C62828', alpha=0.85, label='p<0.01'), Patch(color='#FF8F00', alpha=0.85, label='p<0.05'), Patch(color='#FDD835', alpha=0.85, label='p<0.10'), Patch(color='#9E9E9E', alpha=0.85, label='n.s.'), ], fontsize=9, loc='lower right') plt.tight_layout() fig3.savefig(os.path.join(FIG_DIR, '2022_H5_14_fig3_coef.png'), bbox_inches='tight', dpi=150) plt.close(fig3) print("図3保存: 2022_H5_14_fig3_coef.png") # ================================================================ # ■ 図4: 都市部 vs 地方 消費支出比較（箱ひげ図 + Mann-Whitney U検定） # ================================================================ print("\n■ 図4: 都市部 vs 地方 消費支出比較 箱ひげ図") fig4, axes4 = plt.subplots(1, 2, figsize=(12, 6)) # 左パネル: 箱ひげ図 ax4L = axes4[0] data_urban = urban.values data_rural = rural.values bp = ax4L.boxplot([data_urban, data_rural], labels=['都市部', '地方'], patch_artist=True, widths=0.5, medianprops=dict(color='white', linewidth=2.5), boxprops=dict(linewidth=1.5), whiskerprops=dict(linewidth=1.5), capprops=dict(linewidth=1.5), flierprops=dict(marker='o', markersize=5, alpha=0.6)) bp['boxes'][0].set_facecolor('#C62828') bp['boxes'][0].set_alpha(0.78) bp['boxes'][1].set_facecolor('#1565C0') bp['boxes'][1].set_alpha(0.78) # 個別データ点をジッタで表示 rng = np.random.default_rng(seed=42) # 表示用のジッタのみ（データ生成ではない） for i, (data, color) in enumerate(zip([data_urban, data_rural], ['#C62828', '#1565C0'])): jitter = rng.uniform(-0.12, 0.12, size=len(data)) ax4L.scatter(np.ones(len(data)) * (i + 1) + jitter, data, color=color, alpha=0.5, s=30, zorder=3) # 平均値をダイヤモンドでプロット for i, data in enumerate([data_urban, data_rural]): ax4L.scatter(i + 1, data.mean(), marker='D', color='white', s=60, zorder=5, edgecolors='#333', linewidths=1.5) # Mann-Whitney U検定の結果を注記 sig_txt = '***' if u_pval < 0.001 else '**' if u_pval < 0.01 else '*' if u_pval < 0.05 else 'n.s.' ax4L.text(1.5, max(data_urban.max(), data_rural.max()) * 1.01, f'Mann-Whitney U検定\nU={u_stat:.0f}, p={u_pval:.4f} {sig_txt}', ha='center', va='bottom', fontsize=9, color='#333', bbox=dict(boxstyle='round,pad=0.3', facecolor='#FFF9C4', edgecolor='#F9A825')) ax4L.set_ylabel('消費支出（万円/月）', fontsize=11) ax4L.set_title('都市部 vs 地方\n消費支出の分布（2022年）', fontsize=12, fontweight='bold') ax4L.grid(axis='y', alpha=0.3) # 右パネル: 年度別 格差推移（都市部平均 - 地方平均） ax4R = axes4[1] gap_ts = ( df_all.groupby(['年度', '地域区分'])['消費支出_万円'] .mean() .unstack('地域区分') ) if '都市部' in gap_ts.columns and '地方' in gap_ts.columns: gap_ts['格差'] = gap_ts['都市部'] - gap_ts['地方'] years = gap_ts.index.tolist() ax4R.fill_between(years, gap_ts['格差'], alpha=0.3, color='#C62828') ax4R.plot(years, gap_ts['格差'], marker='o', linewidth=2.2, markersize=5.5, color='#C62828', label='都市部 - 地方（格差）') ax4R.axhline(0, color='gray', linestyle='--', linewidth=1.0) ax4R.set_xlabel('年度', fontsize=11) ax4R.set_ylabel('消費支出の格差（万円/月）', fontsize=11) ax4R.set_title('都市農村間 消費支出格差の推移\n（都市部平均 − 地方平均）', fontsize=12, fontweight='bold') ax4R.set_xticks(years) ax4R.xaxis.set_tick_params(rotation=45) ax4R.legend(fontsize=9.5) ax4R.grid(True, alpha=0.3) plt.tight_layout() fig4.savefig(os.path.join(FIG_DIR, '2022_H5_14_fig4_boxplot.png'), bbox_inches='tight', dpi=150) plt.close(fig4) print("図4保存: 2022_H5_14_fig4_boxplot.png") # ================================================================ # ■ 完了メッセージ # ================================================================ print("\n" + "=" * 65) print("■ 全図生成完了（4枚）") print(f" fig1_timeseries.png : 消費支出の時系列推移（地方別平均）") print(f" fig2_scatter.png : 住宅地価 vs 消費支出 散布図（2022年）") print(f" fig3_coef.png : OLS回帰係数プロット") print(f" fig4_boxplot.png : 都市部 vs 地方 消費支出比較") print("=" * 65) print(f"\n OLS結果サマリ（2022年, N={len(df_reg)}）:") print(f" R²={ols_model.rsquared:.3f}, adj.R²={ols_model.rsquared_adj:.3f}") for label, coef, pval in zip(REG_LABELS, coefs3, pvals3): sig = '***' if pval < 0.001 else '**' if pval < 0.01 else '*' if pval < 0.05 else 'n.s.' print(f" {label}: β={coef:.3f} (p={pval:.4f}) {sig}") print(f"\n ジニ係数（消費支出, 2022年）: 全国={g_all:.4f}, 都市部={g_urban:.4f}, 地方={g_rural:.4f}")