学校教育・部活動と地域の子どもの
体力・健康の関係分析

import os
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

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

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

# ── データ読み込み ──────────────────────────────────────────────
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)

# 2022年断面データ（47都道府県）
df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)
assert len(df) == 47, f"都道府県数が{len(df)}（47を想定）"

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

df['地域'] = df['都道府県'].map(region_map)
df['地域色'] = df['地域'].map(region_colors)

# ── 派生変数の計算 ─────────────────────────────────────────────
df['15歳未満人口比率']   = df['15歳未満人口'] / df['総人口'] * 100
df['65歳以上人口比率']   = df['65歳以上人口'] / df['総人口'] * 100
df['高校進学率']          = df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数'] * 100
df['小中学校合計']        = df['小学校数'] + df['中学校数']
df['小中児童生徒数']      = df['小学校児童数'] + df['中学校生徒数']
df['学校当たり児童生徒数'] = df['小中児童生徒数'] / df['小中学校合計']
df['保育所密度']          = df['保育所等数'] / df['総人口'] * 10000   # 人口1万人当たり
df['教育費']              = df['教育費（二人以上の世帯）']            # 円/月
df['消費支出']            = df['消費支出（二人以上の世帯）']          # 円/月
df['教育費比率']          = df['教育費'] / df['消費支出'] * 100       # %

# 都道府県名を短縮（ラベル用）
df['都道府県_短'] = df['都道府県'].str.replace('都', '都').str.replace('道', '道') \
                                   .str.replace('府', '府').str.replace('県', '県')

# ── 統計確認 ─────────────────────────────────────────────────
print("=== 2022年度 教育費（都道府県別） ===")
ranking = df.sort_values('教育費', ascending=False)[['都道府県', '教育費', '消費支出', '地域']].reset_index(drop=True)
ranking.index += 1
print(ranking.to_string())
print()

fig1, ax1 = plt.subplots(figsize=(10, 7))

for region, color in region_colors.items():
    mask = df['地域'] == region
    ax1.scatter(df.loc[mask, '消費支出'] / 1000,
                df.loc[mask, '教育費'] / 1000,
                c=color, label=region, s=60, alpha=0.85, zorder=3)

# 都道府県ラベル
for _, row in df.iterrows():
    ax1.annotate(row['都道府県'].replace('県', '').replace('都', '都').replace('府', '府').replace('道', '道'),
                 (row['消費支出'] / 1000, row['教育費'] / 1000),
                 fontsize=6.5, ha='center', va='bottom',
                 xytext=(0, 3), textcoords='offset points', color='#333333')

# 回帰線
x_vals = df['消費支出'] / 1000
y_vals = df['教育費'] / 1000
slope, intercept, r_val, p_val, _ = stats.linregress(x_vals, y_vals)
x_range = np.linspace(x_vals.min(), x_vals.max(), 100)
ax1.plot(x_range, slope * x_range + intercept,
         color='#555', linewidth=1.5, linestyle='--', zorder=2,
         label=f'回帰線 (r = {r_val:.3f}, p {"< 0.001" if p_val < 0.001 else f"= {p_val:.3f}"})')

ax1.set_xlabel('消費支出（千円/月）', fontsize=12)
ax1.set_ylabel('教育費（千円/月）', fontsize=12)
ax1.set_title('都道府県別 消費支出と教育費の関係（2022年度）', fontsize=14, fontweight='bold')
ax1.legend(fontsize=9, loc='upper left', framealpha=0.9)
ax1.grid(True, alpha=0.3)
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)

fig1.tight_layout()
path1 = os.path.join(FIG_DIR, '2020_H5_2_fig1.png')
fig1.savefig(path1, bbox_inches='tight')
plt.close(fig1)
print(f"図1 保存: {path1}")

df_sorted = df.sort_values('教育費', ascending=True).reset_index(drop=True)

fig2, ax2 = plt.subplots(figsize=(10, 13))

colors_bar = [region_colors[r] for r in df_sorted['地域']]
bars = ax2.barh(range(len(df_sorted)), df_sorted['教育費'] / 1000,
                color=colors_bar, alpha=0.85, edgecolor='white', linewidth=0.5)

