我が国における人口増減の決定要因

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)

# ===== ステップ1: SSDSE-B-2026 の読み込みと人口増減率の計算 =====
df = pd.read_csv(DATA_B, encoding='cp932', header=1)
df = df[df['地域コード'].str.match(r'^R\d{5}', na=False)].copy()
df['年度'] = df['年度'].astype(int)

# 人口増減率 = 自然増減率 + 社会増減率（原論文 3節の定義そのまま）
#   自然増減率 = (出生数 - 死亡数) / 日本人人口 * 1000
#   社会増減率 = (転入者数 - 転出者数) / 日本人人口 * 1000
df['自然増減率'] = (df['出生数'] - df['死亡数']) / df['日本人人口'] * 1000
df['社会増減率'] = (df['転入者数（日本人移動者）'] -
                    df['転出者数（日本人移動者）']) / df['日本人人口'] * 1000
df['人口増減率'] = df['自然増減率'] + df['社会増減率']

# 説明変数（SSDSE-B-2026 で計算可能な11指標。医師数などは未収録）
df['老年人口割合'] = df['65歳以上人口'] / df['総人口'] * 100
for col in ['幼稚園数', '小学校数', '中学校数', '高等学校数',
            '大学数', '一般病院数', '一般診療所数', '歯科診療所数']:
    df[col + '_10万対'] = df[col] / df['総人口'] * 100000  # 人口10万人あたり
df['中学校卒業者進学率'] = df['中学校卒業者のうち進学者数'] / df['中学校卒業者数'] * 100
df['高等学校卒業者進学率'] = df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数'] * 100

print(f"【SSDSE-B-2026】47都道府県 × {df['年度'].nunique()}年度"
      f"（{df['年度'].min()}〜{df['年度'].max()}）＝ {len(df)}観測")
print('  ※ 原論文の SSDSE-2019B（2005〜2016年度）とは収録期間が異なる')

# ===== ステップ2: 図1 主要指標の経年変化（全国平均±標準偏差、実データ追試） =====
# 原論文 図1-(a)〜(o) と同じ「全国平均値の推移＋標準偏差」の描き方を、
# SSDSE-B-2026 で計算できる4指標について 2012〜2023年度で追試する。
items = [('老年人口割合', '老年人口割合 [%]'),
         ('高等学校卒業者進学率', '高等学校卒業者進学率 [%]'),
         ('一般診療所数_10万対', '一般診療所数（人口10万人あたり）'),
         ('人口増減率', '人口増減率 [人口千対]')]
g = df.groupby('年度')
fig, axes = plt.subplots(2, 2, figsize=(10, 7))
for ax, (col, label) in zip(axes.ravel(), items):
    m, s = g[col].mean(), g[col].std()
    ax.errorbar(m.index, m.values, yerr=s.values, color='#455A64',
                ecolor='#C62828', elinewidth=1, capsize=3, lw=1.5)
    ax.set_title(label, fontsize=11)
    ax.set_xlabel('年度')
    ax.grid(alpha=0.3)
fig.suptitle('図1（追試）主要指標の経年変化：47都道府県の平均±標準偏差\n'
             '（SSDSE-B-2026、2012〜2023年度。原論文 図1 は 2005〜2016年度）',
             fontsize=12)
fig.tight_layout(rect=[0, 0, 1, 0.93])
fig.savefig(f'{FIG_DIR}/2019_U2_fig1.png', bbox_inches='tight')
plt.close(fig)
print('【図1】全国平均の変化（2012年度 → 2023年度、実データ）')
for col, label in items:
    print(f'  {label}: {g[col].mean().iloc[0]:.2f} → {g[col].mean().iloc[-1]:.2f}')

# ===== ステップ3: 図2 原論文 表2（ランダム効果モデル）の報告値の可視化 =====
# 以下の数値は原論文 表2 の転記（回帰係数・95%信頼区間・P値）であり、再計算ではない。
paper_re = [  # (変数名, 係数, CI下限, CI上限, P値表記)
    ('老年人口割合 [%]',              -0.656, -0.772, -0.540, '<0.01'),
    ('幼稚園数 (3-5歳人口10万対)',     0.002,  0.000,  0.004, '0.04'),
    ('小学校数 (6-11歳人口10万対)',   -0.002, -0.008,  0.004, '0.53'),
    ('中学校数 (12-14歳人口10万対)',   0.001, -0.007,  0.009, '0.77'),
    ('高等学校数 (15-17歳人口10万対)', -0.024, -0.036, -0.011, '<0.01'),
    ('大学数 (人口10万対)',           -0.074, -1.833,  1.685, '0.94'),
    ('中学校卒業者進学率',            -0.284, -0.557, -0.011, '0.04'),
    ('高等学校卒業者進学率',           0.125,  0.068,  0.182, '<0.01'),
    ('一般病院数 (人口10万対)',       -0.119, -0.360,  0.122, '0.33'),
    ('一般診療所数 (人口10万対)',     -0.090, -0.134, -0.046, '<0.01'),
    ('歯科診療所数 (人口10万対)',      0.110,  0.010,  0.209, '0.03'),
    ('精神科病院数 (人口10万対)',      0.449, -0.812,  1.711, '0.49'),
    ('医師数 (人口10万対)',            0.030,  0.014,  0.045, '<0.01'),
    ('歯科医師数 (人口10万対)',       -0.024, -0.068,  0.019, '0.28'),
    ('看護師・准看護師数 (人口10万対)', 0.000, -0.003,  0.003, '0.92'),
]
fig, ax = plt.subplots(figsize=(9.6, 7.2))
ys = np.arange(len(paper_re))[::-1]
for y, (name, b, lo, hi, p) in zip(ys, paper_re):
    sig = p.startswith('<') or float(p) < 0.05   # 原論文の報告P値で5%有意判定
    color = ('#1565C0' if b > 0 else '#C62828') if sig else '#9E9E9E'
    ax.plot([lo, hi], [y, y], color=color, lw=2)
    ax.plot(b, y, 'o', color=color, ms=7)
    ax.text(1.80, y, f'P{p}' if p.startswith('<') else f'P={p}',
            va='center', fontsize=9, color='#333')
