マルチレベル分析を用いた
市町村大学等進学率の決定要因分析

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
import statsmodels.formula.api as smf

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

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

# ===== ステップ2: 図1 都道府県別の高等学校卒業者進学率（SSDSE-B 実データ） =====
# 原論文 図1 は「都道府県別大学等進学率 × 全国学力調査正答率」の散布図
# （相関係数 0.226 ＝ 中学時代の学力では進学率格差をほぼ説明できない）。
# 正答率は SSDSE 未収録のため散布図は原論文参照とし、ここでは格差の存在
# そのもの（研究の出発点）を最新の実データで確認する。
dfB = pd.read_csv(DATA_B, encoding='cp932', header=1)
dfB = dfB[dfB['地域コード'].str.match(r'^R\d{5}', na=False)].copy()
latest_y = dfB['年度'].max()
b = dfB[dfB['年度'] == latest_y].copy()
b['高等学校卒業者進学率'] = b['高等学校卒業者のうち進学者数'] / b['高等学校卒業者数'] * 100
b = b.sort_values('高等学校卒業者進学率')
fig, ax = plt.subplots(figsize=(8.5, 11))
cols = ['#C62828' if v >= b['高等学校卒業者進学率'].nlargest(5).min()
        else ('#1565C0' if v <= b['高等学校卒業者進学率'].nsmallest(5).max() else '#90A4AE')
        for v in b['高等学校卒業者進学率']]
ax.barh(range(len(b)), b['高等学校卒業者進学率'], color=cols, height=0.65)
ax.set_yticks(range(len(b)))
ax.set_yticklabels(b['都道府県'], fontsize=9)
mean_v = b['高等学校卒業者進学率'].mean()
ax.axvline(mean_v, color='#555', lw=1, ls='--')
ax.text(mean_v + 0.3, 0.5, f'47都道府県平均 {mean_v:.1f}%', fontsize=9, color='#555')
ax.set_xlabel('高等学校卒業者進学率 [%]（大学・短大等への進学）')
ax.set_title(f'図1：都道府県別 高等学校卒業者進学率（{latest_y}年度、SSDSE-B-2026 実データ）\n'
             '赤＝上位5 / 青＝下位5。原論文の問題意識である「進学率の地域格差」は現在も大きい\n'
             '（原論文 図1〔進学率×全国学力調査正答率の散布図、相関係数0.226〕は原論文参照）',
             fontsize=11)
ax.grid(axis='x', alpha=0.3)
fig.tight_layout()
fig.savefig(f'{FIG_DIR}/2019_U3_fig1.png', bbox_inches='tight')
plt.close(fig)
gap = b['高等学校卒業者進学率'].max() - b['高等学校卒業者進学率'].min()
print(f'【図1】{latest_y}年度の高等学校卒業者進学率（実データ）: '
      f'最高 {b.iloc[-1]["都道府県"]} {b["高等学校卒業者進学率"].max():.1f}% / '
      f'最低 {b.iloc[0]["都道府県"]} {b["高等学校卒業者進学率"].min():.1f}% '
      f'（差 {gap:.1f}ポイント）')

# ===== ステップ5: 集団平均中心化を実データでやってみる（SSDSE-A） =====
# 原論文は市町村レベル変数を「集団平均中心化」（各市町村の値から所属都道府県の
# 平均を引く）してから投入した。都道府県間の差を取り除き、係数を「同じ県の中で
# 値が高い市町村ほど…」という市町村レベルの効果として解釈できるようにする操作。
dfA['知識率_県平均'] = dfA.groupby('都道府県')['知識集約型産業従事者率'].transform('mean')
dfA['知識率_集団平均中心化'] = dfA['知識集約型産業従事者率'] - dfA['知識率_県平均']
chk = dfA.groupby('都道府県')['知識率_集団平均中心化'].mean().abs().max()
print('【集団平均中心化】知識集約型産業従事者率（SSDSE-A 実データ）')
print(f'  中心化後の都道府県平均の最大絶対値 = {chk:.1e} （どの県でも平均0になった）')
print(f"  例）東京都 千代田区: 元の値 {dfA.loc[dfA['市区町村']=='千代田区','知識集約型産業従事者率'].iloc[0]:.1f}%"
      f" → 中心化後 {dfA.loc[dfA['市区町村']=='千代田区','知識率_集団平均中心化'].iloc[0]:+.1f}")

