子どもの体力・運動能力と
地域環境・家庭環境の関係分析

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年断面
df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)
assert len(df) == 47, f"47都道府県が揃っていません: {len(df)}行"

# ── 変数作成 ─────────────────────────────────────────
# 子ども関連指標（体力・運動環境の代理変数）
df['学校当たり児童数']     = df['小学校児童数'] / df['小学校数']       # 学校規模（大→都市型・活発）
df['15歳未満人口比率']     = df['15歳未満人口'] / df['総人口'] * 100   # 子ども人口の豊かさ
df['教育費比率']           = df['教育費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）'] * 100
df['教養娯楽費比率']       = df['教養娯楽費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）'] * 100

# 地域・家庭環境指標
df['保育所密度']           = df['保育所等数'] / df['総人口'] * 10000   # 1万人当たり保育所数
df['保育所定員率']         = df['保育所等定員数'] / df['15歳未満人口'] * 100
df['出生率']               = df['合計特殊出生率']
df['消費支出_万円']        = df['消費支出（二人以上の世帯）'] / 10000   # 万円単位

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

print("=" * 60)
print("SSDSE-B 2022年断面 47都道府県 基本統計")
print("=" * 60)

target = '学校当たり児童数'
features = ['15歳未満人口比率', '教育費比率', '教養娯楽費比率',
            '保育所密度', '保育所定員率', '出生率', '消費支出_万円']

analysis_cols = [target] + features
print(df[analysis_cols].describe().round(2))

# ── 相関分析 ──────────────────────────────────────────
print("\n── Pearson相関（目的変数: 学校当たり児童数）──")
for feat in features:
    r, p = stats.pearsonr(df[target], df[feat])
    sig = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else ''))
    print(f"  {feat:20s}: r={r:+.3f}, p={p:.4f} {sig}")

# ── 重回帰分析 ───────────────────────────────────────
print("\n── OLS重回帰（目的変数: 学校当たり児童数）──")
X = sm.add_constant(df[features])
model = sm.OLS(df[target], X).fit()
print(model.summary())

# 標準化偏回帰係数
df_z = df[analysis_cols].apply(lambda s: (s - s.mean()) / s.std())
Xz = sm.add_constant(df_z[features])
model_z = sm.OLS(df_z[target], Xz).fit()
beta = model_z.params.drop('const')
pval = model_z.pvalues.drop('const')

print("\n── 標準化偏回帰係数 ──")
for feat in features:
    print(f"  {feat:20s}: β={beta[feat]:+.3f}, p={pval[feat]:.4f}")

print(f"\nR²={model.rsquared:.3f}, 調整済みR²={model.rsquared_adj:.3f}, F={model.fvalue:.2f}, p={model.f_pvalue:.4f}")

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

for region, grp in df.groupby('地域区分'):
    col = region_colors[region]
    ax1.scatter(grp['15歳未満人口比率'], grp['学校当たり児童数'],
                color=col, label=region, s=70, alpha=0.85, zorder=3, edgecolors='white', linewidths=0.5)
    for _, row in grp.iterrows():
        ax1.annotate(row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', ''),
                     (row['15歳未満人口比率'], row['学校当たり児童数']),
                     fontsize=6.5, ha='left', va='bottom',
                     xytext=(2, 2), textcoords='offset points', alpha=0.8)

# 回帰線
x_vals = df['15歳未満人口比率'].values
y_vals = df['学校当たり児童数'].values
slope, intercept, r_val, p_val, _ = stats.linregress(x_vals, y_vals)
x_line = np.linspace(x_vals.min(), x_vals.max(), 100)
ax1.plot(x_line, intercept + slope * x_line, 'k--', linewidth=1.5,
         label=f'回帰線 (r={r_val:.2f}, p={p_val:.3f})', zorder=2)

ax1.set_xlabel('15歳未満人口比率（%）', fontsize=12)
ax1.set_ylabel('学校当たり児童数（人）', fontsize=12)
ax1.set_title('散布図：子ども人口比率と学校規模の関係\n（2022年・47都道府県）', fontsize=13, fontweight='bold')
ax1.legend(loc='upper left', fontsize=9, framealpha=0.9)
ax1.grid(True, alpha=0.3, linestyle='--')
ax1.spines[['top', 'right']].set_visible(False)

fig1.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2020_U5_5_fig1.png'), dpi=150, bbox_inches='tight')
plt.close(fig1)
print("\n図1 保存完了: 2020_U5_5_fig1.png")

