LaLondeデータセットに対し傾向スコア（IPW）をPythonで実装して因果推論

はじめに

最近（2022年）のデータサイエンス周りでは、因果推論がホットな話題の一つになっています。私は現時点で高度な因果推論を実務で行うことはあまりないのですが、今後施策の効果検証などを高度な統計モデル・機械学習モデルを使って実施する場面がありそうなので勉強してみました。
因果推論でよく使われるデータセットの一つであるLaLondeデータセットを使って、Pythonで実際に因果推論をやってみました。記事が長くなるので、今回は使うデータの説明と傾向スコアを使った逆確率重み付き推定を実施します。
Doubly Robust 推定というちょっと高度な因果推論は別の記事にしたいと思います。

この記事で想定する読者

この記事では、以下の読者を主なターゲットして書かれてます。

• 因果推論という言葉を聞いたことはある
• 傾向スコア、そのほかの高度な因果推論で手を動かしたことはない
  因果推論を全く知らない方、あるいはECounMLやCausalMLをある程度使ったことがある方は対象としていませんので、ご了承ください。

就労支援の因果推論 - NSWデータ -

この記事ではデータとして一貫して就労支援に関するデータを扱います。労働市場に参加できない人々に対して、カウンセリングや就労経験を与えるプログラムを実施した場合、実際にどの程度年収増加をもたらすのか測定したいとします。

この効果を測定する場合、もし可能ならばランダムサンプリングにより就労支援を実施する介入群、および就労支援を実施しない対照群を作って、この２つの群に属する人々の収入が後からどの程度増加したかを比較するのがよいでしょう。こうしてランダム化することで、２つの群の年齢、教育年数などの分布が等しくなっていれば、収入の差は就労支援という介入によって生じたと結論付けられるでしょう。

実際に、このような実験が1970年代のアメリカで行われました。National Support Work(NSW)というプログラムです。
条件を満たす希望者から、ランダムに選択して就労支援を実施し、その後1978年の収入をデータとして記録しています。

• 結果変数 $Y$ (re78)： '78の収入
• 介入ダミー変数 $Z$ (treat) ： 1:援助を受ける, 0:受けない
• 共変量　：　 年齢(age)、教育年数(education)、黒人(black)、ヒスパニック(hispanic)など

ランダムサンプリングですので、年齢や教育年数の分布に大きな差はないと期待されます。実際に比較したのが下図になります。

こうして比較すると、分布の差がないとは言えないと思います。ただし、後述するLaLondeデータセットと比較すると差は小さい方なのです。ここではランダムサンプリングが適切に実施されているとして、就労支援の効果を確認します。
以下の線形モデルで効果を推定します。