# ===== ステップ1: SSDSE-A（市区町村）から原論文と同型の説明変数を作る =====
# 原論文 表3 の加工方法に従う（経済センサスの産業分類は同一。年次は異なる）。
#   知識集約型産業従事者率 = (情報通信＋金融保険＋不動産物品賃貸＋
#                             学術研究専門技術＋教育学習支援) / 従業者総数 ×100
#   人口密度 [人/km2] = 総人口 / 総面積（SSDSE-A の面積単位は ha → /100 で km2）
#   65歳以上人口比率 = 65歳以上人口 / 総人口 ×100
dfA = pd.read_csv(DATA_A, encoding='cp932', skiprows=2)
knowledge = ['従業者数（民営）（情報通信業）', '従業者数（民営）（金融業、保険業）',
             '従業者数（民営）（不動産業、物品賃貸業）',
             '従業者数（民営）（学術研究、専門・技術サービス業）',
             '従業者数（民営）（教育、学習支援業）']
dfA = dfA[dfA['従業者数（民営）'] > 0].copy()
dfA['知識集約型産業従事者率'] = dfA[knowledge].sum(axis=1) / dfA['従業者数（民営）'] * 100
dfA['人口密度'] = dfA['総人口'] / (dfA['総面積（北方地域及び竹島を除く）'] / 100)
dfA['６５歳以上人口比率'] = dfA['65歳以上人口'] / dfA['総人口'] * 100
print(f"【SSDSE-A-2025】{len(dfA)}市区町村 × {dfA['都道府県'].nunique()}都道府県の階層データ")
print('  ※ 原論文は高等学校卒業生が存在する 1304市町村（学校基本調査ベース）')
print('  ※ 被説明変数「市町村別大学等進学率」と「一人当たり課税対象所得」は')
print('     SSDSE-A 未収録 → 原論文のモデルの再推定は不可（表7・表8 は報告値で提示）')
q = dfA[['知識集約型産業従事者率', '人口密度', '６５歳以上人口比率']].describe().round(2)
print(q.loc[['min', '25%', '50%', 'mean', '75%', 'max']].to_string())

# ===== ステップ6: 図4 級内相関(ICC)を実データで体感する（手法のデモ） =====
# 市区町村の「知識集約型産業従事者率」が都道府県というグループでどれだけ
# 塊になっているかを、(1)箱ひげ図 と (2)ヌルモデル（ランダム切片のみの
# マルチレベルモデル）の ICC で確かめる。原論文のモデル1（被説明変数は
# 大学等進学率、ICC=0.244）と同じ手順を、SSDSE-A で使える別の変数に適用
# した「手法のデモ」であり、原論文の数値の検算ではない。
null_model = smf.mixedlm('知識集約型産業従事者率 ~ 1', dfA, groups=dfA['都道府県'])
res = null_model.fit(reml=True)
var_u = res.cov_re.iloc[0, 0]          # 都道府県レベルの分散
var_e = res.scale                      # 市区町村レベルの残差分散
icc = var_u / (var_u + var_e)
print('【図4】ヌルモデル（ランダム切片のみ）による級内相関 — SSDSE-A 実データ')
print(f'  都道府県レベル分散 = {var_u:.2f} / 市区町村レベル分散 = {var_e:.2f}')
print(f'  ICC = {var_u:.2f} / ({var_u:.2f} + {var_e:.2f}) = {icc:.3f}')
print(f'  → 知識集約型産業従事者率のばらつきの約{icc*100:.1f}%は都道府県間の差')
print('  （原論文の ICC=0.244 は「市町村別大学等進学率」についての報告値。別の変数）')