heatmap_cols = ['学校当たり児童数', '15歳未満人口比率', '教育費比率',
                '教養娯楽費比率', '保育所密度', '保育所定員率', '出生率', '消費支出_万円']
heatmap_labels = ['学校当たり\n児童数', '15歳未満\n人口比率', '教育費\n比率',
                  '教養娯楽費\n比率', '保育所\n密度', '保育所\n定員率', '出生率', '消費支出\n(万円)']

corr_mat = df[heatmap_cols].corr()

fig2, ax2 = plt.subplots(figsize=(9, 7))
n = len(heatmap_cols)

import matplotlib.colors as mcolors
cmap = plt.cm.RdBu_r
norm = mcolors.Normalize(vmin=-1, vmax=1)

im = ax2.imshow(corr_mat.values, cmap=cmap, norm=norm, aspect='auto')
plt.colorbar(im, ax=ax2, fraction=0.046, pad=0.04, label='Pearson r')

ax2.set_xticks(range(n))
ax2.set_yticks(range(n))
ax2.set_xticklabels(heatmap_labels, fontsize=9)
ax2.set_yticklabels(heatmap_labels, fontsize=9)

for i in range(n):
    for j in range(n):
        val = corr_mat.values[i, j]
        r, p = stats.pearsonr(df[heatmap_cols[i]], df[heatmap_cols[j]])
        star = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else ''))
        text_color = 'white' if abs(val) > 0.6 else 'black'
        ax2.text(j, i, f'{val:.2f}{star}', ha='center', va='center',
                 fontsize=8.5, color=text_color, fontweight='bold' if p < 0.05 else 'normal')

ax2.set_title('相関ヒートマップ：子ども・地域環境関連変数\n（2022年・47都道府県、*p<0.05, **p<0.01, ***p<0.001）',
              fontsize=12, fontweight='bold')
fig2.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2020_U5_5_fig2.png'), dpi=150, bbox_inches='tight')
plt.close(fig2)
print("図2 保存完了: 2020_U5_5_fig2.png")

beta_vals = beta.values
pval_vals = pval.values
feat_labels_short = ['15歳未満\n人口比率', '教育費\n比率', '教養娯楽費\n比率',
                     '保育所\n密度', '保育所\n定員率', '出生率', '消費支出\n(万円)']

# 信頼区間（95%）
ci_lower = model_z.conf_int(alpha=0.05).loc[features, 0].values
ci_upper = model_z.conf_int(alpha=0.05).loc[features, 1].values

# 大きい順に並び替え
sort_idx = np.argsort(beta_vals)
beta_sorted = beta_vals[sort_idx]
pval_sorted = pval_vals[sort_idx]
labels_sorted = [feat_labels_short[i] for i in sort_idx]
ci_lo_sorted = ci_lower[sort_idx]
ci_hi_sorted = ci_upper[sort_idx]

colors_bar = ['#e05c5c' if v > 0 else '#4e9af1' for v in beta_sorted]
edge_colors = ['#c62828' if p < 0.05 else '#999' for p in pval_sorted]

fig3, ax3 = plt.subplots(figsize=(9, 6))
bars = ax3.barh(range(len(beta_sorted)), beta_sorted, color=colors_bar, alpha=0.75,
                edgecolor=edge_colors, linewidth=1.5)

# エラーバー
xerr_lo = beta_sorted - ci_lo_sorted
xerr_hi = ci_hi_sorted - beta_sorted
ax3.errorbar(beta_sorted, range(len(beta_sorted)),
             xerr=[xerr_lo, xerr_hi], fmt='none', color='black', capsize=4, linewidth=1)

# p値アノテーション
for i, (b, p) in enumerate(zip(beta_sorted, pval_sorted)):
    star = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else 'n.s.'))
    x_pos = b + (xerr_hi[i] + 0.02) if b >= 0 else b - xerr_lo[i] - 0.02
    ha = 'left' if b >= 0 else 'right'
    ax3.text(x_pos, i, star, ha=ha, va='center', fontsize=10,
             color='#c62828' if p < 0.05 else '#999')