ax.axvline(0, color='#555', lw=0.8, ls='--')
ax.set_yticks(ys)
ax.set_yticklabels([r[0] for r in paper_re], fontsize=10)
ax.set_xlim(-2.0, 2.35)
ax.set_xlabel('回帰係数（人口増減率［人口千対］への影響、95%信頼区間）')
ax.set_title('図2：ランダム効果モデルによる解析結果\n'
             '（原論文 表2 の報告値の可視化。再計算ではない）\n'
             '青＝正で有意 / 赤＝負で有意 / 灰＝有意でない（原論文の報告P値で判定）', fontsize=11)
ax.grid(axis='x', alpha=0.3)
fig.tight_layout()
fig.savefig(f'{FIG_DIR}/2019_U2_fig2.png', bbox_inches='tight')
plt.close(fig)
is_sig = lambda p: p.startswith('<') or float(p) < 0.05
sig_pos = [r[0].split(' ')[0] for r in paper_re if is_sig(r[4]) and r[1] > 0]
sig_neg = [r[0].split(' ')[0] for r in paper_re if is_sig(r[4]) and r[1] < 0]
print('【図2】原論文 表2（ランダム効果モデル）の報告値')
print('  正で有意:', '、'.join(sig_pos))
print('  負で有意:', '、'.join(sig_neg))

# ===== ステップ6: 追試 — 11変数によるパネルデータ分析（FE・RE、期間・変数が異なる） =====
# 原論文と同型のモデル Y_it = α + Σβj X_j,it + u_i + ε_it を、SSDSE-B-2026 で
# 計算可能な11変数・2012〜2023年度で推定する。医師数など4変数は未収録のため入らず、
# 学校数の分母も対象年齢人口ではなく総人口。原論文の係数の「検算」ではない点に注意。
xvars = ['老年人口割合', '幼稚園数_10万対', '小学校数_10万対', '中学校数_10万対',
         '高等学校数_10万対', '大学数_10万対', '中学校卒業者進学率',
         '高等学校卒業者進学率', '一般病院数_10万対', '一般診療所数_10万対',
         '歯科診療所数_10万対']
pdat = df[['都道府県', '年度', '人口増減率'] + xvars].dropna()
y, X = pdat['人口増減率'], sm.add_constant(pdat[xvars])

# (1) プーリングOLS と 固定効果（LSDV）→ 個体効果の存在を F 検定
pool = sm.OLS(y, X).fit()
Xd = pd.concat([X, pd.get_dummies(pdat['都道府県'], drop_first=True, dtype=float)], axis=1)
lsdv = sm.OLS(y, Xd).fit()
dfn, dfd = pdat['都道府県'].nunique() - 1, int(lsdv.df_resid)
F = ((lsdv.rsquared - pool.rsquared) / dfn) / ((1 - lsdv.rsquared) / dfd)
print(f'【追試】個体効果の F 検定: F({dfn}, {dfd}) = {F:.2f}, '
      f'P = {stats.f.sf(F, dfn, dfd):.1e}（原論文: F(46,502)=13.82, P<0.001）')

# (2) Breusch-Pagan の LM 検定（プーリングOLS残差から。ランダム効果の存在を検定）
e = pool.resid
T = pdat['年度'].nunique()
S1 = pdat.assign(e=e).groupby('都道府県')['e'].sum().pow(2).sum()
LM = len(pdat) * T / (2 * (T - 1)) * (S1 / (e ** 2).sum() - 1) ** 2
print(f'  Breusch-Pagan LM 検定: LM = {LM:.1f}, '
      f'P = {stats.chi2.sf(LM, 1):.1e}（原論文: P<0.001 → ランダム効果あり）')