ax2.set_yticks(range(len(df_sorted)))
ax2.set_yticklabels(df_sorted['都道府県'], fontsize=8.5)
ax2.set_xlabel('教育費（千円/月）', fontsize=11)
ax2.set_title('都道府県別 教育費ランキング（2022年度）', fontsize=14, fontweight='bold')

# 凡例
from matplotlib.patches import Patch
legend_elems = [Patch(facecolor=c, label=r) for r, c in region_colors.items()]
ax2.legend(handles=legend_elems, fontsize=9, loc='lower right', framealpha=0.9)

# 全国平均線
avg = df['教育費'].mean() / 1000
ax2.axvline(avg, color='#e05c5c', linewidth=1.5, linestyle='--', label=f'全国平均 {avg:.1f}千円')
ax2.text(avg + 0.2, len(df_sorted) - 1, f'全国平均\n{avg:.1f}千円', fontsize=8, color='#e05c5c', va='top')

ax2.grid(True, axis='x', alpha=0.3)
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)

fig2.tight_layout()
path2 = os.path.join(FIG_DIR, '2020_H5_2_fig2.png')
fig2.savefig(path2, bbox_inches='tight')
plt.close(fig2)
print(f"図2 保存: {path2}")

corr_vars = {
    '教育費': df['教育費'],
    '消費支出': df['消費支出'],
    '15歳未満\n人口比率': df['15歳未満人口比率'],
    '65歳以上\n人口比率': df['65歳以上人口比率'],
    '高校進学率': df['高校進学率'],
    '合計特殊\n出生率': df['合計特殊出生率'],
    '保育所密度': df['保育所密度'],
    '学校当たり\n児童生徒数': df['学校当たり児童生徒数'],
}
df_corr = pd.DataFrame(corr_vars)
corr_matrix = df_corr.corr()
n_vars = len(corr_matrix)

# p値行列
p_matrix = np.zeros((n_vars, n_vars))
cols_list = list(corr_vars.keys())
for i in range(n_vars):
    for j in range(n_vars):
        if i != j:
            _, p_val = stats.pearsonr(df_corr.iloc[:, i].dropna(),
                                      df_corr.iloc[:, j].dropna())
            p_matrix[i, j] = p_val

fig3, ax3 = plt.subplots(figsize=(9, 7))
cmap = plt.get_cmap('RdBu_r')
im = ax3.imshow(corr_matrix.values, cmap=cmap, vmin=-1, vmax=1, aspect='auto')

ax3.set_xticks(range(n_vars))
ax3.set_yticks(range(n_vars))
ax3.set_xticklabels(cols_list, fontsize=9)
ax3.set_yticklabels(cols_list, fontsize=9)

# セル内に相関係数を表示
for i in range(n_vars):
    for j in range(n_vars):
        val = corr_matrix.values[i, j]
        sig = '*' if (i != j and p_matrix[i, j] < 0.05) else ''
        text_color = 'white' if abs(val) > 0.6 else 'black'
        ax3.text(j, i, f'{val:.2f}{sig}', ha='center', va='center',
                 fontsize=8, color=text_color, fontweight='bold' if sig else 'normal')

plt.colorbar(im, ax=ax3, shrink=0.8, label='Pearson r')
ax3.set_title('教育・子ども関連変数の相関行列（2022年度・47都道府県）\n※ * は p<0.05', fontsize=12, fontweight='bold')

fig3.tight_layout()
path3 = os.path.join(FIG_DIR, '2020_H5_2_fig3.png')
fig3.savefig(path3, bbox_inches='tight')
plt.close(fig3)
print(f"図3 保存: {path3}")

pred_vars = {
    '消費支出': df['消費支出'],
    '15歳未満人口比率': df['15歳未満人口比率'],
    '65歳以上人口比率': df['65歳以上人口比率'],
    '高校進学率': df['高校進学率'],
    '合計特殊出生率': df['合計特殊出生率'],
    '保育所密度': df['保育所密度'],
}

# 標準化
y_raw = df['教育費'].values
X_raw = pd.DataFrame(pred_vars).values