$ Y = \beta_1 \times Z + \beta_0$
('78の収入 = $\beta_1 \times$ 介入ダミー変数 + $\beta_0$）

NSWのデータを読み込んで、線形モデルの評価を実施するPythonコードです。線形モデルを作るだけならscikit-learnでも実施できますが、推計された係数の有意性も確認したいのでstatsmodelsを使います。

import statsmodels.api as sm
import pandas as pd

nsw_data = pd.read_stata("https://users.nber.org/~rdehejia/data/nsw_dw.dta")

nsw_X0 = nsw_data[['treat']]
nsw_y = nsw_data['re78']

nsw_X0 = sm.add_constant(nsw_X0)

model0 = sm.OLS(nsw_y, nsw_X0)
result0 = model0.fit() 

result0.summary()

最後の行でsummary()を実行し、推定結果を表示します。

就労支援を受けた効果で、収入が1,794ドル増加していると推定されるます。係数も有意な結果となっています。

LaLondeデータセット　－就労支援に関する”観察データ”－

NSWのデータでは、就労支援という介入を実施する人をランダムに選んでいて、実験となっていました。しかし、実際にはランダムに選ぶことはできない状況もあり得ます。このような場合でも、因果推論の手法を使って適切に就労支援の効果を推定できるでしょうか？
LaLondeデータセットはこのような状況を、NSWデータともう一つ別の調査データであるCPSデータを使って作り出しています。

• NSWのうち、就労支援を実施した人だけを抽出し、介入群のデータとします。
• 別の調査データであるCurrent Population Survey(CPS)を対照群とします。

先ほどと同じように年齢、教育年数の分布を確認してみます。

上のNSWデータだけの場合と比較して、介入群（NSW）と対照群（CPS）の差異が大きくなっています。NSWの方が年齢は若く、教育年数は短くなっています。若くて教育年数が短ければ、就労支援の有無とは関係なしに収入は低くなるというバイアスが存在しそうで、正しく効果を推定できないのではないかと想像できます。

線形モデルで効果を推定します。上記の通り、年齢や教育年数などの共変量の分布の違いを考慮する必要があります。このため、今回は線形モデルの式に年齢や教育年数なども入れて推定します。

cps1_data = pd.read_stata("https://users.nber.org/~rdehejia/data/cps_controls.dta")
nsw_cps1_data = pd.concat( [nsw_data.query("treat==1"), cps1_data], axis=0 )

nsw_cps1_X0 = nsw_cps1_data[['treat', 'age', 'education', 'black', 'hispanic', 'married', 'nodegree', 're74', 're75']]

nsw_cps1_y = nsw_cps1_data['re78']

nsw_cps1_X0 = sm.add_constant(nsw_cps1_X0)

model_nc0 = sm.OLS(nsw_cps1_y, nsw_cps1_X0)
result_nc0 = model_nc0.fit()

result_nc0.summary()

結果は以下の通りとなります

就労支援を受けた効果で、収入が699ドル増加していると推定されます。（実際は1,794ドル）
係数も有意ではなく、就労支援の効果があったとは結論付けられなくなります。

傾向スコアとIPW

介入群と対照群の間で、年齢、教育年数などの共変量の分布に差異があったため、効果を正しく推定できませんでした。この差異を調整することができれば観察データでも効果量を正しく推定できるのではないか。この調整する方法として、傾向スコアがあります。因果推論の世界では有名な手法ですので、詳細は参考文献を見ていただきたいです。ここでは簡単に定義などを記します。

傾向スコアは、下式のように共変量を与えた時に介入群に割り当てられる確率として定義されます。
$e = P(Z=1|age, education, black, ...)$
傾向スコアを推定するときは、ロジスティック回帰モデルがよく使われます。

傾向スコアを用いて、共変量の分布の差異を調整する方法としては

• 傾向スコアマッチング
• 逆確率重み付き推定（IPW）

があります。今回は後者のIPWを使います。

傾向スコアを正しくモデル化できて、共変量が与えられた時に介入群に割り当てられる確率が正しく求められると、介入効果は以下の推定量で求めることができます。

$\frac{1}{N} \Sigma_i \frac{Z_i y_{i}}{e_i} - \frac{1}{N} \Sigma_i \frac{(1-Z_i) y_{i}}{1-e_i} $

ここで、

• $N$ : データの標本数
• $Z_i$ : $i$番目の人の介入ダミー変数。就労支援を受けていれば1、受けていなければ0
• $y_i$ : $i$番目の人の'78における収入
• $e_i$ : $i$番目の人の推定された傾向スコア

それでは、実際にPythonで傾向スコアを求めてみます。ロジスティック回帰を実行するために、statsmodelsの一般線形モデルを扱うライブラリであるGLMを使います。

X_IPW = nsw_cps1_data[['age', 'education', 'black', 'hispanic', 'married', 'nodegree', 're74', 're75']]
y_IPW = nsw_cps1_data.treat

X_IPW = sm.add_constant(X_IPW)

logit = sm.GLM(y_IPW, X_IPW, family=sm.families.Binomial())

logit_result = logit.fit()

logit_result.summary()

求めた傾向スコアのモデルを使って、各サンプルの傾向スコア値を求め、IPWを計算します。そして、IPWを重みとして線形回帰モデルを評価して、就労支援の介入効果を推定します。重み付きの線形回帰を実行するために、statsmodelsのWLSを使います。

nsw_cps1_data['PS_Score'] = logit_result.predict(X_IPW)
nsw_cps1_data["weight"] = nsw_cps1_data.apply( lambda x : 1/x.PS_Score if x.treat==1 else 1/(1-x.PS_Score), axis=1 )

X1 = nsw_cps1_data[['treat']]
y1 = nsw_cps1_data[['re78']]

X1 = sm.add_constant(X1)

lmdl = sm.WLS(y1, X1, weights=nsw_cps1_data['weight'])
lresult = lmdl.fit()

lresult.summary()

結果です。

介入の効果で、収入が8497ドル減少するという逆の結果が出てしまいました。しかも、係数は有意になっています。NSWでの結果を知らなければ、就労支援は逆効果であると結論付けてしまうところです。

なぜこのようなことが起きたのでしょうか。求めた傾向スコアの分布を介入群と対照群ごとに確認してみます。

青色の対照群(treat 0.0)は、傾向スコアが0付近に集中しています。CPSデータの方がデータ件数が多い不均衡データとなっており、CPSデータから構成される対照群の傾向スコアが総じて低く評価されています。0に近い値で割って重みを求めていますので、ちょっとした誤差で重みが大きく変わってしまいます。このため、傾向スコアを使っても介入効果を正しく推定できていないと考えられます。

最後に

今回は就労支援の介入効果検証でよく用いられるLaLondeデータセットを説明し、このデータに対して線形モデルとIPWをPythonで実装して効果を推定してみました。この方法では介入群と対照群の年齢、教育年数など共変量の分布差異をうまく調整することができず、適切に評価することはできませんでした。
因果推論では、Double Robust推定という手法や、最近の機械学習モデルを取り込む手法が出ていて、Pythonでもそれほど難しくなく実装できるようになっています。次回の記事ではこのテーマで記事を描きたいと思います。

参考文献

安井翔太著　効果検証入門
岩波データサイエンスVol.3