""" 教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [大学生・一般の部] ================================================================= 論文タイトル：欠測値を含むパネルデータを用いた健康寿命の要因分析 受賞：審査員奨励賞（大学生・一般の部） 【分析概要】 データ：SSDSE-B-2026.csv（都道府県別パネルデータ, 2012〜2023年度） 対象：全47都道府県 × 12年（2012〜2023） 健康代理変数：保健医療費割合 = 保健医療費 / 消費支出 × 100 (%) （注：健康寿命の直接データはSSDSE-Bに含まれないため、 　　　保健医療費割合を健康意識・健康投資の代理指標として使用） 分析の流れ ・欠測値可視化（教育的：実データにおける各変数の欠測状況の把握） ・保健医療費割合の時系列推移（全国平均 + 高低5県） ・パネルOLS 固定効果モデルによる保健医療費割合の決定要因推定 ・高齢化率 vs 保健医療費割合 散布図（2022年） 【被説明変数】 保健医療費割合(%) = 保健医療費（二人以上の世帯）/ 消費支出（二人以上の世帯）× 100 【説明変数】 高齢化率(%) = 65歳以上人口 / 総人口 × 100 食料費割合(%) = 食料費 / 消費支出 × 100 （エンゲル係数に相当） 病院密度 = 一般病院数 / 総人口 × 100000 （人口10万対病院数） 【推定手法】 linearmodels PanelOLS: entity_effects=True, time_effects=True, cov_type='clustered', cluster_entity=True 【データ出典】 SSDSE-B-2026.csv: 社会・人口統計体系（都道府県別データ） 統計センターより公表の実データ（2012〜2023年度） 【データサイエンス学習ポイント】 1. 欠測値の種類（MCAR/MAR/MNAR）と可視化の重要性 2. 平均代入 vs 多重代入：欠測処理の選択基準 3. 固定効果（FE）パネルモデルの within 変換 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_U5_8_shorei.py # ============================================================ import matplotlib matplotlib.use('Agg') import numpy as np import pandas as pd import matplotlib.pyplot as plt import matplotlib.patches as mpatches from matplotlib import rcParams from linearmodels import PanelOLS import warnings warnings.filterwarnings('ignore') import os # ────────────────────────────────────────────────────────────── # 共通設定 # ────────────────────────────────────────────────────────────── rcParams['font.family'] = 'Hiragino Sans' rcParams['axes.unicode_minus'] = False rcParams['figure.dpi'] = 150 _dir = os.path.dirname(os.path.abspath(__file__)) DATA_B = os.path.join(_dir, '..', 'data', 'raw', 'SSDSE-B-2026.csv') FIG_DIR = os.path.join(_dir, '..', 'html', 'figures') os.makedirs(FIG_DIR, exist_ok=True) # ================================================================ # ■ 実データ読み込み（SSDSE-B-2026: 都道府県別パネルデータ） # ================================================================ print("=" * 65) print("■ SSDSE-B-2026 データ読み込み") print("=" * 65) df = pd.read_csv(DATA_B, encoding='cp932', header=1) df = df[df['地域コード'].str.match(r'^R\d{5}$', na=False)].copy() print(f"都道府県数: {df['都道府県'].nunique()}") print(f"年度範囲 : {df['年度'].min()}〜{df['年度'].max()}") print(f"総行数 : {len(df)} (47都道府県 × 12年)") df = df.sort_values(['都道府県', '年度']).reset_index(drop=True) # ================================================================ # ■ 変数の構築 # ================================================================ # 健康代理変数（目的変数） df['保健医療費割合'] = df['保健医療費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）'] * 100 # 説明変数 df['高齢化率'] = df['65歳以上人口'] / df['総人口'] * 100 df['食料費割合'] = df['食料費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）'] * 100 df['病院密度'] = df['一般病院数'] / df['総人口'] * 100000 print("\n■ 主要変数 基本統計") stats_cols = ['保健医療費割合', '高齢化率', '食料費割合', '病院密度'] print(df[stats_cols].describe().round(3).to_string()) # ================================================================ # ■ 図1：欠測値棒グラフ（教育的：各変数の欠測状況） # ================================================================ print("\n" + "=" * 65) print("■ 図1：欠測値可視化") print("=" * 65) # 分析に関連する変数群の欠測カウント vis_cols = [ '保健医療費（二人以上の世帯）', '消費支出（二人以上の世帯）', '食料費（二人以上の世帯）', '住居費（二人以上の世帯）', '光熱・水道費（二人以上の世帯）', '教育費（二人以上の世帯）', '教養娯楽費（二人以上の世帯）', '総人口', '65歳以上人口', '15歳未満人口', '一般病院数', '一般診療所数', '保育所等数', '合計特殊出生率', '標準価格（平均価格）（住宅地）', '標準価格（平均価格）（商業地）', '大学数', '保育所等利用待機児童数', ] miss_counts = df[vis_cols].isnull().sum() n_total = len(df) # 変数名ラベル（短縮版） short_labels = { '保健医療費（二人以上の世帯）' : '保健医療費', '消費支出（二人以上の世帯）' : '消費支出', '食料費（二人以上の世帯）' : '食料費', '住居費（二人以上の世帯）' : '住居費', '光熱・水道費（二人以上の世帯）' : '光熱・水道費', '教育費（二人以上の世帯）' : '教育費', '教養娯楽費（二人以上の世帯）' : '教養娯楽費', '総人口' : '総人口', '65歳以上人口' : '65歳以上人口', '15歳未満人口' : '15歳未満人口', '一般病院数' : '一般病院数', '一般診療所数' : '一般診療所数', '保育所等数' : '保育所等数', '合計特殊出生率' : '合計特殊出生率', '標準価格（平均価格）（住宅地）' : '住宅地価格', '標準価格（平均価格）（商業地）' : '商業地価格', '大学数' : '大学数', '保育所等利用待機児童数' : '待機児童数', } labels_short = [short_labels[c] for c in vis_cols] miss_pcts = miss_counts.values / n_total * 100 # 色：欠測ゼロ＝青緑（完全）、欠測あり＝オレンジ bar_colors_miss = ['#78909C' if v == 0 else '#E65100' for v in miss_counts.values] fig1, ax1 = plt.subplots(figsize=(11, 5.5)) bars1 = ax1.bar(np.arange(len(vis_cols)), miss_pcts, color=bar_colors_miss, alpha=0.85, edgecolor='white', width=0.65) ax1.set_xticks(np.arange(len(vis_cols))) ax1.set_xticklabels(labels_short, rotation=45, ha='right', fontsize=9) ax1.set_ylabel('欠測率（%）', fontsize=11) ax1.set_title('SSDSE-B-2026（都道府県別パネルデータ）各変数の欠測状況\n' f'全{n_total}観測（47都道府県 × 12年）', fontsize=12, fontweight='bold') ax1.set_ylim(0, max(miss_pcts.max() + 5, 15)) ax1.axhline(5, color='#C62828', linestyle='--', linewidth=1.0, alpha=0.7, label='欠測率 5% 目安') ax1.grid(axis='y', alpha=0.3) # 欠測数・割合ラベル for i, (v, pct) in enumerate(zip(miss_counts.values, miss_pcts)): label = f'{v}件\n({pct:.1f}%)' if v > 0 else '0件\n(完全)' ax1.text(i, max(pct + 0.3, 0.5), label, ha='center', va='bottom', fontsize=7.5, color='#333') # 凡例 patch_ok = mpatches.Patch(color='#78909C', alpha=0.85, label='欠測なし（完全観測）') patch_ng = mpatches.Patch(color='#E65100', alpha=0.85, label='欠測あり') ax1.legend(handles=[patch_ok, patch_ng], fontsize=9, loc='upper right') # 教育的注釈 ax1.text(0.01, 0.97, '★ 欠測パターンの確認は前処理の第一歩（MCAR/MAR/MNARの判断に必要）', transform=ax1.transAxes, fontsize=8.5, color='#1565C0', va='top', style='italic') plt.tight_layout() fig1_path = os.path.join(FIG_DIR, '2022_U5_8_fig1_missing.png') fig1.savefig(fig1_path, bbox_inches='tight', dpi=150) plt.close(fig1) print(f"図1保存: {fig1_path}") # ================================================================ # ■ 図2：保健医療費割合 時系列（全国平均 + 高低5県） # ================================================================ print("\n■ 図2：保健医療費割合 時系列") nat_avg_ts = df.groupby('年度')['保健医療費割合'].mean() # 2022年の値で上位5・下位5県を選定 df_2022 = df[df['年度'] == 2022].copy() top5_prefs = df_2022.nlargest(5, '保健医療費割合')['都道府県'].tolist() bot5_prefs = df_2022.nsmallest(5, '保健医療費割合')['都道府県'].tolist() # 色マップ cmap_top = plt.cm.Reds(np.linspace(0.5, 0.9, 5)) cmap_bot = plt.cm.Blues(np.linspace(0.5, 0.9, 5)) color_map = {} for i, p in enumerate(top5_prefs): color_map[p] = cmap_top[i] for i, p in enumerate(bot5_prefs): color_map[p] = cmap_bot[i] highlight_prefs = top5_prefs + bot5_prefs fig2, ax2 = plt.subplots(figsize=(11, 5.5)) # 全47都道府県（グレー背景） for pref in df['都道府県'].unique(): if pref not in highlight_prefs: tmp = df[df['都道府県'] == pref].sort_values('年度') ax2.plot(tmp['年度'], tmp['保健医療費割合'], color='#BDBDBD', linewidth=0.6, alpha=0.40, zorder=1) # ハイライト：高位5県（赤系） for i, pref in enumerate(top5_prefs): tmp = df[df['都道府県'] == pref].sort_values('年度') val = df_2022[df_2022['都道府県'] == pref]['保健医療費割合'].values[0] ax2.plot(tmp['年度'], tmp['保健医療費割合'], color=color_map[pref], linewidth=1.8, marker='o', markersize=3.5, label=f'{pref}（{val:.1f}%）', zorder=4) # ハイライト：低位5県（青系） for i, pref in enumerate(bot5_prefs): tmp = df[df['都道府県'] == pref].sort_values('年度') val = df_2022[df_2022['都道府県'] == pref]['保健医療費割合'].values[0] ax2.plot(tmp['年度'], tmp['保健医療費割合'], color=color_map[pref], linewidth=1.8, marker='s', markersize=3.5, label=f'{pref}（{val:.1f}%）', zorder=4) # 全国平均（太線・黒破線） ax2.plot(nat_avg_ts.index, nat_avg_ts.values, color='#212121', linewidth=2.5, linestyle='--', label=f'全国平均（2022: {nat_avg_ts[2022]:.1f}%）', zorder=5) ax2.set_xlabel('年度', fontsize=12) ax2.set_ylabel('保健医療費割合（%）', fontsize=12) ax2.set_title('保健医療費割合の時系列推移\n全47都道府県（2012〜2023年度）— 上位5県（赤）・下位5県（青）・全国平均（黒破線）', fontsize=12, fontweight='bold') ax2.set_xticks(sorted(df['年度'].unique())) ax2.tick_params(axis='x', rotation=45) ax2.grid(True, alpha=0.3) ax2.legend(fontsize=8, loc='upper left', ncol=2, framealpha=0.9) # COVID注釈 ymin, ymax = ax2.get_ylim() ax2.axvspan(2019.5, 2021.5, color='#FFECB3', alpha=0.35, zorder=0) ax2.text(2020.5, ymin + (ymax - ymin) * 0.04, 'COVID-19', ha='center', fontsize=8, color='#795548', style='italic') plt.tight_layout() fig2_path = os.path.join(FIG_DIR, '2022_U5_8_fig2_ts.png') fig2.savefig(fig2_path, bbox_inches='tight', dpi=150) plt.close(fig2) print(f"図2保存: {fig2_path}") # ================================================================ # ■ パネルOLS 固定効果モデル推定 # ================================================================ print("\n" + "=" * 65) print("■ パネルOLS 固定効果モデル推定（双方向固定効果）") print("=" * 65) YVAR = '保健医療費割合' XVARS = ['高齢化率', '食料費割合', '病院密度'] df_panel = df.copy().set_index(['都道府県', '年度']) panel_data = df_panel[[YVAR] + XVARS].dropna() print(f"推定サンプル: {len(panel_data)} obs（{panel_data.index.get_level_values(0).nunique()}都道府県）") mod = PanelOLS( dependent=panel_data[YVAR], exog=panel_data[XVARS], entity_effects=True, time_effects=True, ) res = mod.fit(cov_type='clustered', cluster_entity=True) print(res.summary) fe_params = res.params fe_stderr = res.std_errors fe_pvalues = res.pvalues fe_rsq = res.rsquared print("\n■ 推定結果まとめ") print(f" Within R²: {fe_rsq:.4f}") var_labels = { '高齢化率' : '高齢化率（%）', '食料費割合' : '食料費割合（%）', '病院密度' : '病院密度（人口10万対）', } for var in XVARS: sig = ('***' if fe_pvalues[var] < 0.001 else '**' if fe_pvalues[var] < 0.01 else '*' if fe_pvalues[var] < 0.05 else 'n.s.') print(f" {var:<12} β={fe_params[var]:+.4f} SE={fe_stderr[var]:.4f} " f"p={fe_pvalues[var]:.4f} {sig}") # ================================================================ # ■ 図3：FE係数プロット # ================================================================ print("\n■ 図3：FE係数プロット") coefs = [fe_params[v] for v in XVARS] ses = [fe_stderr[v] for v in XVARS] pvals = [fe_pvalues[v] for v in XVARS] labels = [var_labels[v] for v in XVARS] bar_colors_fe = [] for p in pvals: if p < 0.01: bar_colors_fe.append('#C62828') elif p < 0.05: bar_colors_fe.append('#FF8F00') else: bar_colors_fe.append('#9E9E9E') fig3, ax3 = plt.subplots(figsize=(9, 4.5)) ypos = np.arange(len(XVARS)) ax3.barh(ypos, coefs, color=bar_colors_fe, alpha=0.82, edgecolor='white', height=0.5) ax3.errorbar(coefs, ypos, xerr=1.96 * np.array(ses), fmt='none', color='#333', capsize=5, linewidth=1.8) ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0) ax3.set_yticks(ypos) ax3.set_yticklabels(labels, fontsize=12) ax3.set_xlabel('固定効果モデル 推定係数（±1.96 SE）', fontsize=11) ax3.set_title(f'保健医療費割合の決定要因（パネルOLS 双方向固定効果モデル）\n' f'Within R²={fe_rsq:.3f}、N=47都道府県×12年（クラスタリング標準誤差）', fontsize=12, fontweight='bold') ax3.invert_yaxis() ax3.grid(axis='x', alpha=0.3) # p値注釈 x_range = max(coefs) - min(coefs) if len(coefs) > 1 else abs(coefs[0]) * 2 margin = max(x_range * 0.04, 0.001) for i, (c, p, se) in enumerate(zip(coefs, pvals, ses)): sig = ('***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.') x_ann = c + 1.96 * se + margin if c >= 0 else c - 1.96 * se - margin ha = 'left' if c >= 0 else 'right' ax3.text(x_ann, i, sig, va='center', ha=ha, fontsize=10, color=bar_colors_fe[i], fontweight='bold') legend_patches = [ mpatches.Patch(color='#C62828', alpha=0.82, label='p < 0.01'), mpatches.Patch(color='#FF8F00', alpha=0.82, label='p < 0.05'), mpatches.Patch(color='#9E9E9E', alpha=0.82, label='n.s. (p ≥ 0.05)'), ] ax3.legend(handles=legend_patches, fontsize=9, loc='lower right') plt.tight_layout() fig3_path = os.path.join(FIG_DIR, '2022_U5_8_fig3_fe_coef.png') fig3.savefig(fig3_path, bbox_inches='tight', dpi=150) plt.close(fig3) print(f"図3保存: {fig3_path}") # ================================================================ # ■ 図4：高齢化率 vs 保健医療費割合 散布図（2022年） # ================================================================ print("\n■ 図4：高齢化率 vs 保健医療費割合 散布図（2022年）") df_2022 = df[df['年度'] == 2022].copy() # 地方区分カラー region_map = { '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北', '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北', '福島県': '北海道・東北', '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東', '千葉県': '関東', '東京都': '関東', '神奈川県': '関東', '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部', '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部', '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿', '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿', '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国', '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国', '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国', '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄', '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄', '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄', } region_colors_map = { '北海道・東北': '#1565C0', '関東': '#C62828', '中部': '#2E7D32', '近畿': '#E65100', '中国・四国': '#6A1B9A', '九州・沖縄': '#00838F', } df_2022['地方'] = df_2022['都道府県'].map(region_map) fig4, ax4 = plt.subplots(figsize=(10, 7)) # 地方別散布 for region, grp in df_2022.groupby('地方'): ax4.scatter(grp['高齢化率'], grp['保健医療費割合'], color=region_colors_map.get(region, '#999'), s=70, alpha=0.75, zorder=3, label=region, edgecolors='white', linewidths=0.8) # 都道府県ラベル（全47） for _, row in df_2022.iterrows(): ax4.annotate( row['都道府県'].replace('県', '').replace('都', '').replace('道', '').replace('府', ''), xy=(row['高齢化率'], row['保健医療費割合']), xytext=(3, 3), textcoords='offset points', fontsize=6.5, color='#444', zorder=5, ) # 回帰直線 x_v = df_2022['高齢化率'].values y_v = df_2022['保健医療費割合'].values z = np.polyfit(x_v, y_v, 1) xs = np.linspace(x_v.min(), x_v.max(), 200) ax4.plot(xs, np.poly1d(z)(xs), color='#333', linewidth=1.8, linestyle='--', alpha=0.8, zorder=4, label='回帰直線') from scipy import stats as sp_stats r_val, p_val = sp_stats.pearsonr(x_v, y_v) ax4.set_xlabel('高齢化率（%）', fontsize=12) ax4.set_ylabel('保健医療費割合（%）', fontsize=12) ax4.set_title(f'高齢化率 vs 保健医療費割合（2022年度・47都道府県）\n' f'ピアソン相関: r = {r_val:.3f}（p = {p_val:.4f}）', fontsize=12, fontweight='bold') ax4.legend(fontsize=9, loc='upper right', framealpha=0.9) ax4.grid(True, alpha=0.3) plt.tight_layout() fig4_path = os.path.join(FIG_DIR, '2022_U5_8_fig4_scatter.png') fig4.savefig(fig4_path, bbox_inches='tight', dpi=150) plt.close(fig4) print(f"図4保存: {fig4_path}") # ================================================================ # ■ 完了サマリ # ================================================================ print("\n" + "=" * 65) print("■ 全図の生成完了（4枚）") print("=" * 65) print(f" 2022_U5_8_fig1_missing.png : 欠測値棒グラフ（変数別欠測率）") print(f" 2022_U5_8_fig2_ts.png : 保健医療費割合 時系列（全国 + 高低5県）") print(f" 2022_U5_8_fig3_fe_coef.png : FE係数プロット（保健医療費割合の決定要因）") print(f" 2022_U5_8_fig4_scatter.png : 高齢化率 vs 保健医療費割合 散布図（2022年）") print(f"\n パネルOLS結果（双方向固定効果、クラスタリング標準誤差）：") for var in XVARS: sig = ('***' if fe_pvalues[var] < 0.001 else '**' if fe_pvalues[var] < 0.01 else '*' if fe_pvalues[var] < 0.05 else 'n.s.') print(f" {var_labels[var]:<22} β={fe_params[var]:+.4f} " f"p={fe_pvalues[var]:.4f} {sig}") print(f"\n Within R² = {fe_rsq:.4f}")