何を確かめたいか
2010→2015年に保育所等に通う子どもが増えた市区町村ほど、女性の労働力率は上がったかを、全国1,611市区町村で調べました。この記事では、その分析を Python でどう実装したかに絞って書きます。結果の読み方は note の記事にまとめています。
使った方法は「ダブル機械学習（DML）」です。そろえたい条件（コントロール変数）が多いときに、結果と処置の両方を機械学習で予測し、予測からの外れ（残差）どうしで回帰して効果を推定します。ライブラリは DoubleML を使いました。
コードはすべて GitHub の articles/F16_childcare_dml/ にあり、ここに載せるのはそこから抜き出したものです。
全体の流れ
pip install -r requirements.txt
cd src
python 02_整形.py # data/raw → data/panel_long.csv
python 03_分析.py # results/results.json・estimates.csv（2コアで約30分）
python 04_図.py # figures/
| ファイル | やること |
|---|---|
01_取得.md |
e-Stat「社会・人口統計体系 市区町村データ」の取り方のメモ |
02_整形.py |
取得した表を「市区町村×年度」の縦長の表にする |
common.py |
結果・処置・80の条件を作る関数 |
03_分析.py |
単純な回帰・重回帰・DML（3種類の学習器）、頑健性、プラセボ |
04_図.py |
図6枚 |
requirements.txt は pandas・numpy・statsmodels・scikit-learn・lightgbm・doubleml・matplotlib です。
データの取り方
1. e-Stat の市区町村データを取る
データは e-Stat の地域統計（市区町村データ）の表示画面から取りました。API キーの登録は要りません。画面で一度「表示」したあと、そのとき送られる api_get_result（POST）のパラメータのうち cols（項目の一覧）と tops（年度）だけを書き換えて送り直しています。どちらも JSON を gzip して base64 にした文字列で、inputNumberOfRows を 2000 にすると全市区町村が1回で返ります（詳しくは src/01_取得.md）。
項目は1回に25項目までなので、3つのまとまり（b1〜b3）に分けて取り、{"2015b1": [[見出し…], [地域, 値…], …], …} の形の JSON に保存しました。
2. 「市区町村×年度」の表にする
02_整形.py では、まとまりごとの表を年度ごとに横に結合し、最後に縦に積みます。
def tidy(rows, y):
h = rows[0]; out = []
for r in rows[1:]:
code, name = r[0].split('_', 1)
rec = {'code': code, 'name': name, 'year': y}
for c, v in zip(h[1:], r[1:]):
k = c.split('_')[0]; v = v.replace(',', '')
try: rec[k] = float(v)
except ValueError: rec[k] = np.nan
out.append(rec)
return pd.DataFrame(out)
def is_ward(code):
# 政令指定都市の区（例 01101）は除き、市全体（例 01100）を残す。東京23区は残す
return code[2] == '1' and code[4] != '0' and code[:2] != '13'
値は「1,234」のようにカンマ入りの文字列で返ってくるので、カンマを外して float にし、変換できないもの（「-」など）は欠損にしています。政令市は市全体と区の両方が入っているため、区を落として市全体だけを残します。東京23区は区ごとに残します。
3. 結果・処置・条件を作る
結果（女性の労働力率）と処置（保育所等に通う子ども）は common.py の関数で作ります。
def lfpr(g):
"""女性の労働力率（%）＝労働力人口（女）÷（労働力人口（女）＋非労働力人口（女））"""
return 100 * g.F110102 / (g.F110102 + g.F110802)
def enroll(g, kodomoen=False):
"""保育所等に通う子ども（0〜5歳100人あたり）"""
n = g.J2506.copy()
if kodomoen:
n = n + g.J2545.fillna(0) + g.J2546.fillna(0)
return 100 * n / g.A1405
80の条件は、2010年の地域の性格（人口・密度・年齢構成・世帯の形・産業・通勤・所得・財政力・学歴・雇用の形など）、2005→10年の変化（前の傾向）、同じ期間の年齢構成の変化、都道府県のダミー（46列）です。一部を抜き出します。
X['log人口'] = np.log(g.A1101)
X['log人口密度'] = np.log(g.A1101 / g.B1103)
X['女性65歳以上の割合'] = 100 * g.A130302 / g.A110102
# …（中略）…
X['前の5年の労働力率の変化'] = lfpr(g) - lfpr(g5)
X['前の5年の在所児の変化'] = enroll(g) - enroll(g5)
# 財政力指数は東京23区などで空欄。中央値で埋め、空欄だったことを示す列を足す
X['財政力指数が空欄'] = X['財政力指数'].isna().astype(float)
X['財政力指数'] = X['財政力指数'].fillna(X['財政力指数'].median())
4. 推定に使う表をまとめる
03_分析.py の build で、差（2015年−2010年）をとった Y・D と条件を1つの表にします。効果の単位は「在所児が100人あたり10人増えたとき」にしたいので、D を10で割っています。
UNIT = 10 # 効果は「0〜5歳100人あたりの在所児が10人増えたとき」で表す
MIN_CHILD = 100 # 0〜5歳が100人未満の町村は、割合が大きく振れるので除く
def build(a=2010, b=2015, y_fn=lfpr, d_a=None, d_b=None, x_pre=True, extra_drop=None):
X = covariates(w)
# …（中略）…
Y = y_fn(w[b]) - y_fn(w[a])
D = ((d_b if d_b is not None else enroll(w[b])) - (d_a if d_a is not None else enroll(w[a]))) / UNIT
df = pd.concat([Y.rename('Y'), D.rename('D'), X, C, P], axis=1)
df = df[(w[2010].A1405 >= MIN_CHILD).reindex(df.index).fillna(False)]
df = df.replace([np.inf, -np.inf], np.nan).dropna()
これで 1,611 市区町村・80列の表になります。
手法の実装
DML の本体
DoubleML の DoubleMLPLR（部分線形モデル）を使います。ml_l が Y を X から予測するモデル、ml_m が D を X から予測するモデルです。score='partialling out' にすると、「Y の残差を D の残差で回帰する」という形の推定になります。
def _dml_fit(df, kind='rf', folds=5, reps=5, seed=SEED):
xcols = [c for c in df.columns if c not in ('Y', 'D')]
data = dml.DoubleMLData(df, y_col='Y', d_cols='D', x_cols=xcols)
ml_l, ml_m = learners(kind)
np.random.seed(seed)
m = dml.DoubleMLPLR(data, ml_l=ml_l, ml_m=ml_m, n_folds=folds, n_rep=reps, score='partialling out')
m.fit(store_predictions=True)
ci = m.confint().iloc[0]
out = dict(coef=float(m.coef[0]), se=float(m.se[0]), lo=float(ci.iloc[0]), hi=float(ci.iloc[1]), p=float(m.pval[0]), n=len(df),
reps=[float(v) for v in np.ravel(m.all_coef)])
# 予測の当てはまり（それぞれの部分の R²）
pl = m.predictions['ml_l'][:, 0, 0]; pm = m.predictions['ml_m'][:, 0, 0]
out['r2_y'] = float(1 - np.mean((df.Y - pl) ** 2) / np.var(df.Y))
out['r2_d'] = float(1 - np.mean((df.D - pm) ** 2) / np.var(df.D))
return out
n_folds=5 が交差適合（データを5つに分け、4つで学習して残り1つの残差を出す）、n_rep=5 がその分け方のくり返し回数です。くり返した推定値は all_coef に入り、coef はその中央値になります。store_predictions=True にしておくと、交差適合での予測が残るので、2つの予測の当てはまり（R²）も出せます。
学習器を3種類用意する
学習器を変えても答えが大きく動かないかを見るため、ランダムフォレスト・LightGBM・Lasso を用意しました。
def learners(kind):
if kind == 'rf':
m = RandomForestRegressor(n_estimators=200, min_samples_leaf=5, max_features=0.33, n_jobs=-1, random_state=SEED)
return m, m
if kind == 'lgbm':
m = LGBMRegressor(n_estimators=300, learning_rate=0.04, num_leaves=15, min_child_samples=20, subsample=0.8, subsample_freq=1, colsample_bytree=0.7, verbose=-1, random_state=SEED)
return m, m
if kind == 'lasso':
m = make_pipeline(StandardScaler(), LassoCV(cv=5, random_state=SEED, max_iter=20000))
return m, m
Lasso は変数の大きさに左右されるので、StandardScaler と組み合わせています。
比べる相手：単純な回帰と重回帰
DML の値だけでは「80の条件をそろえたことで何が変わったか」が分からないので、statsmodels の OLS（不均一分散に強い HC3 の標準誤差）で、何もそろえない回帰と、同じ80列を直線で入れた重回帰も出します。
def ols(df, controls=True):
cols = ['D'] + ([c for c in df.columns if c not in ('Y', 'D')] if controls else [])
r = sm.OLS(df.Y, sm.add_constant(df[cols])).fit(cov_type='HC3')
lo, hi = r.conf_int().loc['D']
return dict(coef=float(r.params['D']), se=float(r.bse['D']), lo=float(lo), hi=float(hi), p=float(r.pvalues['D']), n=len(df))
頑健性とプラセボ
同じ dml_fit に渡す表や引数を変えるだけで、頑健性の推定がまとめて書けます。
rob['分割 2'] = dml_fit(df, 'rf', folds=2, key='r_分割 2')
rob['乱数を変える'] = dml_fit(df, 'rf', seed=SEED + 99, key='r_乱数を変える')
rob['人口1万人未満を除く'] = dml_fit(df[np.exp(df['log人口']) >= 10000], 'rf', key='r_人口1万人未満を除く')
rob['前の傾向を入れない'] = dml_fit(df.drop(columns=['前の5年の労働力率の変化', '前の5年の在所児の変化']), 'rf', key='r_前の傾向を入れない')
プラセボは、Y を「前の5年（2005→10年）の労働力率の変化」に差し替えます。保育所が増える前の変化まで後の在所児の変化で説明できてしまうなら、「もともと伸びていた町で保育所が増えた」だけかもしれない、という確認です。前の傾向の列と、同じ期間の年齢構成の変化の列は外しておきます。
dp = build(x_pre=False)
dp['Y'] = (lfpr(w[2010]) - lfpr(w[2005])).reindex(dp.index)
dp = dp.drop(columns=['女性65歳以上の割合の変化', '女性15〜64歳の割合の変化', '0〜5歳の割合の変化', '人口増減率']).dropna()
R['placebo'] = dml_fit(dp, 'rf', key='placebo')
結果の出し方
推定値は results/results.json にすべて残し、主なものを results/estimates.csv に表にしています。
| 推定 | 係数 | 95%信頼区間 |
|---|---|---|
| 単純な回帰 | +0.048 | −0.021〜+0.117 |
| 重回帰（80列を直線で） | +0.048 | −0.012〜+0.108 |
| DML（Lasso） | +0.047 | −0.010〜+0.103 |
| DML（ランダムフォレスト） | +0.061 | −0.002〜+0.125 |
| DML（LightGBM） | +0.066 | +0.006〜+0.125 |
単位は「在所児が100人あたり10人増えたときの、女性の労働力率の変化（ポイント）」です。プラスだが小さく、区間が0をまたぐものも多い、という結果でした。
推定値の比較図は、DML を青、比べる相手を灰色にして、errorbar で95%信頼区間を描いています（04_図.py）。
交差適合の分け方を20通り変えた推定値（n_rep=20）は、all_coef をそのままヒストグラムにしました。範囲は +0.054〜+0.081 で、すべてプラスでした。
頑健性とプラセボは、1枚に横向きで並べています。プラセボはオレンジにしました。プラセボの推定は −0.037（−0.119〜+0.047）で、0と区別できませんでした。
つまずいた点と落とし穴
1. 1回の実行が長いので、推定ごとに結果を保存する
DML はランダムフォレストを「5分割×5回×2つの予測」で学習するので、頑健性まで含めると全体で約30分かかります（2コア）。途中で止まってもやり直さずに済むよう、推定ごとに結果を JSON に書き出し、同じキーがあれば読み込むだけにしました。
CACHE_F = RES / '_cache.json'
CACHE = json.load(open(CACHE_F)) if CACHE_F.exists() else {}
def dml_fit(df, kind='rf', folds=5, reps=5, seed=SEED, key=None):
if key and key in CACHE:
return CACHE[key]
out = _dml_fit(df, kind, folds, reps, seed)
if key:
CACHE[key] = out; json.dump(CACHE, open(CACHE_F, 'w'), ensure_ascii=False)
print('done', key, round(out['coef'], 4), round(out['lo'], 4), round(out['hi'], 4), flush=True)
return out
設定を変えたときは、キャッシュを消してから流し直す必要があります。
2. 空欄を落とすと、東京23区がまるごと消える
財政力指数は東京23区などで空欄です。dropna() をそのまま使うと23区がすべて落ちるので、中央値で埋めたうえで「空欄だった」ことを示す列を足しました。
3. 年度によって、ない項目がある
保育所等在所児数の詳細票（J2506）は2017年度まで、0〜5歳人口（A1405）は国勢調査の年だけです。長い差（2010→20年）の確認では、2017年の在所児を2015年の0〜5歳人口で割っています。
g17 = w[2017].copy(); g17['A1405'] = w[2015].A1405 # 2017年の0〜5歳人口はないので2015年の値で割る
d3 = build(a=2010, b=2020, d_a=enroll(w[2010]), d_b=enroll(g17, kodomoen=True))
また、2015年は認定こども園への移行が始まった年で、J2506 に含まれない施設に移った分は「減った」ように見えることがあります。長い差では認定こども園（J2545・J2546）を足して確かめました。
4. 小さな町村では割合が大きく振れる
0〜5歳が数十人の町村では、在所児が数人増えるだけで「100人あたり」の値が大きく動きます。0〜5歳が100人未満の町村は除きました（MIN_CHILD = 100）。人口1万人未満を除く頑健性では推定値が +0.028 まで小さくなるので、小さな町の扱いで結果が動くことは意識しておく必要があります。
5. 重要度は、県のダミーをまとめて見る
「何が予測に効いたか」をランダムフォレストの重要度で見るとき、46列の県ダミーは1列ずつだと小さく見えます。合計して1行にまとめてから並べました。
s = pd.Series(rf.feature_importances_, index=xcols)
pref = s[[c for c in xcols if c.startswith('県_')]].sum()
s = s[[c for c in xcols if not c.startswith('県_')]]
s['都道府県（46のダミーの合計）'] = pref
Y の予測にいちばん効いたのは「女性65歳以上の割合の変化」（高齢化）、D の予測には「前の5年の在所児の変化」でした。なお、この重要度は全データで当てはめたもので、DML の推定そのものとは別の参考値です。
言えないこと
- 80の条件で取り除けない要因（共働きを望む若い世帯の転入、企業の求人など）があれば、推定は偏ります。DML でも、観察されていない交絡は消えません
- 結果は15歳以上の女性全体の労働力率で、子育て中の母親だけの値ではありません
- 処置は「定員」ではなく「在所児数」で、需要にも左右されます
参考
- note の記事（結果の読み方）：https://note.com/takayukimatsuura/n/ncbe1abba4999
- コードと集計結果：https://github.com/takayuki-0127/japan-open-data-notes/tree/main/articles/F16_childcare_dml
- 出典：政府統計の総合窓口(e-Stat)（https://www.e-stat.go.jp/）「社会・人口統計体系 市区町村データ」（総務省統計局）を加工して作成