若者の大都市から地方への移動要因を探る
―修正重力モデルによる分析―

print('=== SSDSE-A 読み込みと説明変数の算出【実再現】 ===')
raw = pd.read_csv(DATA_A, encoding='cp932', header=2)
raw = raw.rename(columns={raw.columns[0]: '地域コード',
                          raw.columns[1]: '都道府県',
                          raw.columns[2]: '市区町村'})
df = raw[raw['地域コード'].astype(str).str.match(r'^R\d{5}$', na=False)].copy()

num_cols = ['総人口', '就業者数', '第１次産業就業者数', '第２次産業就業者数',
            '第３次産業就業者数', '完全失業者数', '幼稚園数', '小学校数',
            '中学校数', '高等学校数', '公民館数', '図書館数', '一般病院数',
            '保育所等数（基本票）']
for c in num_cols:
    df[c] = pd.to_numeric(df[c], errors='coerce')
# 総人口0・就業者0の市町村（2020年国勢調査で人口0の避難区域等）は率が定義できないため除外
df = df[(df['総人口'] > 0) & (df['就業者数'] > 0)].copy()

# 原論文と同じ変数の作り方（3.3 使用した説明変数）を SSDSE-A-2025 に適用
df['第1次産業就業者率'] = df['第１次産業就業者数'] / df['就業者数'] * 100
df['第2次産業就業者率'] = df['第２次産業就業者数'] / df['就業者数'] * 100
df['第3次産業就業者率'] = df['第３次産業就業者数'] / df['就業者数'] * 100
df['完全失業率'] = df['完全失業者数'] / (df['完全失業者数'] + df['就業者数']) * 100
df['学校数_10万人'] = (df['幼稚園数'] + df['小学校数'] + df['中学校数'] + df['高等学校数']) / df['総人口'] * 1e5
df['福祉施設数_10万人'] = (df['公民館数'] + df['図書館数'] + df['一般病院数'] + df['保育所等数（基本票）']) / df['総人口'] * 1e5

# 大都市（移動元）＝東京23特別区＋政令指定都市20市、地方（移動先）＝総人口20万人未満
seirei = ['札幌市', '仙台市', 'さいたま市', '千葉市', '横浜市', '川崎市', '相模原市',
          '新潟市', '静岡市', '浜松市', '名古屋市', '京都市', '大阪市', '堺市',
          '神戸市', '岡山市', '広島市', '北九州市', '福岡市', '熊本市']
is_ward = (df['都道府県'] == '東京都') & df['市区町村'].str.endswith('区', na=False)
is_seirei = df['市区町村'].isin(seirei)
df['地域区分'] = np.where(is_ward | is_seirei, '大都市（移動元）',
                       np.where(df['総人口'] < 200000, '地方（移動先）', 'その他（中核市等）'))

n_big = (df['地域区分'] == '大都市（移動元）').sum()
n_local = (df['地域区分'] == '地方（移動先）').sum()
print(f'市区町村総数: {len(df)}')
print(f'大都市（移動元）: {n_big} 市区  ／  地方（移動先）: {n_local} 市町村  ／  その他（中核市等）: {(df["地域区分"]=="その他（中核市等）").sum()}')
print('※ 原論文（平成27年国勢調査）では 43 市区 / 1606 市町村。収録版は2020年国勢調査ベースのため件数は近似。')

show_vars = ['総人口', '第1次産業就業者率', '第2次産業就業者率', '第3次産業就業者率',
             '完全失業率', '学校数_10万人', '福祉施設数_10万人']
for grp in ['大都市（移動元）', '地方（移動先）']:
    sub = df[df['地域区分'] == grp]
    print(f'\n【{grp}】 N={len(sub)}')
    stat = sub[show_vars].agg(['mean', 'std', 'min', 'max']).T
    stat.columns = ['平均', '標準偏差', '最小', '最大']
    print(stat.round(2).to_string())
print()

corr_vars = ['総人口', '第1次産業就業者率', '第2次産業就業者率', '第3次産業就業者率',
             '完全失業率', '学校数_10万人', '福祉施設数_10万人']