# (3) ランダム効果（変量切片）モデル: MixedLM で推定
re = sm.MixedLM(y, X, groups=pdat['都道府県']).fit(reml=True)
ci = re.conf_int()
print('\n【追試】ランダム効果モデル（SSDSE-B-2026、2012〜2023年度、11変数）')
print('  ※ 期間・変数定義が異なるため原論文 表2 と数値は一致しない')
paper_map = {r[0].split(' ')[0]: r for r in paper_re}
print(f"  {'変数':　<11s}{'係数':>8s} {'95%信頼区間':>16s} {'P値':>6s}   原論文表2の係数")
for v in xvars:
    key = v.replace('_10万対', '')
    pb = paper_map[key]
    ptxt = f'P{pb[4]}' if pb[4].startswith('<') else f'P={pb[4]}'
    star = '*' if re.pvalues[v] < 0.05 else ' '
    print(f'  {key:　<11s}{re.params[v]:+8.3f} '
          f'[{ci.loc[v, 0]:+.3f}, {ci.loc[v, 1]:+.3f}] {re.pvalues[v]:6.3f}{star}  '
          f'{pb[1]:+.3f} ({ptxt})')
print('  （* は追試で5%有意。右端の欄は原論文 表2 の報告値）')

# ===== ステップ4: 図3 人口増加が続いた5都道府県の高等学校卒業者進学率（実データ追試） =====
# 原論文 図2: 2005〜2016年度に人口増減率が常に正だった5都道府県
# （埼玉・東京・神奈川・愛知・沖縄）の高等学校卒業者進学率の推移。
# ここでは同じ5都道府県を SSDSE-B-2026（2012〜2023年度）で追試する。
five = ['埼玉県', '東京都', '神奈川県', '愛知県', '沖縄県']
colors5 = ['#1565C0', '#C62828', '#2E7D32', '#6A1B9A', '#E65100']
fig, ax = plt.subplots(figsize=(9, 5))
for pref, c in zip(five, colors5):
    sub = df[df['都道府県'] == pref].sort_values('年度')
    ax.plot(sub['年度'], sub['高等学校卒業者進学率'], marker='o', ms=4, color=c, label=pref)
ax.set_xlabel('年度')
ax.set_ylabel('高等学校卒業者進学率 [%]')
ax.set_title('図3（追試）人口増加が続いた5都道府県の高等学校卒業者進学率\n'
             '（SSDSE-B-2026、2012〜2023年度。原論文 図2 は 2005〜2016年度）', fontsize=11)
ax.legend(loc='lower right', fontsize=10)
ax.grid(alpha=0.3)
fig.tight_layout()
fig.savefig(f'{FIG_DIR}/2019_U2_fig3.png', bbox_inches='tight')
plt.close(fig)
print('【図3】5都道府県の高等学校卒業者進学率（2012→2023年度、実データ）')
for pref in five:
    sub = df[df['都道府県'] == pref].sort_values('年度')
    print(f"  {pref}: {sub['高等学校卒業者進学率'].iloc[0]:.1f}% → "
          f"{sub['高等学校卒業者進学率'].iloc[-1]:.1f}%")

# ===== ステップ5: 図4 2023年度の人口増減率ランキング＋「常に正」の検証（実データ） =====
always_pos = [p for p, s in df.groupby('都道府県')['人口増減率'] if (s > 0).all()]
print('【検証】2012〜2023年度に人口増減率が常に正だった都道府県:',
      '、'.join(always_pos) if always_pos else 'なし')
latest = df[df['年度'] == 2023].sort_values('人口増減率')
fig, ax = plt.subplots(figsize=(8.5, 11))
cols = ['#1565C0' if v > 0 else ('#E65100' if p in five else '#B0BEC5')
        for v, p in zip(latest['人口増減率'], latest['都道府県'])]
ax.barh(range(len(latest)), latest['人口増減率'], color=cols, height=0.65)
ax.set_yticks(range(len(latest)))
ax.set_yticklabels(latest['都道府県'], fontsize=9)
ax.axvline(0, color='#555', lw=0.8)
ax.set_xlabel('人口増減率 [人口千対]（自然増減率＋社会増減率）')
ax.set_title('図4：2023年度の都道府県別人口増減率（SSDSE-B-2026 から計算した実データ）\n'
             '青＝プラス / オレンジ＝原論文の「12年間常に正」だった5都道府県のうちマイナスの県', fontsize=11)
ax.grid(axis='x', alpha=0.3)
fig.tight_layout()
fig.savefig(f'{FIG_DIR}/2019_U2_fig4.png', bbox_inches='tight')
plt.close(fig)
top3 = latest.tail(3)[['都道府県', '人口増減率']].iloc[::-1]
print('【図4】2023年度の人口増減率 上位3:',
      ' / '.join(f'{r.都道府県} {r.人口増減率:+.2f}' for r in top3.itertuples()))
print('  原論文の5都道府県の2023年度値:',
      ' / '.join(f"{p} {latest.set_index('都道府県').loc[p, '人口増減率']:+.2f}" for p in five))