何を確かめたいか

「所得が高い街ほど、人が流入しているのか？」を、全国1,740市区町村の e-Stat データで確かめました。結論は、そのまま比べると関係があるように見えるが、年齢構成などをそろえると0と区別できなくなる、というものです（結果の読み方は note の記事に書いています）。

この記事では、その確かめ方の実装に絞ります。ポイントは次の3つです。

説明変数を標準化して、「1標準偏差あたり何人違うか」でそろえる
条件を1つずつ足した4本の回帰を回し、同じ係数がどう動くかを表にする
表をそのまま係数の階段図にして、0をまたいだら灰色にする

使ったライブラリは pandas・numpy・statsmodels・matplotlib です。

データの取り方と整形

データは e-Stat の「社会・人口統計体系市区町村データ」と「統計でみる市区町村のすがた2026」を、項目コードの列名のまま1つの CSV（data/municipal_base.csv）にまとめたものです。1行が1市区町村で、A1101（総人口）、C120110（課税対象所得）のような項目コードが列名になっています。

使った変数は次のとおりです。

変数	定義（項目コード）	年
転入超過	（転入者数 A5103 − 転出者数 A5104）÷ 住民基本台帳人口 A2301 × 1,000	2024
納税者1人あたり所得（対数）	課税対象所得 C120110 ÷ 納税義務者数 C120120	2024年度
そのほかの条件	65歳以上の割合・大卒の割合・人口密度・第1次／第2次産業の割合・単独世帯の割合	2020

読み込みと変数づくりは、ほかの記事と共通の common.py の load() にまとめています。比率を作ったあと、条件に使う変数を全部 z 値（平均0・標準偏差1）にしているのが大事なところです。

LABEL = {"log_inc": "納税者1人あたり所得（対数）", "old": "65歳以上の割合", "univ": "大学・大学院卒の割合",
         "log_density": "人口密度（可住地、対数）", "mfg": "第2次産業で働く人の割合", "primary": "第1次産業で働く人の割合",
         "single": "単独世帯の割合"}
CONTROLS = ["log_inc", "old", "univ", "log_density", "mfg", "primary", "single"]


def load():
    d = pd.read_csv(ROOT / "data" / "municipal_base.csv")
    d = d[(d["A1101"] > 0) & (d["C120120"] > 0) & (d["A2301"] > 0)].copy()
    d["inc"] = d["C120110"] / d["C120120"] * 100                 # 万円（2024年度＝2023年の所得）
    d["log_inc"] = np.log(d["inc"])
    d["old"] = d["A1303"] / d["A1101"]
    d["univ"] = d["E9106"] / d["E9101"]
    d["mfg"] = d["F2211"] / d["F1102"]
    d["primary"] = d["F2201"] / d["F1102"]
    d["single"] = d["A810105"] / d["A710101"]
    d["ward23"] = d["code"].between(13101, 13123)
    for c in CONTROLS:
        d["z_" + c] = (d[c] - d[c].mean()) / d[c].std()
    return d

標準化しておくと、単位の違う変数（割合・対数・人数）の係数を「1標準偏差あたり、転入超過が1,000人あたり何人違うか」という同じ物差しで並べられます。

被説明変数の転入超過は、分析スクリプト側で1行で作っています。

d = load()
d["net_mig"] = (d["A5103"] - d["A5104"]) / d["A2301"] * 1000        # 住民1,000人あたりの転入超過（2024年）

全国1,740市区町村のうち、転入が転出を上回ったのは26%（share_positive 0.257）で、中央値は1,000人あたり −3.7人でした。

手法の実装：条件を1つずつ足す4本の回帰

「条件をそろえると答えが変わる」を見せるために、同じ係数（ここでは所得）を、条件を足しながら4回推定します。

そのまま（何もそろえない）
＋都道府県（C(pref_code) で都道府県の固定効果）
＋年齢構成（65歳以上の割合）
＋ほかの条件すべて

これを関数 steps() にしています。statsmodels の式（formula）を文字列で組み立てて、cov_type="HC3" で不均一分散に頑健な標準誤差にしています。市区町村は人口の規模がばらばらで、残差のばらつきがそろっている前提を置きにくいためです。