short = {'総人口': '総人口', '第1次産業就業者率': '第1次産業', '第2次産業就業者率': '第2次産業',
         '第3次産業就業者率': '第3次産業', '完全失業率': '完全失業率',
         '学校数_10万人': '学校数', '福祉施設数_10万人': '福祉施設数'}
cmat = loc[corr_vars].corr()
cmat_disp = cmat.rename(index=short, columns=short)
print('（地方＝移動先市町村 N={} での相関）'.format(len(loc)))
print(cmat_disp.round(2).to_string())
# 最も強い相関ペア
pairs = []
cols = corr_vars
for i in range(len(cols)):
    for j in range(i + 1, len(cols)):
        pairs.append((abs(cmat.iloc[i, j]), short[cols[i]], short[cols[j]], cmat.iloc[i, j]))
pairs.sort(reverse=True)
print('\n相関の強い変数ペア上位3組（多重共線性の目安）:')
for a, x, y, v in pairs[:3]:

reported = [
    ('総人口＿x',                              0.061648),
    ('総人口＿y',                              0.069374),
    ('移動間距離',                             -0.568857),
    ('第1次産業就業者率＿x',                    0.134868),
    ('第1次産業就業者率＿y',                   -0.029745),
    ('第2次産業就業者率＿x',                   -0.060645),
    ('第2次産業就業者率＿y',                   -0.130457),
    ('第3次産業就業者率＿x',                    0.078206),
    ('完全失業率＿x',                          0.015544),
    ('完全失業率＿y',                          -0.032623),
    ('学校数＿x',                              0.114437),
    ('学校数＿y',                              0.107113),
    ('福祉施設数＿x',                          0.134630),
    ('福祉施設数＿y',                          0.041808),
    ('子供を対象とした活動＿x',                -0.044284),
    ('子供を対象とした活動＿y',                 0.010989),
    ('スポーツ文化芸術学術活動＿x',            -0.076732),
    ('スポーツ文化芸術学術活動＿y',            -0.031711),
    ('安全な生活のための活動＿x',               0.062168),
    ('安全な生活のための活動＿y',              -0.023933),
    ('課税対象所得＿x',                        -0.086024),
    ('課税対象所得＿y',                        -0.042550),
]
rep = pd.DataFrame(reported, columns=['変数名', '回帰係数'])
print('（原論文 表2 の報告値。x＝移動元、y＝移動先。※移動フロー・距離は SSDSE 未収録のため再計算不可）')
rep_sorted = rep.reindex(rep['回帰係数'].abs().sort_values(ascending=False).index)
for _, r in rep_sorted.iterrows():
    print(f'  {r["変数名"]:<24} {r["回帰係数"]:+.6f}')
print(f'\n絶対値最大：移動間距離 {rep.loc[rep["変数名"]=="移動間距離","回帰係数"].iloc[0]:+.6f}（長距離移動が人口移動のネック）')

pairs_xy = [
    ('総人口', 0.061648, 0.069374),
    ('第1次産業就業者率', 0.134868, -0.029745),
    ('第2次産業就業者率', -0.060645, -0.130457),
    ('完全失業率', 0.015544, -0.032623),
    ('学校数', 0.114437, 0.107113),
    ('福祉施設数', 0.134630, 0.041808),
    ('子供を対象とした活動', -0.044284, 0.010989),
    ('安全な生活のための活動', 0.062168, -0.023933),
    ('課税対象所得', -0.086024, -0.042550),
]
pxy = pd.DataFrame(pairs_xy, columns=['変数', '移動元x', '移動先y'])
print('（同じ指標でも移動元と移動先で符号が変わる項目に注目。原論文 表2 の報告値）')
for _, r in pxy.iterrows():
    flip = '  ← 符号反転' if r['移動元x'] * r['移動先y'] < 0 else ''
    print(f'  {r["変数"]:<20} x={r["移動元x"]:+.4f}  y={r["移動先y"]:+.4f}{flip}')