日本で暮らす外国人の動向から見た多民族化

import pandas as pd

# SSDSE-A: 1行目=項目コード、2行目=データ年度、3行目=日本語項目名
df_a = pd.read_csv('data/raw/SSDSE-A-2025.csv', encoding='cp932', header=2)
df_a = df_a[df_a['市区町村'].notna()].copy()

# 派生変数（原論文と同じ定義）
df_a['外国人比率'] = df_a['外国人人口'] / df_a['総人口'] * 10000   # 1万人あたり
df_a['生産年齢人口比率'] = df_a['15～64歳人口'] / df_a['総人口'] * 100
df_a['人口密度'] = df_a['総人口'] / (df_a['総面積（北方地域及び竹島を除く）'] / 100)

# 都道府県レベル＝市区町村を合算
df_p = df_a.groupby('都道府県', as_index=False)[
    ['総人口', '外国人人口', '15～64歳人口', '事業所数（民営）']].sum()

# 相関係数（Excel の CORREL と同じ Pearson の r）
r_muni = df_a['総人口'].corr(df_a['外国人人口'])
r_pref = df_p['総人口'].corr(df_p['外国人人口'])
print(f"総人口×外国人人口: 市区町村 r={r_muni:.3f} / 都道府県 r={r_pref:.3f}")
print("原論文の報告値   : 市町村   r=0.879 / 都道府県 r=0.932")

import numpy as np
import matplotlib.pyplot as plt

def scatter_two(dm, dp, x, y, xlab, ylab, title, fname, xdiv=1000, ydiv=1000):
    fig, axes = plt.subplots(1, 2, figsize=(12.5, 5.4))
    for ax, d, name in [(axes[0], dm, '市区町村'), (axes[1], dp, '都道府県')]:
        dd = d[[x, y]].dropna()
        xs, ys = dd[x] / xdiv, dd[y] / ydiv
        r = dd[x].corr(dd[y])
        ax.scatter(xs, ys, s=26, color='#3F6BB5', alpha=0.55)
        m, b = np.polyfit(xs, ys, 1)              # 回帰直線（1次式）
        xr = np.linspace(xs.min(), xs.max(), 100)
        ax.plot(xr, m * xr + b, 'k--', linewidth=1.4)
        ax.set_title(f'{name}（n={len(dd)}）  r = {r:.3f},  R² = {r**2:.3f}')
        ax.set_xlabel(xlab); ax.set_ylabel(ylab)
    fig.suptitle(title, fontweight='bold')
    fig.savefig(fname, bbox_inches='tight')

scatter_two(df_a, df_p, '事業所数（民営）', '外国人人口',
            '事業所数（千事業所）', '外国人人口（千人）',
            '図2（再現）：事業所数と外国人人口', '2019_H4_fig2.png')

top10 = df_a.nlargest(10, '外国人比率')[
    ['都道府県', '市区町村', '外国人比率', '総人口', '人口密度']]
for i, (_, row) in enumerate(top10.iterrows(), 1):
    print(f"{i:>2}. {row['都道府県']} {row['市区町村']:<8} "
          f"{row['外国人比率']:>8.2f} /1万人")

# 散布図の上位5点に市町村名を注記
dd = df_a[['人口密度', '外国人比率', '市区町村']].dropna()
top5 = dd.nlargest(5, '外国人比率')
for _, row in top5.iterrows():
    ax.annotate(row['市区町村'], (row['人口密度'], row['外国人比率']),
                xytext=(6, 2), textcoords='offset points',
                fontsize=9, color='#C62828', fontweight='bold')

df_b = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', header=1)
df_b = df_b[df_b['都道府県'] != '全国'].copy()

# 原論文と同じ算出方法：外国人人口 ＝ 総人口 − 日本人人口
df_b['外国人人口'] = df_b['総人口'] - df_b['日本人人口']
df_b['外国人比率'] = df_b['外国人人口'] / df_b['総人口'] * 10000

# 年×都道府県の表に変換して47本の折れ線を描く
pivot = df_b.pivot_table(index='年度', columns='都道府県', values='外国人比率')

# 全国計でスパイクを数値確認
nat = df_b.groupby('年度')[['総人口', '日本人人口']].sum()
nat['外国人比率'] = (nat['総人口'] - nat['日本人人口']) / nat['総人口'] * 10000
for y in nat.index:
    mark = ' ←国勢調査年' if y in (2015, 2020) else ''
    print(f"  {y}: {nat.loc[y, '外国人比率']:>6.1f}{mark}")

df_p['外国人比率'] = df_p['外国人人口'] / df_p['総人口'] * 10000
rank = df_p.sort_values('外国人比率', ascending=False).reset_index(drop=True)

for i in range(6):
    print(f"{i+1}位 {rank.loc[i, '都道府県']} {rank.loc[i, '外国人比率']:.1f}")
hiro = rank.index[rank['都道府県'] == '広島県'][0] + 1
print(f"…{hiro}位 広島県、47位 {rank.loc[46, '都道府県']}")

# ランキング棒グラフ（塗り分けマップの代替表現）
fig, ax = plt.subplots(figsize=(13.5, 6))
ax.bar(range(47), rank['外国人比率'])
ax.set_xticks(range(47))
ax.set_xticklabels(rank['都道府県'], rotation=60, ha='right', fontsize=8.5)

prefs3 = ['東京都', '群馬県', '広島県']
orig_keijo = {'東京都': -0.365, '群馬県': -0.268, '広島県': -0.164}  # 原論文の報告値
orig_saimu = {'東京都': -0.440, '群馬県': -0.391, '広島県': -0.011}  # 原論文の報告値

for pref in prefs3:
    d = df_a[df_a['都道府県'] == pref]
    for xcol, orig in [('経常収支比率（市町村財政）', orig_keijo),
                       ('実質公債費比率（市町村財政）', orig_saimu)]:
        dd = d[[xcol, '外国人人口']].dropna()
        r = dd[xcol].corr(dd['外国人人口'])
        print(f"  {pref} {xcol[:7]}×外国人人口: "
              f"原論文 {orig[pref]:+.3f} → 再現 {r:+.3f} (n={len(dd)})")