def steps(d, y, key, label_key, extra=None):
    """key の係数が、条件を足すとどう変わるか（標準化した1標準偏差あたり）"""
    others = [c for c in CONTROLS if c != key]
    specs = [("そのまま（何もそろえない）", f"{y} ~ z_{key}"),
             ("＋都道府県", f"{y} ~ z_{key} + C(pref_code)"),
             ("＋年齢構成", f"{y} ~ z_{key} + z_old + C(pref_code)" if key != "old" else f"{y} ~ z_{key} + C(pref_code)"),
             ("＋ほかの条件すべて", f"{y} ~ z_{key} + " + " + ".join("z_" + c for c in others) + " + C(pref_code)")]
    rows = []
    for name, f in specs:
        r = smf.ols(f, d).fit(cov_type="HC3")
        k = "z_" + key
        rows.append({"model": name, "coef": r.params[k], "lo": r.conf_int().loc[k, 0], "hi": r.conf_int().loc[k, 1],
                     "p": r.pvalues[k], "r2": r.rsquared, "n": int(r.nobs)})
    return pd.DataFrame(rows)

この関数は「どの変数の係数を追うか（key）」と「被説明変数（y）」を引数にしているので、同じリポジトリのほかの記事（出生と所得、昼夜間人口と所得など）でもそのまま使い回しています。

全部の条件を同時に入れたモデルの係数表は full_model() で作ります。

def full_model(d, y):
    f = f"{y} ~ " + " + ".join("z_" + c for c in CONTROLS) + " + C(pref_code)"
    r = smf.ols(f, d).fit(cov_type="HC3")
    out = pd.DataFrame({"var": [LABEL[c] for c in CONTROLS], "coef": [r.params["z_" + c] for c in CONTROLS],
                        "lo": [r.conf_int().loc["z_" + c, 0] for c in CONTROLS], "hi": [r.conf_int().loc["z_" + c, 1] for c in CONTROLS],
                        "p": [r.pvalues["z_" + c] for c in CONTROLS]})
    out["r2"], out["n"] = r.rsquared, int(r.nobs)
    return out

呼び出し側はこれだけです。

t = steps(d, "net_mig", "log_inc", "所得")
t.round(3).to_csv(RES / "a3_steps.csv", index=False, encoding="utf-8-sig")
fm = full_model(d, "net_mig"); fm.round(3).to_csv(RES / "a3_full.csv", index=False, encoding="utf-8-sig")

結果（`results/a3_steps.csv`）

モデル	所得の係数	95%区間	決定係数
そのまま（何もそろえない）	+2.086	1.518〜2.653	0.050
＋都道府県	+1.904	1.143〜2.664	0.146
＋年齢構成	+0.263	−0.646〜1.171	0.181
＋ほかの条件すべて	−0.577	−1.734〜0.579	0.234

都道府県をそろえても係数はほとんど変わらず、年齢構成を足したところで一気に0付近まで落ちています。「同じ県の中で比べても関係はあるが、若い人が多いかどうかで説明できてしまう」という形です。

結果の出し方：散布図と係数の階段図

散布図（対数目盛・人口で円の大きさ）

所得は右に長い分布なので、横軸を対数目盛にしています。目盛を「250, 300, 400…」と読める数字にするため、set_xticks でラベルを指定し、補助目盛のラベルは NullFormatter で消しています。東京23区だけ色を変えています。

fig, ax = plt.subplots(figsize=(8, 5))
o = d[~d.ward23]; w = d[d.ward23]
ax.scatter(o.inc, o.net_mig, s=np.clip(o.A2301 / 5000, 4, 120), color=GRAY, alpha=0.5, lw=0)
ax.scatter(w.inc, w.net_mig, s=np.clip(w.A2301 / 5000, 4, 120), color=BLUE, alpha=0.8, lw=0, label="東京23区")
ax.axhline(0, color=INK, lw=0.8); ax.set_xscale("log")
ax.set_xticks([250, 300, 400, 500, 700, 1000, 1500], ["250", "300", "400", "500", "700", "1,000", "1,500"])
from matplotlib.ticker import NullFormatter; ax.xaxis.set_minor_formatter(NullFormatter())
ax.set_ylim(-40, 40); ax.set_xlabel("納税者1人あたり所得（万円、対数目盛）"); ax.set_ylabel("転入超過（住民1,000人あたり、2024年）")

