傾向スコアと効果推定精度についての考察

記事の概要

　因果推論分野で登場する傾向スコアを勉強していて少し気になったことがあったので実験してみました。気になったことというのは「傾向スコアの算出」にどれだけ&どのように拘るべきかという点です。今回はわざとバイアスのあるデータを作り、

1. 真の傾向スコアを用いた場合
2. データセット生成と同じモデルで傾向スコアを算出した場合（ロジスティック回帰）
3. 介入の有無の予測精度を重視した場合（LightGBM）

の３通りで効果推定の精度を比較してみました。なお、この記事では「効果」といえば「ATE ( Average Treatment Effect )」を指します。

　初投稿ということでドキドキですが、何かまずい点がございましたらご教示くださると幸いです。

バイアスのあるデータセットの用意

　特徴量$x_0, x_1, x_2$に基づいてロジスティック回帰と同じモデルで介入$z \in \{0, 1 \}$を行う確率（＝真の傾向スコア）を決めます。$z$による真の効果（effect）はテキトーに1.3としており、効果の対象となるtargetは$x_0, x_1, x_2, z$および誤差項の線形結合で計算されるものとします。線形結合の係数もテキトーに決めています。

import numpy as np
import pandas as pd

# sample size
n = 10000

# features
np.random.seed(0)
x_0 = np.random.randn(n)

np.random.seed(1)
x_1 = np.random.randn(n)

np.random.seed(2)
x_2 = np.random.randn(n)

# treatment
z = np.zeros_like(x_0)
true_P_score = np.zeros_like(x_0)

# true effect of z
effect = 1.3

for i in range(n):
    # those who have high value of x_0, x_1 and x_2 tend to be z = 1.
    true_P_score[i] = 1 / (1 + np.exp(- x_0[i] - x_1[i] - x_2[i]))

    if np.random.rand(1) < true_P_score[i]:
        z[i] = 1
    else:
        z[i] = 0
        
# error
np.random.seed(3)
error = 0.01 * np.random.randn(n)

# generate target y
target = 2*x_0 + 3*x_1 + 1.5*x_2 + effect*z + error

# make dataframe
data = pd.DataFrame({'x_0': x_0,
                     'x_1': x_1,
                     'x_2': x_2,
                     'z': z,
                     'target': target,
                     'true_P_score': true_P_score})

　一応、無作為抽出を仮定した場合にどのくらいbiasの影響を受けてしまうのか確認してみると、見かけ上かなり大きく見積もってしまっていることがわかります。

真の効果: 1.3
見かけの効果: 5.5

# confirm the bias
print("the true effect of z = ", effect)
print('the pseudo effect of z = ',
      np.mean(data[data.z == 1].target) - np.mean(data[data.z == 0].target))

傾向スコアを用いたIPW

　傾向スコアを重みとしてIPW（逆確率重み付き推定）で効果を推定してみます。原理的には傾向スコアを$P$として、$z=1$のサンプル$i$には重み$Z_i/P_i$、$z=0$のサンプル$i$には重み$(1-Z_i)/(1-P_i)$を乗せて$z=1, 0$それぞれのtargetの期待値を計算し、その差が効果の推定量になります。

　しかしここではライブラリにｔ検定までやってほしいので少し工夫してやると、すべてのサンプルに対して重み$[ Z_i/P_i + (1-Z_i)/(1-P_i) ]$を乗せた、$z$による$target$の重み付き線形単回帰を考えてやればよいことがわかります。

※　説明のためにコードの順番が前後していますが、実際にはこのコードより前に傾向スコアを算出し、'Z/P', '(1-Z)/(1-P)'のカラムを用意しておく必要があります。

from statsmodels.regression.linear_model import WLS # weighted least square
from statsmodels.tools import add_constant

X = add_constant(data['z'])
y = data['target']
weights = data['Z/P'] + data['(1-Z)/(1-P)']
wls = WLS(y, X, weights=weights).fit()
wls.summary()

　蛇足かもしれませんが、以下のコードで真の効果と推定結果が出力できます。

print("the true effect of z = ", effect)
print("effect of z is estimated as", wls.params[1])

傾向スコアによる効果推定の比較

　さて、真の効果はデータセット生成時に決めたように1.3ですが、何通りかの傾向スコアを用いてその推定を行いました。傾向スコアというものが本当にデータセットに存在するのであれば、効果の推定量は不偏であるので、データセット生成時に利用した介入確率をそのまま傾向スコアとして用いた場合が最も精度が高くなるはずです。

　また、ロジスティック回帰のモデルは、データセット生成機構と同じものであるので、算出される傾向スコアも真の値に近くなるはずです。したがって効果の推定精度も結構良いのでは？と予想されます。ちなみに介入$z$の"予測精度"という意味では約76%でした。（いわゆる訓練データしかないので予測と呼ぶのかは謎です。）また、傾向スコアは真の値から平均的に$\pm 0.4$%程度（誤差の比率ではなく差の絶対値の平均）でした。

　最後にLightBGMで介入$z$の"予測精度"を向上させてみました。介入$z$の予測精度は約80%です。ハイパーパラメータはデフォルトのままです。傾向スコアは真の値から平均的に$\pm 6$%程度（誤差の比率ではなく差の絶対値の平均）でした。

これら傾向スコアを用いた効果推定結果が以下のようになっています。

真の効果: 1.3
真の傾向スコアを用いた場合の推定結果: 1.3345...
ロジスティック回帰を用いた場合の推定結果: 1.4002...
LightGBMを用いた場合の推定結果: 2.7601...

　LightGBMは介入$z$の予測精度こそ高いものの、ずいぶん大きく効果を推定してしまっていることがわかります。

結論

　当然ですが、真の傾向スコアを用いた推定が最も精度が高くなりました。つまり大事なのは、**「介入$z$の予測精度」ではなく「傾向スコアの精度（検証不可能）」**ということなのだと思います。とはいえ基本的に真の傾向スコアは知ることができません。したがって実務で傾向スコアを算出するようなときには以下のようなことに気を付ければよいのではないかと思いました。

1. 介入するか否かの意思決定者から基準を詳しくヒアリングする。****（＝データ生成機構に少しでも近いモデルを作るため）

2. 介入$z$の予測精度の向上に多大な労力を費やすのではなく、標準化平均差などを用いて共変量のバランスをモニタリングする。（今回は省略していますが、参考文献にはその重要性が説かれています。）
   　
   　最後までお読みくださりありがとうございました。

参考文献

• 効果検証入門　正しい比較のための因果推論/計量経済学の基礎（安井翔太　著）ISBN-10 : 4297111179　ISBN-13 : 978-4297111175