order = dfA.groupby('都道府県')['知識集約型産業従事者率'].median().sort_values().index
data = [dfA.loc[dfA['都道府県'] == p, '知識集約型産業従事者率'].values for p in order]
fig, ax = plt.subplots(figsize=(9, 11))
bp = ax.boxplot(data, vert=False, patch_artist=True, showfliers=False, widths=0.6)
for patch in bp['boxes']:
    patch.set_facecolor('#BBDEFB'); patch.set_edgecolor('#1565C0')
for med in bp['medians']:
    med.set_color('#C62828')
ax.set_yticklabels(order, fontsize=9)
ax.set_xlabel('市区町村の知識集約型産業従事者率 [%]（原論文 表3 と同じ定義で計算）')
ax.set_title('図4：級内相関を体感する — 市区町村の知識集約型産業従事者率の都道府県別分布\n'
             f'（SSDSE-A-2025 実データによる手法デモ。ヌルモデルの ICC = {icc:.3f}）\n'
             '箱の位置が県ごとに系統的にずれる＝都道府県レベルの要因が存在する',
             fontsize=11)
ax.grid(axis='x', alpha=0.3)
fig.tight_layout()
fig.savefig(f'{FIG_DIR}/2019_U3_fig4.png', bbox_inches='tight')
plt.close(fig)
top3 = dfA.groupby('都道府県')['知識集約型産業従事者率'].median().nlargest(3)
print('  都道府県中央値の上位3:', ' / '.join(f'{k} {v:.1f}%' for k, v in top3.items()))

# ===== ステップ3: 図2 原論文 表7 の係数の可視化（報告値。再計算ではない） =====
# (変数名, モデル3係数, SE, モデル4係数, SE, レベル)  … 原論文 表7 の転記
paper_t7 = [
    ('一人当たり課税対象所得(十万円)', 0.46, 0.14, 1.16, 0.24, '市町村'),
    ('知識集約型産業従事者率(%)',      0.67, 0.14, 0.63, 0.14, '市町村'),
    ('人口密度',                       0.09, 0.03, 0.09, 0.03, '市町村'),
    ('65歳以上人口比率(%)',           -0.64, 0.09, -0.48, 0.09, '市町村'),
    ('都道府県別大学収容率(%)',        0.16, 0.04, 0.14, 0.03, '都道府県'),
    ('min(東京までの距離,京都までの距離)', -0.02, 0.00, -0.02, 0.00, '都道府県'),
]
fig, ax = plt.subplots(figsize=(9.6, 5.6))
ys = np.arange(len(paper_t7))[::-1]
for y, (name, b3, se3, b4, se4, lv) in zip(ys, paper_t7):
    for b_, se_, dy, c, m in [(b3, se3, 0.16, '#1565C0', 'o'), (b4, se4, -0.16, '#E65100', 's')]:
        ax.plot([b_ - 1.96 * se_, b_ + 1.96 * se_], [y + dy, y + dy], color=c, lw=2)
        ax.plot(b_, y + dy, m, color=c, ms=7)
ax.axvline(0, color='#555', lw=0.8, ls='--')
ax.set_yticks(ys)
ax.set_yticklabels([f"〔{r[5]}〕{r[0]}" for r in paper_t7], fontsize=10)
ax.plot([], [], 'o-', color='#1565C0', label='モデル3（ランダム切片）')
ax.plot([], [], 's-', color='#E65100', label='モデル4（ランダム切片・傾き）')
ax.legend(loc='lower right', fontsize=10)
ax.set_xlabel('大学等進学率への係数（点＝係数、横線＝係数±1.96×標準誤差）')
ax.set_title('図2：マルチレベル分析の係数（原論文 表7 の報告値の可視化。再計算ではない）\n'
             'いずれの係数もモデル3・4で5%水準より強い有意（所得は**、他は***）\n'
             '市町村レベル変数は集団平均中心化後、都道府県レベル変数は全体平均中心化後の値',
             fontsize=11)
