""" 教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [高校生の部] ================================================================= 論文タイトル：健康寿命の地域差：生活習慣・医療環境の重回帰分析 受賞：審査員奨励賞（高校生の部） 【分析概要】 データ：SSDSE-B-2026.csv（都道府県別データ, 2012〜2023年度） 対象 ：47都道府県のパネルデータ 方法 ：相関分析、OLS重回帰分析、地域比較、時系列分析 目的変数（代理）： 保健医療費比率 = 保健医療費 / 消費支出 × 100（医療依存度の代理） 説明変数： ・高齢化率 = 65歳以上人口 / 総人口 × 100 ・食料費割合 = 食料費 / 消費支出 × 100 ・医療機関密度 = 一般病院数 / 総人口 × 10000 ・消費支出_log = log(消費支出) 【変数説明】 L322106 保健医療費（二人以上の世帯）: 月間の保健医療費支出（円） L3221 消費支出（二人以上の世帯） : 月間の消費支出合計（円） L322101 食料費（二人以上の世帯） : 月間の食料費支出（円） A1303 65歳以上人口 : 都道府県の高齢者人口（人） A1101 総人口 : 都道府県の総人口（人） I510120 一般病院数 : 都道府県の一般病院総数（施設） 【学習ポイント】 1. 代理変数（Proxy Variable）の設計：直接測定できない「健康寿命」を 「保健医療費比率」で代理する手法 2. 相関分析（Pearson）と散布図による変数間関係の把握 3. OLS重回帰分析の実装・係数の解釈（標準化 vs 非標準化） 4. 地域ブロック別分類と箱ひげ図による地域差の可視化 【データ出典】 SSDSE-B-2026.csv: 社会・人口統計体系（都道府県別データ） 統計センターより公表の実データ ================================================================= """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・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_2_shorei.py # ============================================================ import os import warnings 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 warnings.filterwarnings('ignore') # ────────────────────────────────────────────────────────────── # パス設定 # ────────────────────────────────────────────────────────────── FIG_DIR = 'html/figures' DATA_B = 'data/raw/SSDSE-B-2026.csv' os.makedirs(FIG_DIR, exist_ok=True) # ────────────────────────────────────────────────────────────── # matplotlib 設定 # ────────────────────────────────────────────────────────────── plt.rcParams['font.family'] = 'Hiragino Sans' plt.rcParams['axes.unicode_minus'] = False plt.rcParams['figure.dpi'] = 150 # ────────────────────────────────────────────────────────────── # 地域ブロック定義（8ブロック） # ────────────────────────────────────────────────────────────── REGION_MAP = { '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北', '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北', '福島県': '北海道・東北', '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東', '千葉県': '関東', '東京都': '関東', '神奈川県': '関東', '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部', '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部', '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿', '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿', '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国', '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国', '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国', '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄', '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄', '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄', } REGION_ORDER = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄'] REGION_COLORS = { '北海道・東北': '#1565C0', '関東': '#E65100', '中部': '#2E7D32', '近畿': '#6A1B9A', '中国・四国': '#00838F', '九州・沖縄': '#C62828', } # ================================================================ # ■ データ読み込み（SSDSE-B-2026: 都道府県別パネルデータ） # ================================================================ df_raw = pd.read_csv(DATA_B, encoding='cp932', header=1) df_b = df_raw[df_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].copy() # 使用列の数値変換 num_cols = [ '総人口', '65歳以上人口', '消費支出（二人以上の世帯）', '食料費（二人以上の世帯）', '保健医療費（二人以上の世帯）', '一般病院数', ] for c in num_cols: df_b[c] = pd.to_numeric(df_b[c], errors='coerce') df_b['年度'] = pd.to_numeric(df_b['年度'], errors='coerce') # ────────────────────────────────────────────────────────────── # 派生変数の計算 # ────────────────────────────────────────────────────────────── df_b['保健医療費比率'] = df_b['保健医療費（二人以上の世帯）'] / df_b['消費支出（二人以上の世帯）'] * 100 df_b['食料費割合'] = df_b['食料費（二人以上の世帯）'] / df_b['消費支出（二人以上の世帯）'] * 100 df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100 df_b['医療機関密度'] = df_b['一般病院数'] / df_b['総人口'] * 10000 df_b['消費支出_log'] = np.log(df_b['消費支出（二人以上の世帯）']) # 地域ブロック付与 df_b['地域ブロック'] = df_b['都道府県'].map(REGION_MAP) years = sorted(df_b['年度'].dropna().unique()) print("=" * 60) print("■ データ概要") print("=" * 60) print(f" 総観測数 : {len(df_b)} 行 ({df_b['都道府県'].nunique()} 都道府県 × {len(years)} 年度)") print(f" 期間 : {int(years[0])}〜{int(years[-1])} 年度") print() # ─── 2022年断面データ ─────────────────────────────────────────── df22 = df_b[df_b['年度'] == 2022].dropna( subset=['保健医療費比率', '高齢化率', '食料費割合', '医療機関密度', '消費支出_log'] ).copy() print("■ 2022年の基本統計") for col in ['保健医療費比率', '高齢化率', '食料費割合', '医療機関密度']: m = df22[col].mean() s = df22[col].std() print(f" {col:<14}: mean={m:.2f}, std={s:.2f}") print() # ─── 相関分析（2022年） ───────────────────────────────────────── y_col = '保健医療費比率' x_cols = ['高齢化率', '食料費割合', '医療機関密度', '消費支出_log'] print("■ 相関分析（2022年, n=47）") corr_results = {} for xc in x_cols: r, p = stats.pearsonr(df22[xc].dropna(), df22[y_col][df22[xc].notna()]) corr_results[xc] = {'r': r, 'p': p} sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.' print(f" {xc:<16}: r = {r:+.3f} p = {p:.4f} {sig}") print() # ─── OLS 重回帰（2022年） ────────────────────────────────────── X_ols = sm.add_constant(df22[x_cols]) model = sm.OLS(df22[y_col], X_ols).fit() print("■ OLS 重回帰結果（2022年）") print(f" R² = {model.rsquared:.4f}, 調整済みR² = {model.rsquared_adj:.4f}, n = {len(df22)}") print(f" F統計量 = {model.fvalue:.2f}, p(F) = {model.f_pvalue:.4f}") print() print(f" {'変数':<16} {'係数':>8} {'標準誤差':>8} {'t値':>7} {'p値':>8} 有意性") print(" " + "-" * 62) for nm, coef, se, t, p in zip( model.params.index, model.params.values, model.bse.values, model.tvalues.values, model.pvalues.values): sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.' print(f" {nm:<16} {coef:>8.4f} {se:>8.4f} {t:>7.3f} {p:>8.4f} {sig}") print() # ================================================================ # ■ 図の生成（4枚） # ================================================================ # ──────────────────────────────────────────────────────────────── # 図1: 保健医療費比率の時系列推移（地域ブロック別） # ──────────────────────────────────────────────────────────────── fig1, ax1 = plt.subplots(figsize=(12, 6)) for region in REGION_ORDER: df_reg = df_b[df_b['地域ブロック'] == region].copy() ts = df_reg.groupby('年度')['保健医療費比率'].mean().sort_index() ax1.plot(ts.index, ts.values, color=REGION_COLORS[region], linewidth=2.2, marker='o', markersize=4.5, label=region) ax1.set_xlabel('年度', fontsize=12) ax1.set_ylabel('保健医療費比率（保健医療費 / 消費支出 × 100）', fontsize=11) ax1.set_title('保健医療費比率の時系列推移（地域ブロック別平均）\n' 'SSDSE-B 2012〜2023年度, 二人以上世帯', fontsize=13, fontweight='bold') ax1.set_xticks(years) ax1.legend(fontsize=9.5, loc='upper left', ncol=2) ax1.grid(True, alpha=0.3, axis='y') ax1.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f'{x:.1f}%')) plt.tight_layout() fig1.savefig(os.path.join(FIG_DIR, '2022_H5_2_fig1_ts_health_ratio.png'), bbox_inches='tight', dpi=150) plt.close(fig1) print("図1保存: 2022_H5_2_fig1_ts_health_ratio.png") # ──────────────────────────────────────────────────────────────── # 図2: 高齢化率 vs 保健医療費比率 散布図（2022年） # ──────────────────────────────────────────────────────────────── fig2, ax2 = plt.subplots(figsize=(9, 7)) for region in REGION_ORDER: sub = df22[df22['地域ブロック'] == region] ax2.scatter(sub['高齢化率'], sub['保健医療費比率'], color=REGION_COLORS[region], s=65, alpha=0.82, label=region, zorder=3) # 回帰直線（全体） x_all = df22['高齢化率'].values y_all = df22['保健医療費比率'].values m_fit, b_fit, r_val, p_val, _ = stats.linregress(x_all, y_all) x_line = np.linspace(x_all.min(), x_all.max(), 100) ax2.plot(x_line, m_fit * x_line + b_fit, color='#333333', linewidth=1.8, linestyle='--', label=f'回帰直線 (r={r_val:+.3f}, p={p_val:.3f})', zorder=2) # 代表的な都道府県にラベル for _, row in df22.iterrows(): pref = row['都道府県'] if pref in ['秋田県', '東京都', '沖縄県', '高知県', '神奈川県', '愛知県']: ax2.annotate(pref, xy=(row['高齢化率'], row['保健医療費比率']), xytext=(5, 3), textcoords='offset points', fontsize=8.5, color='#333') ax2.set_xlabel('高齢化率（65歳以上人口 / 総人口 × 100）', fontsize=12) ax2.set_ylabel('保健医療費比率（保健医療費 / 消費支出 × 100）', fontsize=11) ax2.set_title('高齢化率 vs 保健医療費比率（2022年度, 47都道府県）', fontsize=13, fontweight='bold') ax2.legend(fontsize=9, loc='upper left') ax2.grid(True, alpha=0.3) plt.tight_layout() fig2.savefig(os.path.join(FIG_DIR, '2022_H5_2_fig2_scatter_aging.png'), bbox_inches='tight', dpi=150) plt.close(fig2) print("図2保存: 2022_H5_2_fig2_scatter_aging.png") # ──────────────────────────────────────────────────────────────── # 図3: OLS 回帰係数プロット（95%信頼区間付き） # ──────────────────────────────────────────────────────────────── fig3, ax3 = plt.subplots(figsize=(9, 5)) # 標準化係数を計算（比較のため） df22_std = df22[x_cols + [y_col]].copy() for c in df22_std.columns: df22_std[c] = (df22_std[c] - df22_std[c].mean()) / df22_std[c].std() X_std = sm.add_constant(df22_std[x_cols]) model_std = sm.OLS(df22_std[y_col], X_std).fit() vars_plot = x_cols coefs = model_std.params[vars_plot].values ci_low = model_std.conf_int().loc[vars_plot, 0].values ci_high = model_std.conf_int().loc[vars_plot, 1].values pvals = model_std.pvalues[vars_plot].values colors_bar = [] for p in pvals: if p < 0.001: colors_bar.append('#C62828') elif p < 0.05: colors_bar.append('#E65100') else: colors_bar.append('#9E9E9E') y_pos = np.arange(len(vars_plot)) ax3.barh(y_pos, coefs, color=colors_bar, alpha=0.82, height=0.55, zorder=3) ax3.errorbar(coefs, y_pos, xerr=[coefs - ci_low, ci_high - coefs], fmt='none', color='#333', capsize=5, linewidth=1.5, zorder=4) ax3.axvline(0, color='black', linewidth=1.2, linestyle='-') ax3.set_yticks(y_pos) ax3.set_yticklabels(vars_plot, fontsize=11) ax3.set_xlabel('標準化回帰係数（95%信頼区間）', fontsize=12) ax3.set_title('OLS重回帰：標準化回帰係数プロット（2022年度）\n' '目的変数：保健医療費比率', fontsize=13, fontweight='bold') # 凡例パッチ import matplotlib.patches as mpatches patches = [ mpatches.Patch(color='#C62828', label='p < 0.001 ***'), mpatches.Patch(color='#E65100', label='p < 0.05 *'), mpatches.Patch(color='#9E9E9E', label='n.s.'), ] ax3.legend(handles=patches, fontsize=9, loc='lower right') # p値アノテーション for i, (coef, p) in enumerate(zip(coefs, pvals)): sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.' ax3.text(ci_high[i] + 0.005, i, f' {sig}', va='center', fontsize=10) ax3.grid(True, alpha=0.3, axis='x') plt.tight_layout() fig3.savefig(os.path.join(FIG_DIR, '2022_H5_2_fig3_coef.png'), bbox_inches='tight', dpi=150) plt.close(fig3) print("図3保存: 2022_H5_2_fig3_coef.png") # ──────────────────────────────────────────────────────────────── # 図4: 地域ブロック別 医療機関密度の箱ひげ図（全年度データ） # ──────────────────────────────────────────────────────────────── fig4, ax4 = plt.subplots(figsize=(11, 6)) data_by_region = [] labels_region = [] colors_box = [] for region in REGION_ORDER: vals = df_b[df_b['地域ブロック'] == region]['医療機関密度'].dropna().values data_by_region.append(vals) labels_region.append(region) colors_box.append(REGION_COLORS[region]) bp = ax4.boxplot(data_by_region, labels=labels_region, patch_artist=True, notch=False, widths=0.5, medianprops=dict(color='white', linewidth=2.5), whiskerprops=dict(linewidth=1.5), capprops=dict(linewidth=1.5), flierprops=dict(marker='o', markersize=4, alpha=0.6)) for patch, color in zip(bp['boxes'], colors_box): patch.set_facecolor(color) patch.set_alpha(0.78) for i, (vals, region) in enumerate(zip(data_by_region, labels_region)): jitter = np.random.default_rng(42).uniform(-0.18, 0.18, size=len(vals)) ax4.scatter(np.full(len(vals), i + 1) + jitter, vals, color=REGION_COLORS[region], alpha=0.45, s=20, zorder=3) ax4.set_xlabel('地域ブロック', fontsize=12) ax4.set_ylabel('医療機関密度（一般病院数 / 総人口 × 10,000）', fontsize=11) ax4.set_title('地域ブロック別 医療機関密度の分布\n（SSDSE-B 2012〜2023年度, 47都道府県×12年分）', fontsize=13, fontweight='bold') ax4.grid(True, alpha=0.3, axis='y') plt.xticks(rotation=15, ha='right') plt.tight_layout() fig4.savefig(os.path.join(FIG_DIR, '2022_H5_2_fig4_boxplot_density.png'), bbox_inches='tight', dpi=150) plt.close(fig4) print("図4保存: 2022_H5_2_fig4_boxplot_density.png") print() print("=" * 60) print("■ 全図の生成完了") print(f" 保存先: {FIG_DIR}") print("=" * 60)