相関は0.22です。

係数の階段図

steps() の表を、そのまま横向きの区間つき点の図にします。95%区間が0をまたいだ段（p ≥ 0.05）は灰色にして、「はっきりしない」と書き添えます。数字のラベルが区間の右端に重ならないよう、x 軸の右側を広めに取っています。

def plot_steps(t, title, xlabel, path, color=BLUE):
    fig, ax = plt.subplots(figsize=(8, 3.6))
    y = np.arange(len(t))[::-1]
    for yi, (_, r) in zip(y, t.iterrows()):
        c = color if r.p < 0.05 else GRAY
        ax.hlines(yi, r.lo, r.hi, color=c, lw=2.5); ax.scatter(r.coef, yi, color=c, zorder=3, s=40)
        ax.text(max(r.hi, r.coef) + (t.hi.max() - t.lo.min()) * 0.03, yi, f"{r.coef:+.2f}" + ("" if r.p < 0.05 else "（はっきりしない）"), va="center", fontsize=9)
    ax.axvline(0, color=INK, lw=0.8); ax.set_yticks(y, t.model); ax.set_xlabel(xlabel)
    ax.set_xlim(min(t.lo.min(), 0) - (t.hi.max() - t.lo.min()) * 0.1, t.hi.max() + (t.hi.max() - t.lo.min()) * 0.45)
    ax.set_title(title, loc="left")
    fig.text(0.01, -0.06, SOURCE + "／重回帰・HC3 頑健標準誤差。灰色は95%区間が0をまたぐ", fontsize=8, color="#8a8984")
    fig.savefig(path); plt.close(fig)

全部入りモデルの係数図

full_model() の表は、係数の大きさ順に並べ、プラスで有意なら青、マイナスで有意ならオレンジ、0をまたげば灰色にしています。

def plot_coefs(t, title, xlabel, path):
    t = t.sort_values("coef")
    fig, ax = plt.subplots(figsize=(8, 4.2)); y = np.arange(len(t))
    col = [(BLUE if v > 0 else ORANGE) if p < 0.05 else GRAY for v, p in zip(t.coef, t.p)]
    ax.hlines(y, t.lo, t.hi, color=col, lw=2.2); ax.scatter(t.coef, y, color=col, zorder=3)
    ax.axvline(0, color=INK, lw=0.8); ax.set_yticks(y, t["var"]); ax.set_xlabel(xlabel); ax.set_title(title, loc="left")
    fig.text(0.01, -0.05, SOURCE + "／都道府県固定効果・HC3。灰色は95%区間が0をまたぐ", fontsize=8, color="#8a8984")
    fig.savefig(path); plt.close(fig)

はっきり残ったのは、大学・大学院卒の割合（+2.364）、人口密度（+2.743）、単独世帯の割合（−1.458）でした。決定係数は0.234です。

つまずいた点と落とし穴

係数の物差しをそろえる：標準化しないと、割合（0〜1）と対数の係数が並んでも大きさを比べられません。load() の最後で z 値の列（z_ 付き）を作っておき、式では必ず z_ の列を使うようにしました
年のずれ：所得は2024年度（2023年の所得）、転入超過は2024年、ほかの条件は主に2020年の国勢調査です。年がそろっていないことは、結果を読む前提として残ります
分母の違い：転入超過は住民基本台帳人口（A2301）で割り、65歳以上の割合などは国勢調査の人口（A1101）で割っています。どの分母を使ったかは README の表に書いて、後から追えるようにしました
散布図の外れ値：転入超過は人口の少ない町村で大きく振れるので、縦軸を ±40 で切り、そのことを図の注に書いています
「消えた」の読み方：係数が0と区別できなくなっても、「関係がない」の証明ではありません。「年齢構成などで説明できる」という意味です。また観察データの横断分析なので、因果（所得が高いから人が来る）は言えません

参考

note の記事：人が集まっている街は、所得が高い街か？全国1,740市区町村で「条件をそろえる」と答えが変わった
コード・データ・図：GitHub articles/008_migration_income

出典：政府統計の総合窓口(e-Stat)（https://www.e-stat.go.jp/）「社会・人口統計体系市区町村データ」「統計でみる市区町村のすがた2026」（総務省統計局）を加工して作成