y_std = (y_raw - y_raw.mean()) / y_raw.std(ddof=1)
X_std = (X_raw - X_raw.mean(axis=0)) / X_raw.std(axis=0, ddof=1)

X_std_sm = sm.add_constant(X_std)
model = sm.OLS(y_std, X_std_sm).fit()

print("\n=== 重回帰分析結果（標準化） ===")
print(model.summary())

# 係数（定数項除く）
coef_names = list(pred_vars.keys())
coefs = model.params[1:]
pvals = model.pvalues[1:]
conf_int_arr = model.conf_int()  # numpy array shape (n_params, 2)
conf_int_lo = conf_int_arr[1:, 0]
conf_int_hi = conf_int_arr[1:, 1]

colors_coef = ['#e05c5c' if p < 0.05 else '#aaaaaa' for p in pvals]

fig4, ax4 = plt.subplots(figsize=(8, 5))

y_pos = range(len(coef_names))
ax4.barh(y_pos, coefs, color=colors_coef, alpha=0.85, edgecolor='white', linewidth=0.5)

# 95%信頼区間
for i, (lo, hi, coef) in enumerate(zip(conf_int_lo, conf_int_hi, coefs)):
    ax4.plot([lo, hi], [i, i], color='black', linewidth=1.5, zorder=5)
    ax4.plot([lo, lo], [i - 0.1, i + 0.1], color='black', linewidth=1.5, zorder=5)
    ax4.plot([hi, hi], [i - 0.1, i + 0.1], color='black', linewidth=1.5, zorder=5)

ax4.axvline(0, color='black', linewidth=1, linestyle='-')
ax4.set_yticks(y_pos)
ax4.set_yticklabels(coef_names, fontsize=11)
ax4.set_xlabel('標準化偏回帰係数（β）', fontsize=11)
ax4.set_title(f'教育費の決定要因：標準化偏回帰係数\n（R² = {model.rsquared:.3f}, F = {model.fvalue:.2f}, p < 0.001）',
              fontsize=12, fontweight='bold')

# 凡例
from matplotlib.patches import Patch
legend_elems2 = [Patch(facecolor='#e05c5c', label='有意 (p < 0.05)'),
                 Patch(facecolor='#aaaaaa', label='非有意 (p ≥ 0.05)')]
ax4.legend(handles=legend_elems2, fontsize=9, loc='lower right')

# 有意な係数にラベル
for i, (coef, p) in enumerate(zip(coefs, pvals)):
    if p < 0.05:
        ax4.text(coef + (0.02 if coef >= 0 else -0.02), i,
                 f'β={coef:.3f}*', fontsize=8.5, va='center',
                 ha='left' if coef >= 0 else 'right', color='#c0392b')

ax4.grid(True, axis='x', alpha=0.3)
ax4.spines['top'].set_visible(False)
ax4.spines['right'].set_visible(False)

fig4.tight_layout()
path4 = os.path.join(FIG_DIR, '2020_H5_2_fig4.png')
fig4.savefig(path4, bbox_inches='tight')
plt.close(fig4)
print(f"図4 保存: {path4}")

# ── 統計サマリー出力 ─────────────────────────────────────────────
print("\n=== 主要統計サマリー ===")
print(f"教育費 全国平均  : {df['教育費'].mean():.0f} 円/月")
print(f"教育費 最大      : {df['教育費'].max():.0f} 円/月 ({df.loc[df['教育費'].idxmax(), '都道府県']})")
print(f"教育費 最小      : {df['教育費'].min():.0f} 円/月 ({df.loc[df['教育費'].idxmin(), '都道府県']})")
print(f"教育費 中央値    : {df['教育費'].median():.0f} 円/月")
print(f"消費支出 vs 教育費 相関: r = {r_val:.3f}")
print(f"重回帰 R²        : {model.rsquared:.3f}")
print(f"重回帰 自由度調整R²: {model.rsquared_adj:.3f}")
print()
print("偏回帰係数:")
for name, coef, p in zip(coef_names, coefs, pvals):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {name:20s}: β={coef:+.4f}  p={p:.4f} {sig}")

print("\n全図生成完了")