市区町村別に見た
少子化の要因分析

df = pd.read_csv(DATA_B, encoding='cp932', header=0, skiprows=[1])
df.columns = [c.strip() for c in df.columns]
YEAR = 'SSDSE-B-2026'
d = df[df[YEAR] == 2023].copy()          # 最新年（2023年）の47都道府県を抽出
for c in ['A1101', 'A1303', 'A4101', 'A4103', 'A9101',
          'A510102', 'A510202', 'E6302']:
    d[c] = pd.to_numeric(d[c], errors='coerce')

# 目的変数・指標を作成（原論文の指標を都道府県粒度で近似）
d['TFR']       = d['A4103']                          # 合計特殊出生率
d['婚姻率']    = d['A9101'] / d['A1101'] * 1000       # 人口千人当たり婚姻件数
d['高齢化率']  = d['A1303'] / d['A1101'] * 100        # 65歳以上人口割合
d['女性転入出比'] = d['A510102'] / d['A510202']       # 女性転入/女性転出（全年齢の近似）
d['大学収容力'] = d['E6302'] / d['A1101'] * 1000      # 人口千人当たり大学学生数

print('SSDSE-B 読み込み: N =', len(d), '都道府県  （抽出年: 2023年）')
print(f'TFR: 平均 {d["TFR"].mean():.2f}  範囲 {d["TFR"].min():.2f}〜{d["TFR"].max():.2f}')
rm, pm = stats.pearsonr(d['婚姻率'], d['TFR'])
print(f'参考）婚姻率×TFR の相関 r={rm:+.3f}（p={pm:.3f}）'

r1 = d.sort_values('TFR', ascending=True)
print('最低5県:')
for _, r in r1.head(5).iterrows():
    print(f'  {r["Prefecture"]:6s} {r["TFR"]:.2f}')
print('最高5県:')
for _, r in r1.tail(5).iloc[::-1].iterrows():
    print(f'  {r["Prefecture"]:6s} {r["TFR"]:.2f}')

print('=== 図2: TFRと女性転入出比の散布図（実再現・近似）===')
m2 = d[['女性転入出比', 'TFR', 'Prefecture']].dropna()
r2, p2 = stats.pearsonr(m2['女性転入出比'], m2['TFR'])
print(f'女性転入出比 × TFR の相関 r={r2:+.3f}（p={p2:.3f}, N={len(m2)}）')
print('女性転入出比が高い（流入超過）5県:')
for _, r in m2.nlargest(5, '女性転入出比').iterrows():
    print(f'  {r["Prefecture"]:6s} 転入出比 {r["女性転入出比"]:.3f}  TFR {r["TFR"]:.2f}')

print('=== 図3: 大学収容力と女性転入出比の散布図（実再現・近似）===')
m3 = d[['大学収容力', '女性転入出比', 'Prefecture']].dropna()
r3, p3 = stats.pearsonr(m3['大学収容力'], m3['女性転入出比'])
print(f'大学収容力 × 女性転入出比 の相関 r={r3:+.3f}（p={p3:.4f}, N={len(m3)}）')
print('大学収容力が高い5県:')
for _, r in m3.nlargest(5, '大学収容力').iterrows():
    print(f'  {r["Prefecture"]:6s} 大学生/千人 {r["大学収容力"]:.1f}  転入出比 {r["女性転入出比"]:.3f}')

print('=== 報告値: TFRモデル重要特徴量（原論文 図6・報告値の可視化）===')
CAT_COLOR = {'地理': '#1565C0', '人口構成': '#EF6C00',
             '家族': '#6A1B9A', '就業・教育': '#2E7D32'}
TFR_FEATS = [
    ('経度', '地理'), ('緯度', '地理'),
    ('東京からの距離', '地理'), ('大阪からの距離', '地理'),
    ('15歳未満人口割合', '人口構成'), ('15〜64歳人口割合', '人口構成'),
    ('男性比率', '人口構成'), ('20〜24歳人口割合', '人口構成'),
    ('単独世帯割合', '家族'), ('核家族世帯割合', '家族'),
    ('男性就業率', '就業・教育'), ('最終学歴人口（大学・大学院）', '就業・教育'),
    ('教育業界の従業者数', '就業・教育'),
]
for name, cat in TFR_FEATS:
    print(f'  [{cat:5s}] {name}')
print('平均決定係数 R^2 = 0.78（上位20特徴量・20回交差検証, 原論文の報告値）')

print('=== 報告値: 若年女性転入出モデル重要特徴量（原論文 図8・報告値の可視化）===')
CAT_COLOR2 = dict(CAT_COLOR)
CAT_COLOR2['社会福祉'] = '#00838F'
CAT_COLOR2['産業'] = '#8D6E63'
YW_FEATS = [
    ('65歳以上人口割合', '人口構成'), ('15歳以上人口割合', '人口構成'),
    ('単独世帯数割合', '家族'),
    ('生徒当たり小学校教員数', '就業・教育'), ('大学（院）卒人口割合', '就業・教育'),
    ('生徒当たり小学校数', '就業・教育'),
    ('0〜3歳人口当たり保育所数', '社会福祉'), ('人口当たり一般診療所数', '社会福祉'),
    ('老人福祉費割合', '社会福祉'),
    ('第3次産業従業者割合', '産業'), ('第2次産業事業所数割合', '産業'),
    ('東京からの距離', '地理'), ('傾斜', '地理'),
]
for name, cat in YW_FEATS:
    print(f'  [{cat:5s}] {name}')
print('平均決定係数 R^2 = 0.63（上位25特徴量・20回交差検証, 原論文の報告値）')
print('※ TFRモデルで非常に高かった地理的要因（経度・緯度・距離）が、'
      'こちらでは東京からの距離・傾斜が「多少」効く程度に後退している点が対照的。')