何を確かめたいか
県境をはさんで接する市区町村どうしで、最低賃金がより大きく上がった側の飲食・小売の従業者の伸びは小さかったかを、2016→2021年の経済センサスで調べました。この記事では、その分析を Python でどう実装したかに絞って書きます。結果の読み方は note の記事にまとめています。
「境界不連続デザイン」は、制度が切り替わる境界のすぐ両側を比べる方法です。隣どうしなら景気や気候はほぼ同じなので、残った違いを制度の差と考えやすくなります。
コードはすべて GitHub の articles/F15_border_minwage/ にあり、ここに載せるのはそこから抜き出したものです。
全体の流れ
pip install -r requirements.txt
# data/raw/ に SSDSE-A-2021.csv・SSDSE-A-2026.csv・minimum_wage_2002_2026.csv を置く（同梱）
cd src
python 01_pairs.py # 県境のペア（TopoJSON を取得）
python 02_census.py # 市区町村×産業の従業者数と最低賃金
python 03_reg.py # 推定・頑健性・並べ替え検定
python 04_図.py # 図5枚
| ファイル | やること |
|---|---|
01_pairs.py |
境界データから「県境をはさんで接するペア」と「中心が20km以内のペア」を作る |
02_census.py |
SSDSE から市区町村×産業の従業者数を読み、最低賃金と合わせて変化率を出す |
03_reg.py |
ペアの中の差の回帰（クラスタ標準誤差）、頑健性、並べ替え検定 |
04_図.py |
地図・分布・散布図・産業別・頑健性の図 |
使ったライブラリは pandas・numpy・statsmodels・matplotlib です。
データの取り方
1. 境界データ（TopoJSON）を取る
市区町村の境界は、国土数値情報の行政区域（2021）をもとにした簡略版の TopoJSON（smartnews-smri/japan-topography、CC BY 4.0）を使いました。初回だけ取得して data/raw/ に置きます。
url = ("https://raw.githubusercontent.com/smartnews-smri/japan-topography/main/"
"data/municipality/topojson/s0010/N03-21_210101.json")
dst = RAW / "N03-21_210101.json"
if not dst.exists():
urllib.request.urlretrieve(url, dst)
topo = json.load(open(dst, encoding="utf-8"))
geoms = list(topo["objects"].values())[0]["geometries"]
sc, tr = topo["transform"]["scale"], topo["transform"]["translate"]
TopoJSON の座標は「弧（arc）」ごとに差分で符号化され、整数に量子化されています。中心の位置を出すために、累積和を取って transform で経緯度に戻します。
arcs_xy = []
for arc in topo["arcs"]:
x = y = 0
pts = []
for dx, dy in arc:
x += dx
y += dy
pts.append((x * sc[0] + tr[0], y * sc[1] + tr[1]))
arcs_xy.append(np.array(pts))
2. SSDSE から従業者数を読む
従業者数は、統計センターの SSDSE-A（教育用標準データセット）の2021年版（経済センサス2016年6月）と2026年版（同2021年6月）から取りました。SSDSE は cp932 で、1行目が項目コード、3行目までが見出しなので、1行目だけを列名にして4行目から読みます。
IND = {"total": "", "food": "47", "retail": "41", "life": "48", "manuf": "37", "medical": "50"}
def load(fn, year_tag):
h = pd.read_csv(RAW / fn, encoding="cp932", header=None, nrows=1).iloc[0].tolist()
d = pd.read_csv(RAW / fn, encoding="cp932", header=None, skiprows=3)
d.columns = h
d = d.rename(columns={h[0]: "code", h[1]: "pref", h[2]: "muni"})
d["code"] = d["code"].astype(str).str[1:].astype(int)
out = d[["code", "pref", "muni"]].copy()
out["pop"] = pd.to_numeric(d["A1101"], errors="coerce")
for k, s in IND.items():
out[f"emp_{k}"] = pd.to_numeric(d[f"C2208{s}"], errors="coerce")
out[f"est_{k}"] = pd.to_numeric(d[f"C2108{s}"], errors="coerce")
out.columns = ["code", "pref", "muni"] + [f"{c}_{year_tag}" for c in out.columns[3:]]
return out
C2208 が従業者数（民営）、C2108 が事業所数（民営）で、後ろに産業大分類の番号がつきます。市区町村コードは先頭に1文字ついているので、str[1:] で落としています。
3. 最低賃金を「調査日に効いていた額」で合わせる
経済センサスの調査日は6月1日で、そのとき効いている最低賃金は前年度10月の改定額です。そのため、2016年には2015年度の額、2021年には2020年度の額を合わせました。
m["pref_s"] = m["pref"].str.replace(r"(都|府|県)$", "", regex=True)
m.loc[m.pref == "北海道", "pref_s"] = "北海道"
m = m.merge(mw[["pref_s", "2015", "2020"]].rename(columns={"2015": "mw15", "2020": "mw20"}), on="pref_s", how="left")
assert m.mw15.notna().all(), m[m.mw15.isna()].pref.unique()
m["dlog_mw"] = np.log(m.mw20) - np.log(m.mw15)
for k in IND:
m[f"dlog_{k}"] = np.log(m[f"emp_{k}_21"]) - np.log(m[f"emp_{k}_16"])
m[f"dlog_est_{k}"] = np.log(m[f"est_{k}_21"]) - np.log(m[f"est_{k}_16"])
m = m.replace([np.inf, -np.inf], np.nan)
県名の結合で1県でも落ちると、その県の市区町村がまるごと抜けます。assert で止めるようにしています。
ペアを作る：「同じ弧を共有する」で隣どうしを決める
隣どうしかを判定するのに、ポリゴンの交差判定は使いませんでした。TopoJSON では隣どうしの市区町村が同じ弧を共有するので、「同じ弧を使っている、県の違う2つの市区町村」を接しているとみなせます。簡略化した境界でも隙間ができる心配がありません。
arc_owner = defaultdict(set)
pts_by = defaultdict(list)
for g in geoms:
p = g.get("properties", {})
c = p.get("N03_007")
if not c or g.get("arcs") is None:
continue
city = to_city(int(c))
for i in flat(g["arcs"]):
k = i if i >= 0 else ~i
arc_owner[k].add(city)
pts_by[city].append(arcs_xy[k])
pairs = set()
for k, s in arc_owner.items():
s = sorted(s)
for i in range(len(s)):
for j in range(i + 1, len(s)):
a, b = s[i], s[j]
if a // 1000 != b // 1000: # 県が違う
pairs.add((a, b))
ポイントは2つです。
- 弧の番号は、向きが逆のとき負の数（ビット反転
~i）で入っています。k = i if i >= 0 else ~iで元の番号に戻します - 市区町村コードの上2桁が都道府県なので、
a // 1000 != b // 1000で「県が違う」を判定できます
これで724組のペア（656市区町村、県の組み合わせ＝境界は87本）ができました。
頑健性チェック用に、接していなくても中心どうしが20km以内のペアも作りました。中心は境界の点の平均で近似し、距離は緯度で経度の長さを補正して出しています。
cent = {c: np.vstack(v).mean(axis=0) for c, v in pts_by.items()}
def dist_km(a, b):
(x1, y1), (x2, y2) = cent[a], cent[b]
kx = 111.32 * np.cos(np.radians((y1 + y2) / 2))
return float(np.hypot((x1 - x2) * kx, (y1 - y2) * 110.57))
手法の実装：ペアの中の差を回帰する
式は次のとおりです。Δ は2016→2021年の対数差です。
(Δlog 従業者_A − Δlog 従業者_B) = β ×（Δlog 最低賃金_A − Δlog 最低賃金_B）+ 誤差
ペアの中で差を取ると、隣どうしに共通の要因は消えます。ペア固定効果を入れた回帰と同じことです。β は弾力性（最低賃金が1%多く上がると従業者が何%変わるか）になります。
ペアごとの差を作る関数です。
def build(pairs, ycol, xcol="dlog_mw"):
a, b = m.loc[pairs.code_a], m.loc[pairs.code_b]
d = pd.DataFrame({
"dy": a[ycol].values - b[ycol].values,
"dx": a[xcol].values - b[xcol].values,
"seg": pairs.seg.values,
"w": np.sqrt(a["pop_21"].values * b["pop_21"].values),
"mwgap": np.abs(a["mw20"].values - b["mw20"].values),
}).replace([np.inf, -np.inf], np.nan).dropna()
return d
推定は statsmodels の WLS で、重みなしのときは全部1にしています。
def fit(d, weight=False, winsor=0.01):
d = d.copy()
if winsor:
lo, hi = d.dy.quantile([winsor, 1 - winsor])
d["dy"] = d.dy.clip(lo, hi)
mod = sm.WLS(d.dy, d[["dx"]], weights=d.w if weight else np.ones(len(d)))
r = mod.fit(cov_type="cluster", cov_kwds={"groups": pd.factorize(d.seg)[0]})
b, se = r.params["dx"], r.bse["dx"]
return {"beta": b, "se": se, "lo": b - 1.96 * se, "hi": b + 1.96 * se,
"p": r.pvalues["dx"], "n_pairs": len(d), "n_seg": d.seg.nunique()}
実装で気をつけたのは次の3点です。
-
定数項を入れない：A と B のどちらを先にするかに意味はないので、
sm.add_constantはしていません。定数項を入れると、ペアの並べ方で結果が変わってしまいます -
県の組でクラスタ：同じ境界のペアは同じ最低賃金の差を持ちます。
cov_type="cluster"に、境界の ID をpd.factorizeで整数にしたものを渡しています - 外れ値を上下1%で丸める：従業者が少ない町では、数人の増減で対数差が大きく振れます。丸めない場合も頑健性チェックで出しています
製造業を対照にした「三重差」は、同じ市区町村の中で飲食＋小売の伸びから製造業の伸びを引いた列を作って、同じ fit に渡すだけです。
m["dlog_fr"] = np.log((m.emp_food_21 + m.emp_retail_21) / (m.emp_food_16 + m.emp_retail_16))
m["dlog_fr_minus_manuf"] = m.dlog_fr - m.dlog_manuf
res["main_foodretail"] = fit(build(touch, "dlog_fr"))
res["main_fr_minus_manuf"] = fit(build(touch, "dlog_fr_minus_manuf"))
頑健性チェックを辞書にまとめる
設定を変えた推定を、名前つきで辞書に入れておくと、表と図にそのまま使えます。
rob = {}
Y = "dlog_fr"
rob["基本（接するペア）"] = res["main_foodretail"]
rob["中心が20km以内のペアも入れる"] = fit(build(P[(P.kind == "touch") | (P.dist_km <= 20)], Y))
rob["大都市の境を除く"] = fit(build(touch[~touch.metro], Y))
rob["人口で重みづけ"] = fit(build(touch, Y), weight=True)
d0 = build(touch, Y)
rob["最低賃金の差が大きいペアだけ（20円以上）"] = fit(d0[d0.mwgap >= 20])
rob["外れ値を丸めない"] = fit(d0, winsor=0)
rob["事業所数で見る（飲食）"] = fit(build(touch, "dlog_est_food"))
rob["製造業を引く（三重差）"] = res["main_fr_minus_manuf"]
並べ替え検定：境界ごとに符号を入れ替える
クラスタ数が87と多くないので、標準誤差だけに頼らず並べ替え検定もしました。ペアごとではなく「境界ごと」に最低賃金の差の符号をランダムに入れ替えて、偶然でどれくらいの傾きが出るかを2,000回見ています。
rng = np.random.default_rng(SEED)
d = build(touch, Y)
lo, hi = d.dy.quantile([0.01, 0.99]); d["dy"] = d.dy.clip(lo, hi)
obs = (d.dx * d.dy).sum() / (d.dx ** 2).sum()
segs = d.seg.unique()
perm = []
for _ in range(2000):
sgn = dict(zip(segs, rng.choice([-1, 1], len(segs))))
s = d.seg.map(sgn)
perm.append(((d.dx * s) * d.dy).sum() / (d.dx ** 2).sum())
perm = np.array(perm)
p_perm = float((np.abs(perm) >= abs(obs)).mean())
定数項なしの単回帰なので、傾きは Σ dx·dy / Σ dx² で直接計算できます。毎回 statsmodels を呼ぶより速く済みます。
結果の出し方
産業別の推定を1枚にする
results/estimates.json から読み、errorbar で点と95%区間を並べます。
r = pd.DataFrame(E["main"]).T.loc[list(lab)]
fig, ax = plt.subplots(figsize=(6.6, 4.2))
yy = np.arange(len(r))[::-1]
ax.errorbar(r.beta, yy, xerr=[r.beta - r.lo, r.hi - r.beta], fmt="o", color="#4C72B0", capsize=3)
ax.axvline(0, color="k", lw=0.8)
ax.set_yticks(yy); ax.set_yticklabels([lab[i] for i in r.index])
飲食＋小売の弾力性は1.41（95%区間 −1.96〜+4.78）、製造業を引いた三重差は −1.12（−4.97〜+2.73）で、どちらも0をまたぎました。並べ替え検定の p 値は0.44です。
頑健性の図に「見つけたい大きさ」の帯を重ねる
頑健性の図には、海外の研究でよく報告される弾力性 0〜−0.3 の範囲を axvspan で重ねました。区間の広さと比べられるようにするためです。
R = pd.DataFrame(E["robust"]).T
fig, ax = plt.subplots(figsize=(6.8, 4.2))
yy = np.arange(len(R))[::-1]
ax.errorbar(R.beta, yy, xerr=[R.beta - R.lo, R.hi - R.beta], fmt="o", color="#DD8452", capsize=3)
ax.axvline(0, color="k", lw=0.8)
ax.axvspan(-0.3, 0, color="#55A868", alpha=0.25, label="0〜−0.3")
つまずいた点と落とし穴
政令市の区と SSDSE の単位がそろわない
境界データは政令市が区ごとに分かれていますが、SSDSE は市単位です。そのままだと区のコードが SSDSE に見つからず、ペアから落ちます。区のコードを「同じ県の中で、そのコードより小さい SSDSE のコードのうち最大のもの」に寄せて、市にまとめました。
def to_city(code: int) -> int:
if code in codes:
return code
cand = [c for c in codes if c < code and c // 1000 == code // 1000]
return max(cand) if cand else code
政令市のコードは区より小さい番号で、区はその直後に並ぶので、この方法で寄せられます。
弧のリストは入れ子になっている
geometry["arcs"] は Polygon と MultiPolygon で入れ子の深さが違います。深さを気にせず番号だけ取り出すために、再帰で平らにする関数を挟みました。
def flat(a):
if isinstance(a, int):
return [a]
out = []
for b in a:
out += flat(b)
return out
「水準の差」と「上がり方の差」は別物
最低賃金の差を最初に見たとき、水準の差（中央値33円、最大175円）は十分あるように見えました。ところが、変化で比べる式に入るのは「上がり方の差」で、2015→20年度の上昇率はどの県も11.1〜13.6%に収まり、ペアの中の差は最大でも1.8ポイントでした。
説明変数の散らばりが小さいと、区間は広くなります。この設計で80%の確率で見つけられる弾力性は約±4.8で、0〜−0.3 程度の効果は見分けられません。推定を回す前に説明変数の分布を図にしておくと、「差が見えなかったときにどう読むか」を先に決めておけます。
0を外れた結果に飛びつかない
20km以内のペアまで広げると、1.46（+0.34〜+2.57）と区間が0を外れました。ただし向きは「最低賃金が多く上がった側で雇用が伸びた」で、対照の製造業も同じ向きに動いています。地域全体の動きの違いを拾っていると考え、効果の証拠としては扱っていません。
言えないこと
- 「最低賃金は雇用に影響しない」とは言えません。ありそうな大きさの効果を見分けられない設計でした
- 経済センサスは2016年と2021年の2時点だけで、前の時期の動きがそろっていたかは確かめられません
- 2021年6月はコロナの最中で、県ごとの営業制限の違いが飲食の雇用に混ざっています
- 通勤で市区町村をまたぐ人の影響、個人の賃金・労働時間・パートの割合は分かりません
参考
- note の記事（結果の読み方）：https://note.com/takayukimatsuura/n/n82a98a9d71cd
- コードと集計結果：https://github.com/takayuki-0127/japan-open-data-notes/tree/main/articles/F15_border_minwage
- 出典：独立行政法人統計センター「SSDSE」SSDSE-A-2021・SSDSE-A-2026（元は総務省・経済産業省「経済センサス‐活動調査」、総務省「国勢調査」）、厚生労働省「地域別最低賃金の全国一覧」、国土交通省「国土数値情報（行政区域データ）」をもとにした smartnews-smri/japan-topography（CC BY 4.0）を加工して作成