ax3.axvline(0, color='black', linewidth=0.8, linestyle='-')
ax3.set_yticks(range(len(labels_sorted)))
ax3.set_yticklabels(labels_sorted, fontsize=10)
ax3.set_xlabel('標準化偏回帰係数 (β)', fontsize=12)
ax3.set_title(f'標準化偏回帰係数：学校当たり児童数の決定要因\n（OLS, R²={model.rsquared:.3f}, 調整済みR²={model.rsquared_adj:.3f}）',
              fontsize=12, fontweight='bold')
ax3.grid(True, axis='x', alpha=0.3, linestyle='--')
ax3.spines[['top', 'right']].set_visible(False)

from matplotlib.patches import Patch
legend_elements = [Patch(facecolor='#e05c5c', alpha=0.75, label='正の効果（学校規模を増大）'),
                   Patch(facecolor='#4e9af1', alpha=0.75, label='負の効果（学校規模を縮小）')]
ax3.legend(handles=legend_elements, loc='lower right', fontsize=9)

fig3.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2020_U5_5_fig3.png'), dpi=150, bbox_inches='tight')
plt.close(fig3)
print("図3 保存完了: 2020_U5_5_fig3.png")

df_rank = df.sort_values('15歳未満人口比率', ascending=True).reset_index(drop=True)

fig4, ax4 = plt.subplots(figsize=(10, 12))
bar_colors = [region_colors[r] for r in df_rank['地域区分']]
bars4 = ax4.barh(range(len(df_rank)), df_rank['15歳未満人口比率'],
                 color=bar_colors, alpha=0.85, edgecolor='white', linewidth=0.5)

# 全国平均線
nat_avg = df['15歳未満人口比率'].mean()
ax4.axvline(nat_avg, color='black', linewidth=1.5, linestyle='--',
            label=f'全国平均: {nat_avg:.2f}%')

# 都道府県名ラベル
ax4.set_yticks(range(len(df_rank)))
ax4.set_yticklabels(df_rank['都道府県'], fontsize=8)

# 値ラベル
for i, val in enumerate(df_rank['15歳未満人口比率']):
    ax4.text(val + 0.05, i, f'{val:.1f}%', va='center', fontsize=7.5, color='#333')

ax4.set_xlabel('15歳未満人口比率（%）', fontsize=12)
ax4.set_title('都道府県別 子ども人口比率ランキング\n（2022年・15歳未満人口 / 総人口×100）',
              fontsize=13, fontweight='bold')
ax4.spines[['top', 'right']].set_visible(False)
ax4.grid(True, axis='x', alpha=0.3, linestyle='--')

# 凡例
legend_handles = [Patch(facecolor=col, label=reg, alpha=0.85)
                  for reg, col in region_colors.items()]
legend_handles.append(plt.Line2D([0], [0], color='black', linewidth=1.5,
                                  linestyle='--', label=f'全国平均: {nat_avg:.2f}%'))
ax4.legend(handles=legend_handles, loc='lower right', fontsize=9, framealpha=0.9)

fig4.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2020_U5_5_fig4.png'), dpi=150, bbox_inches='tight')
plt.close(fig4)
print("図4 保存完了: 2020_U5_5_fig4.png")

# ── 統計サマリー出力（HTML生成用） ────────────────────
print("\n" + "=" * 60)
print("HTML生成用 統計サマリー")
print("=" * 60)
top3_child = df_rank.tail(3)['都道府県'].tolist()
bot3_child = df_rank.head(3)['都道府県'].tolist()
print(f"15歳未満人口比率 上位3位: {top3_child}")
print(f"15歳未満人口比率 下位3位: {bot3_child}")
print(f"全国平均 15歳未満人口比率: {nat_avg:.2f}%")
print(f"全国平均 学校当たり児童数: {df['学校当たり児童数'].mean():.1f}人")
print(f"全国平均 保育所密度: {df['保育所密度'].mean():.2f}（1万人当たり）")
print(f"全国平均 合計特殊出生率: {df['出生率'].mean():.2f}")
print(f"学校当たり児童数 vs 15歳未満人口比率: r={stats.pearsonr(df['学校当たり児童数'], df['15歳未満人口比率'])[0]:.3f}")
print(f"R²={model.rsquared:.3f}, 調整済みR²={model.rsquared_adj:.3f}")
print(f"有意な係数: {', '.join([f for f, p in zip(features, model.pvalues[1:]) if p < 0.05])}")

# 地域別平均
print("\n── 地域別 15歳未満人口比率 ──")
reg_stats = df.groupby('地域区分')['15歳未満人口比率'].agg(['mean', 'min', 'max']).round(2)
print(reg_stats)

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