ax.grid(axis='x', alpha=0.3)
fig.tight_layout()
fig.savefig(f'{FIG_DIR}/2019_U3_fig2.png', bbox_inches='tight')
plt.close(fig)
print('【図2】原論文 表7 の報告値（モデル3 → モデル4）')
for name, b3, se3, b4, se4, lv in paper_t7:
    print(f'  〔{lv}〕{name}: {b3:+.2f} → {b4:+.2f}')
print('  ※ モデル4では所得の傾きにランダム効果（分散0.52、切片との共分散 -1.19）')

# ===== ステップ7: 原論文 表7 のモデル比較指標（報告値の整理） =====
print('\n【原論文 表7 の報告値】モデルの当てはまり（AIC・BIC は小さいほどよい）')
t7 = pd.DataFrame({
    'AIC': [11260.38, 11038.80, 11023.46, 11011.03],
    'BIC': [11275.90, 11075.01, 11070.02, 11067.94],
    'LogLik': [-5627.19, -5512.40, -5502.73, -5494.52],
}, index=['モデル1(ヌル)', 'モデル2(市町村)', 'モデル3(＋都道府県)', 'モデル4(＋傾き)'])
print(t7.to_string())
print('  → モデル4が最良: 都道府県間で「所得の傾き」が異なるモデルの当てはまりが最もよい')
print('  → 切片と所得傾きのランダム効果の共分散 -1.19（負）: 都道府県要因の切片が')
print('     大きい県ほど、所得が進学率に与える影響が小さい（原論文 4.2-1）')

# ===== ステップ4: 図3 都道府県切片分散の比較（原論文 表7・表8 の報告値） =====
# 「採択した都道府県レベル変数が、都道府県要因（切片のばらつき）をどれだけ
#   説明したか」を Var: prefecture (Intercept) の減少で確認する（原論文 4.2-2）。
models = ['モデル1\n(ヌル)', 'モデル2\n(＋市町村変数)', 'モデル3-1\n(＋大学収容率)',
          'モデル3-2\n(＋距離)', 'モデル3\n(＋両方)', 'モデル4\n(＋所得の傾き)']
var_pref = [66.42, 69.24, 42.72, 46.38, 29.10, 29.91]   # 表7・表8の報告値
var_res  = [306.79, 254.34, 254.34, 254.32, 254.40, 247.64]
fig, axes = plt.subplots(1, 2, figsize=(11, 4.6))
colors = ['#78909C', '#78909C', '#42A5F5', '#42A5F5', '#1565C0', '#E65100']
axes[0].bar(models, var_pref, color=colors)
for i, v in enumerate(var_pref):
    axes[0].text(i, v + 1, f'{v:.1f}', ha='center', fontsize=10)
axes[0].set_ylabel('Var: prefecture (Intercept)')
axes[0].set_title('都道府県切片の分散（都道府県要因の残り）', fontsize=11)
axes[1].bar(models, var_res, color=colors)
for i, v in enumerate(var_res):
    axes[1].text(i, v + 3, f'{v:.0f}', ha='center', fontsize=10)
axes[1].set_ylabel('Var: Residual')
axes[1].set_title('市町村レベルの残差分散', fontsize=11)
for ax in axes:
    ax.tick_params(axis='x', labelsize=8.5)
    ax.grid(axis='y', alpha=0.3)
fig.suptitle('図3：分散成分の比較（原論文 表7・表8 の報告値の可視化。再計算ではない）\n'
             '大学収容率と距離を両方入れると都道府県切片の分散は 69.24 → 29.10 に減少（それでも残る）',
             fontsize=11)
fig.tight_layout(rect=[0, 0, 1, 0.90])
fig.savefig(f'{FIG_DIR}/2019_U3_fig3.png', bbox_inches='tight')
plt.close(fig)
print('【図3】都道府県切片の分散（原論文の報告値）: '
      f'モデル2 {var_pref[1]} → モデル3-1 {var_pref[2]} → モデル3-2 {var_pref[3]} → モデル3 {var_pref[4]}')