これはこのノートブックの日本語訳です。
from __future__ import division
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("whitegrid")
sns.set_palette("colorblind")
%matplotlib inline
import datagenerators as dg
import warnings
warnings.filterwarnings('ignore')
この投稿では、優れたCausalInferenceパッケージを使用して、観測データしかない状況について因果関係を推論しようとするフレームワークポテンシャル・アウトカムの使用方法の概要を説明します。著者は、その機能に関する一連のブログ投稿を書いています。
ダウンロードできるほとんどのデータセットは静的であるため、この投稿では、自分の関数を使用してデータを生成します。 これには2つの利点があります。特定のプロパティを持つデータセットを生成できること、および生成すること、そしてデータ生成システムに直接「介入」して、推論が正しいかどうかをチェックできることです。 これらのデータジェネレーターはすべて、いくつかの分布からi.i.d.サンプルを生成し、結果をpandasデータフレームとして返します。 これらのデータセットを生成する関数は、github ここの添付ファイル datagenerators.pyにあります。
まずは、やる気の出る例を見てみましょう。
導入
ある日、チームのリーダーが、チームの一部のメンバーがかっこいい帽子をかぶっていて、チームのこれらのメンバーの生産性が低下する傾向にあることに気付きました。 チームリーダーはデータに基づいて、チームメンバーがクールハットを着用しているかどうか（クールハットの場合は$X=1$、クールハットでない場合は$X=0$）、および生産的かどうか（$Y=1$は生産的、$Y=0$は非生産的）を記録し始めます。
1週間観察した結果、チームは次のようなデータセットを得ました。
observed_data_0 = dg.generate_dataset_0()
observed_data_0.head()
| x | y | |
|---|---|---|
| 0 | 0 | 0 |
| 1 | 1 | 1 |
| 2 | 0 | 0 |
| 3 | 0 | 1 |
| 4 | 0 | 0 |
チームリーダーが最初に尋ねる質問は「クールな帽子をかぶった人は、そうでない人よりも生産的である可能性が高いですか？」です。これは以下の数量を推定することを意味します。
$P(Y=1|X=1) - (Y=1|X=0)$
そして直接データからそれを推定することができます:
def estimate_uplift(ds):
"""
Estiamte the difference in means between two groups.
Parameters
----------
ds: pandas.DataFrame
a dataframe of samples.
Returns
-------
estimated_uplift: dict[Str: float] containing two items:
"estimated_effect" - the difference in mean values of $y$ for treated and untreated samples.
"standard_error" - 90% confidence intervals arround "estimated_effect"
"""
base = ds[ds.x == 0]
variant = ds[ds.x == 1]
delta = variant.y.mean() - base.y.mean()
delta_err = 1.96 * np.sqrt(
variant.y.var() / variant.shape[0] +
base.y.var() / base.shape[0])
return {"estimated_effect": delta, "standard_error": delta_err}
estimate_uplift(observed_data_0)
{'estimated_effect': -0.15738429037833684,
'standard_error': 0.08639322378194732}
クールな帽子をかぶった人は生産性が低いようです。
確かめるために、統計的検定を実行することもできます。
from scipy.stats import chi2_contingency
contingency_table = (
observed_data_0
.assign(placeholder=1)
.pivot_table(index="x", columns="y", values="placeholder", aggfunc="sum")
.values
)
_, p, _, _ = chi2_contingency(contingency_table, lambda_="log-likelihood")
# p-value
p
0.0005626147113456373
これは1つの小さなp値です。統計家は誇りに思うでしょう
この情報を使用して、クールな帽子をかぶっている人がいる場合、その人の確率についてどのように考えるかについて説明できます。以前の観察と「同じ分布から引き出されている」と私たちが信じている限り、同じ相関関係が存在すると予想されます。
問題は、この情報を、チームリーダーが人々にクールな帽子を被ることを強制すべきかどうかの議論として使おうとした場合です。もしチームリーダーがこのようなことをすると、私たちがサンプリングしているシステムを根本的に変えてしまい、以前に観察した相関関係を変更したり、逆にしてしまう可能性があります。
システムの一部の変更の影響を実際に測定する最もクリーンな方法は、ランダム化比較試験 を実行することです。具体的には、誰がクールな帽子を手に入れ、誰がそうしないかをランダム化し、得られる値 $y$ の違いを見ます。これにより、気になる指標に影響を及ぼしている可能性のある交絡変数 の影響が除去されます。
既知のプロセス（この場合は私が作成した関数）からデータセットを生成したので、それに直接介入してA/Bテストの効果を測定できます。
def run_ab_test(datagenerator, n_samples=10000, filter_=None):
"""
Generates n_samples from datagenerator with the value of X randomized
so that 50% of the samples recieve treatment X=1 and 50% receive X=0,
and feeds the results into `estimate_uplift` to get an unbiased
estimate of the average treatment effect.
Returns
-------
effect: dict
"""
n_samples_a = int(n_samples / 2)
n_samples_b = n_samples - n_samples_a
set_X = np.concatenate([np.ones(n_samples_a), np.zeros(n_samples_b)]).astype(np.int64)
ds = datagenerator(n_samples=n_samples, set_X=set_X)
if filter_ != None:
ds = ds[filter_(ds)].copy()
return estimate_uplift(ds)
run_ab_test(dg.generate_dataset_0)
{'estimated_effect': 0.18679999999999997,
'standard_error': 0.019256613018207393}
突然、かっこいい帽子をかぶる効果の方向性が逆転したようです。
何が起こっているのでしょうか？
注意：上記の例とこれに続く全ての例で、標本はi.i.d であり、SUTVA に従っていることを仮定しています。基本的にこれは、ある人が本当にかっこいい帽子を選ぶか、または被ることを強制されたとき、彼らは本当にかっこいい帽子を被っている別の人の選択や効果に影響を与えないことを意味します。構造上、私が使用している合成データジェネレーターはすべてこの性質を持っています。現実では、それはあなたが真実であると仮定しなければならないもう一つのことです。
因果関係の定義
先の例は、昔の統計家の格言を示しています:
"因果関係"は、漠然とした哲学的に響く単語です。 現在の文脈では、「 $X$ の変更が $Y$ に及ぼす影響は何ですか」という意味で使用しています。
正確には、$X$ と $Y$ はランダム変数 であり、知りたい「効果」は、$X$ に特定の値を代入すると $Y$ の分布がどのように変化するかです。 変数に特定の値を強制的に与えるこの行為は、「介入」と呼ばれます。
先ほどの例では、システムに何も介入しなかった場合、 $X$ を観測したことを条件に、 $Y$ の観測分布が得られます。
$P(Y|X)$
私たちが人々にかっこいい帽子を強制するとき、私たちは介入をしています。そして、$Y$の分布は、介入 分布によって与えられます。
$P(Y|\hbox{do}(X))$
一般的には、この二つは同じではありません。
これらのノートが試みて答えようとする質問は、観察データしか入手できない場合に、どのようにして介入分布について推論できるかということです。介入の効果を直接測定するためにA/Bテストを実行することが非現実的、実現不可能、または非倫理的な状況がたくさんあるので、これは有用な質問です。このような状況でも、介入の効果が何であるかについて何かを言えるようにしたいと思います。
ポテンシャル・アウトカム
この問題にアプローチする1つの方法は、2つの新しいランダム変数をシステムに導入することです： $Y_{0}$ と $Y_{1}$ で、ポテンシャル・アウトカムとして知られています。これらの変数は存在し、他のランダム変数と同じように扱うことができると仮定します。
$Y$ は、次のように定義されています。
- $Y = Y_{1}$ when $X=1$
- $Y = Y_{0}$ when $X=0$
これは、問題を介入の下で分布がどのように変化するかについてのものから、欠損値を持ついくつかの基礎となる分布からi.i.d.で描かれたデータについてのものへとシフトさせます。なぜ値が欠損しているのかについてのある仮定の下では、欠損値を推定する方法についてはよく発展した理論があります。
目標
しばしば、我々は完全な介入分布、$P(Y|\hbox{do}(X))$を気にしません、そして、2つのグループ間の平均の差の推定値があれば十分です。これは、平均処置効果(ATE)として知られている数値です。
$\Delta = E[Y_{1} - Y_{0}]$
A/Bテストを実行し、各グループの平均値を比較するとき、これは私たちが測定している数値に直結しています。
この数値を観測分布から推定してみると、以下のものが得られます。
$\Delta_{bad} = E[Y|X=1] - E[Y|X=0] \
= E[Y_{1}|X=1] - E[Y_{0}|X=0] \
\neq \Delta$
これは、一般的には真のATEとは一致しません。なぜならば
$E[Y_{i}|X=i] \neq E[Y_{i}]$
関連する2つの量は
- $ATT = E[Y_{1} - Y_{0}|X=1]$, the "Average Treatment effect of the Treated"
- $ATC = E[Y_{1} - Y_{0}|X=0]$, the "Average Treatment effect of the Control"
だからです。
ATCを解釈する一つの方法は、自然には処置されない標本のみを処理した場合の効果を測定するものであり、ATTの場合はその逆です。使用ケースにもよりますが、それらはあなたが気にしていることのより自然な尺度になるかもしれません。以下のテクニックを使えば、それらをすべて推定することができます。
$\def\ci{\perp!!!\perp}$
仮定をする
A/Bテストを行うときには、 $X$ の代入をランダムにします。これは、 $Y_{1}$ と $Y_{0}$ のどちらの変数を明らかにするかを選択できるようにする効果があります。これにより、結果は $X$ の値とは無関係になります。
これを次のように書きます。
$Y_{1}, Y_{0} \mathrel{\unicode{x2AEB}} X$
つまり、 $X, Y_{0}, Y_{1}$ の分布は、次のように因数分解されます。
$P(X, Y_{0}, Y_{1}) = P(X)P(Y_{0}, Y_{1})$
もし、この独立性が担保されるのであれば、次のようになります。
$E[Y_{1}|X=1] = E[Y_{1}]$
観察データを用いてATEを推定したい場合、標本について持っている他の情報を用いる必要があります。特に各対象の処置の選択を完全に説明するのに十分な追加情報があると仮定する必要があります。
その追加情報をランダム変数 $Z$ と呼ぶと、この仮定は次のように書けます。
$Y_{1}, Y_{0} \mathrel{\unicode{x2AEB}} X , | , Z$
あるいは
$P(X, Y_{0}, Y_{1}| Z) = P(X|Z)P(Y_{0}, Y_{1}|Z)$
これは、標本が受ける観察された治療 $X$ は、$Z$ によって完全に説明されることを意味します。これは、"無視可能性"仮定と呼ばれることもあります。
この例では、クールハットについての動機付けの例では、人の生産性とクールハットをかぶっているかどうかの両方に影響を与える他の要因（「スキル」と呼ぶことにしましょう）があることを意味します。上の例では、熟練した人は生産性が高く、クールハットをかぶっている可能性は低くなります。これらの事実は、A/Bテストを実行したときにクールハットの効果が逆転するように見えた理由を一緒に説明することが できるかも しれません。
もしその人が熟練しているかどうかでデータを分割すると、各サブグループでクールハットの着用と生産性の間に正の関係があることがわかります。
observed_data_0_with_confounders = dg.generate_dataset_0(show_z=True)
print(estimate_uplift(observed_data_0_with_confounders.loc[lambda df: df.z == 0]))
print(estimate_uplift(observed_data_0_with_confounders.loc[lambda df: df.z == 1]))
{'estimated_effect': 0.1647058823529412, 'standard_error': 0.15924048578984673}
{'estimated_effect': 0.14041297935103236, 'standard_error': 0.17328031714939449}
残念ながら、同じ標本で $Y_{0}$ と $Y_{1}$ を観測できないので、以下の仮定を検証することはできません。
$Y_{1}, Y_{0} \mathrel{\unicode{x2AEB}} X , | , Z$
それは、私たちが調査しているシステムの知識を使って評価しなければならないものです。
あなたが行う予測の質は、この仮定がどれだけうまく保持されているかにかかっています。シンプソンのパラドックスは、もし$Z$がすべての交絡変数を含まない場合、我々が行う推論が間違っている可能性があるという事実の極端な例です。Facebook社には、条件付き独立性が保持されていない場合に、どのように効果が過大評価されるかを示す、直接A/Bテストを用いた異なる因果推論アプローチを比較した良い論文を公開しています。
一旦この仮定をしてしまえば、これにアプローチするためのいくつかのテクニックがあります。この記事の残りの部分では、より簡単なアプローチのいくつかを概説しますが、これは現在進行中の研究分野であることを覚えておいてください。
反事実のモデル化
以上のことから、 $Y_{0}$ と $Y_{1}$ がわかれば、ATE を推定できることは明らかです。では、これらを直接モデル化してみてはどうでしょうか？
具体的には、以下のような推定量を作ることができます。
- $\hat{Y}_{0}(Z) = E[Y|Z, X=0]$
- $\hat{Y}_{1}(Z) = E[Y|Z, X=1]$
この２つの量をモデル化することができれば、ATEを次のように推定することができます。
このアプローチの成功は、潜在的な結果をどれだけうまくモデル化できるかにかかっています。実際にそれを見るために、以下のデータ生成プロセスを使用してみましょう。
observed_data_1 = dg.generate_dataset_1()
observed_data_1.plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False);
反事実のモデル化に入る前に、データを見てみましょう。$Y$がどのように分布しているかを見てみると、2つのグループの間には小さな違いがあるように見えます。
sns.kdeplot(observed_data_1.loc[lambda df: df.x == 0].y, label="untreated")
sns.kdeplot(observed_data_1.loc[lambda df: df.x == 1].y, label="treated")
<matplotlib.axes._subplots.AxesSubplot at 0x1a1a6c3590>
このことは、2つのグループ間の平均の差を見ることで確認することができます。
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_1)))
Observed ATE: 0.113 (0.103)
しかし、共分散である$Z$の分布を見てみると、群間で差があることがわかります。
sns.kdeplot(observed_data_1.loc[lambda df: df.x == 0].z, label="untreated")
sns.kdeplot(observed_data_1.loc[lambda df: df.x == 1].z, label="treated")
<matplotlib.axes._subplots.AxesSubplot at 0x1a1a7bced0>
もし $Z$ が $Y$ に何らかの影響を与えると考えるならば、このことは気になるところです。$X$が$Y$に与える影響と$Z$が$Y$に与える影響を区別する方法が必要です。
実際のATEをシミュレーションしたA/Bテストで確認し、それが観測値との差であることを確認することができます。
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(dg.generate_dataset_1)))
Real ATE: -0.479 (0.026)
しかし、このA/Bテストを実行できない場合はどうなるのでしょうか？システムをモデル化する必要があります。
最も単純なタイプのモデルは線形モデルです。具体的には、次のように仮定します。
$Y_{0} = \alpha + \beta Z + \epsilon$
$Y_{1} = Y_{0} + \gamma$
これが正確であれば、データをモデルにフィッティングさせて
$Y = \alpha + \beta Z + \gamma X$
を線形回帰を用いると、ATEの推定値が得られます。
causalinference パッケージはこれを行うためのシンプルなインターフェースを提供してくれます。
from causalinference import CausalModel
cm = CausalModel(
Y=observed_data_1.y.values,
D=observed_data_1.x.values,
X=observed_data_1.z.values)
cm.est_via_ols(adj=1)
print(cm.estimates)
Treatment Effect Estimates: OLS
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE -0.488 0.032 -15.122 0.000 -0.551 -0.425
causalinference はATEの推定値と、その推定値の統計的性質を返します。推定値に対して報告された信頼区間は、モデルが正確に反事実を記述していると仮定した場合の信頼区間であり、モデルが反事実をどの程度よく記述しているかについての信頼区間ではないことを認識することが重要です。
このケースでは、パッケージは正しいATEを識別することに成功しています - これは良いことですが、データ生成プロセスは仮定を満たすように特別に設計されています。失敗する可能性のあるいくつかのケースを見てみましょう。
最初のケースは、効果が単に相加的ではない場合です。
observed_data_2 = dg.generate_dataset_2()
observed_data_2.plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False)
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_2)))
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(dg.generate_dataset_2)))
Observed ATE: 0.689 (0.104)
Real ATE: 0.563 (0.031)
cm = CausalModel(
Y=observed_data_2.y.values,
D=observed_data_2.x.values,
X=observed_data_2.z.values)
cm.est_via_ols(adj=1)
print(cm.estimates)
Treatment Effect Estimates: OLS
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE 0.246 0.088 2.805 0.005 0.074 0.419
通常、これはより強力な推定量を使用することで克服することができます。シンプルでノンパラメトリックなアプローチは、マッチングの手法です。このアイデアは、処置を受けた各標本について、処置を受けなかった類似する標本を見つけ、これらの値を直接比較することです。「類似する」とは何を意味するのかは、あなたの特定の使用例に依存します。
causalinference パッケージは、各ユニットについて、他の処置群から最も類似したユニットを置換で選択し、ATEを計算するためにこれら2つのユニット間の差を使用することで、マッチングを実装しています。デフォルトでは、マッチングの選択は、共変量空間$Z$において最近傍法を用いて選択され、距離は各次元の逆分散で重み付けされます。
比較される単位の数や、一致における各次元の重み付けを変更するオプションがあります。詳細は、ドキュメントを参照してください。
マッチング推定値は、以下のコードで計算できます。
cm = CausalModel(
Y=observed_data_2.y.values,
D=observed_data_2.x.values,
X=observed_data_2.z.values)
cm.est_via_matching()
print(cm.estimates)
Treatment Effect Estimates: Matching
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE 0.480 0.148 3.239 0.001 0.190 0.770
ATC 0.973 0.142 6.854 0.000 0.695 1.251
ATT 0.036 0.211 0.169 0.866 -0.379 0.450
推定値の信頼区間には、真のATEが含まれています。
共変量の不均衡
対応するのがより困難な問題は，使用している共変量が不均衡な場合：共変量空間に，処理された標本または処理されていない標本のみを含む領域がある場合です。ここでは、治療の効果を外挿しなければなりません - これは、我々が使用する仮定モデルに大きく依存します。
以下の例は、これを実証しています。
observed_data_3 = dg.generate_dataset_3()
observed_data_3.plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False)
# actual response curves
z = np.linspace(0,1,100)
y0 = np.where(z >= 0.4, -4*(z - 0.4), 0)
y1 = np.where(z < 0.6, -4*(z - 0.6), 0) + 1
plt.plot(z,y0, "b")
plt.plot(z,y1, "r")
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_3)))
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(dg.generate_dataset_3)))
Observed ATE: 1.362 (0.080)
Real ATE: 2.449 (0.033)
# OLS estimator
cm = CausalModel(
Y=observed_data_3.y.values,
D=observed_data_3.x.values,
X=observed_data_3.z.values)
cm.est_via_ols()
print(cm.estimates)
Treatment Effect Estimates: OLS
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE 1.969 0.055 35.707 0.000 1.861 2.078
ATC 1.971 0.056 35.117 0.000 1.861 2.081
ATT 1.968 0.073 26.923 0.000 1.825 2.112
# Matching estimator
cm = CausalModel(
Y=observed_data_3.y.values,
D=observed_data_3.x.values,
X=observed_data_3.z.values)
cm.est_via_matching()
print(cm.estimates)
Treatment Effect Estimates: Matching
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE 1.913 0.156 12.269 0.000 1.607 2.218
ATC 1.965 0.277 7.093 0.000 1.422 2.509
ATT 1.876 0.177 10.626 0.000 1.530 2.222
OLS推定量では真の効果を捉えることができず、マッチング推定量では少し改善されますが、重なりがない領域に完全に外挿するにはデータ内に十分な情報がないだけです。
この例は矛盾しているように見えるかもしれませんが、一度高次元の共変量を見始めると、この問題はもっと一般的になります。
causalinference は、summary_stats プロパティを用いて変数の重複を素早く評価するための便利なツールを提供します。
print(cm.summary_stats)
Summary Statistics
Controls (N_c=206) Treated (N_t=294)
Variable Mean S.d. Mean S.d. Raw-diff
--------------------------------------------------------------------------------
Y -0.151 0.438 1.210 0.467 1.362
Controls (N_c=206) Treated (N_t=294)
Variable Mean S.d. Mean S.d. Nor-diff
--------------------------------------------------------------------------------
X0 0.267 0.185 0.681 0.205 2.116
ここで、正規化された差は次のように定義されます。
これは厳密な統計的検定ではありませんが、各共変量の間にどれだけの重複があるかを示すいくつかの指標を提供します。1より大きい値は、あまり重複がないことを示唆しています。
傾向スコア
傾向スコアは、共変量が与えられている場合に、対象が処置を受けることになった可能性がどれだけ高いかの推定値です。
$\hat{p}(Z) = P(X|Z)$
これは好きなように推定することができますが、一旦推定してしまえば、いくつかのことができるようになります。
逆傾向スコアによる重み付け
因果推論の測定の問題は、量 $E[Y_{i}]$ を知りたいのですが、 $E[Y_{i}|X=i]$ の標本にしかアクセスできないことを覚えておいてください。
ポテンシャル・アウトカムの確率は、次のように拡張することができます。
$P(Y_{i}) = P(Y_{i}| X = i)P(X = i)$
このことは、真の
$E[Y_{i}] = E[\frac{Y_{i}}{P(X=i|Z)}P(X=i|Z)] = E[\frac{Y_{i}}{P(X=i|Z)}|X=i, Z]$
を推定することができることを示唆しています。
そこで、各点に逆傾向スコアで重みを付ければ、潜在的な結果を回収することができます。その結果が 逆傾向スコア重み推定法 です。
$\Delta_{IPS} = \frac{1}{N}\left(\sum_{i \in 1} \frac{y_{i}}{\hat{p}(z_{i})} - \sum_{i \in 0} \frac{y_{i}}{1 - \hat{p}(z_{i})}\right)$
以前のデータセットの一つをどうするか見てみましょう。
observed_data_1 = dg.generate_dataset_1()
observed_data_1.plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False)
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_1)))
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(dg.generate_dataset_1)))
Observed ATE: 0.147 (0.102)
Real ATE: -0.490 (0.025)
CausalInferenceパッケージのメソッド est_propensity_s や est_propensity を用いると、共変量に対するロジスティック回帰を用いて傾向を推定することができます。
cm = CausalModel(
Y=observed_data_1.y.values,
D=observed_data_1.x.values,
X=observed_data_1.z.values)
cm.est_propensity_s()
propensity = cm.propensity["fitted"]
df = observed_data_1
df["ips"] = np.where(
df.x == 1,
1 / propensity,
1 / (1 - propensity))
df["ipsw"] = df.y * df.ips
ipse = (
df[df.x == 1]["ipsw"].sum()
- df[df.x == 0]["ipsw"].sum()
) / df.shape[0]
ipse
-0.4613297031604915
この例で使用しているデータ生成器では、単純なロジスティック回帰で関係をうまく記述することができます。この例で使用しているデータ生成器では、関係性は単純なロジスティック回帰で十分に記述できますが、もし我々がskleanのロジスティック回帰関数（デフォルトでは正則化が使用されています）を使用して傾向スコアを推定しようとした場合、間違った答えが得られるでしょう。
from sklearn.linear_model import LogisticRegression
lg = LogisticRegression()
X = df.z.values.reshape(-1,1)
y = df.x.values
lg.fit(X,y)
propensity = lg.predict_proba(X)[:,1]
df["ips"] = np.where(
df.x == 1,
1 / propensity,
1 / (1 - propensity))
df["ipsw"] = df.y * df.ips
ipse = (
df[df.x == 1]["ipsw"].sum()
- df[df.x == 0]["ipsw"].sum()
) / df.shape[0]
ipse
-0.3291860760196773
我々の素朴な推定値よりは良いですが、正しくありません。
頑健な重み付き推定器
このように、逆傾向スコアの重み付け推定量と効果量の線形推定量を組み合わせて、どちらかのモデルの欠陥を減らしてみることができます。これは、各ポイントを逆傾向スコアで重み付けした状態で、データに重み付き線形回帰を事前に実行することで行われます。その結果がdoubly robust weighted estimatorです。
考え方としては、観測データの中でどの標本が処置されているかにバイアスがあるので、標本の中でも処置されたが、処置される可能性が低かった標本の方がより重要で、より多くの重みを与えるべきだという点です。
これを応用するには、以下のような方法があります。
observed_data_1 = dg.generate_dataset_1()
observed_data_1.plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False)
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_1)))
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(dg.generate_dataset_1)))
Observed ATE: 0.145 (0.106)
Real ATE: -0.494 (0.026)
cm = CausalModel(
Y=observed_data_1.y.values,
D=observed_data_1.x.values,
X=observed_data_1.z.values)
cm.est_propensity_s()
cm.est_via_weighting()
print(cm.estimates)
Treatment Effect Estimates: Weighting
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE -0.521 0.034 -15.344 0.000 -0.588 -0.455
cm.estimates
{'weighting': {'ate': -0.5214252844257842, 'ate_se': 0.03398288940994425}}
Unconfoundedness(非交絡性)と傾向スコア
前のセクションでは、共変量が与えられているので、結果と処置は独立していると仮定しました。
$Y_{1}, Y_{0} \mathrel{\unicode{x2AEB}} X , | ,Z$
また、もう少し強い仮定もできます。それは、結果が処置とは無関係で、傾向の確率を条件にしているということです。
$Y_{1}, Y_{0} \mathrel{\unicode{x2AEB}} X , | ,\hat{p}(Z)$
この仮定では、交絡変数の次元を減らすことができます。これは、高次元の設定では機能しないかもしれないいくつかの手法を実行することを可能にします。
トリミング
共変量の不均衡が問題を引き起こすことがあることをこれまでに確認しました。簡単な解決策は、良いオーバーラップがある領域でのみ反事実の予測を行うか、または良いオーバーラップがない点を "トリム "することです。高次元データの場合、"良いオーバーラップ "を定義するのは難しいかもしれません - オーバーラップを定義するために傾向スコアを使用することは、これを解決する1つの方法です。
オーバーラップの少ないデータセットを見てみましょう。
observed_data_3 = dg.generate_dataset_3()
observed_data_3.plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False)
# actual response curves
z = np.linspace(0,1,100)
y0 = np.where(z >= 0.4, -4*(z - 0.4), 0)
y1 = np.where(z < 0.6, -4*(z - 0.6), 0) + 1
plt.plot(z,y0, "b")
plt.plot(z,y1, "r")
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_3)))
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(dg.generate_dataset_3)))
Observed ATE: 1.361 (0.080)
Real ATE: 2.437 (0.033)
CausalInferenceは、傾向スコアに基づいてデータをトリミングする方法を提供しています。
# OLS estimator
cm = CausalModel(
Y=observed_data_3.y.values,
D=observed_data_3.x.values,
X=observed_data_3.z.values)
cm.est_propensity_s()
cm.trim_s()
cm.est_via_matching()
print(cm.estimates)
Treatment Effect Estimates: Matching
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE 2.068 0.127 16.350 0.000 1.820 2.316
ATC 2.249 0.229 9.843 0.000 1.802 2.697
ATT 1.920 0.102 18.739 0.000 1.719 2.121
残りのデータを次のように見ることができます。
# mask out data ignored by the
propensity = cm.propensity["fitted"]
cutoff = cm.cutoff
mask = (propensity > cutoff) & (propensity < 1 - cutoff)
# plot the data
observed_data_3[mask].plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False)
# actual response curves
z = np.linspace(0,1,100)
y0 = np.where(z >= 0.4, -4*(z - 0.4), 0)
y1 = np.where(z < 0.6, -4*(z - 0.6), 0) + 1
plt.plot(z,y0, "b")
plt.plot(z,y1, "r")
# filter out data in regions we have trimmed when we calculate the true uplift
filter_ = lambda df: (df.z > 0.2) & (df.z < 0.7)
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_3[mask])))
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(
**run_ab_test(dg.generate_dataset_3, filter_= filter_)))
Observed ATE: 1.708 (0.125)
Real ATE: 1.984 (0.031)
これは、そのような場合にはかなりうまくいきます。
トリミングを適用するとき、共変量空間の一部のサンプルについてのみ因果推論を行うことが可能であることを明示的に述べています。これらの領域外の標本について、ATEについては何も言えません。
層化
傾向スコアのもう1つの用途は、層別化またはブロック化推定量です。 これは、データポイントを同様の傾向のグループにグループ化し、これらのグループ内のATEを推定することで構成されます。 繰り返しますが、 CausalInferenceはこれを実現するための優れたインターフェースを提供します。
stratify（ユーザー定義のstata境界の場合）または stratify_s（自動的に境界を選択する）メソッドを使用して、層を決定します。
observed_data_1 = dg.generate_dataset_1()
cm = CausalModel(
Y=observed_data_1.y.values,
D=observed_data_1.x.values,
X=observed_data_1.z.values)
cm.est_propensity_s()
cm.stratify_s()
print(cm.strata)
Stratification Summary
Propensity Score Sample Size Ave. Propensity Outcome
Stratum Min. Max. Controls Treated Controls Treated Raw-diff
--------------------------------------------------------------------------------
1 0.093 0.260 108 18 0.163 0.182 -0.544
2 0.266 0.327 23 9 0.296 0.292 -0.639
3 0.327 0.411 14 17 0.368 0.376 -0.441
4 0.412 0.555 31 31 0.473 0.490 -0.428
5 0.557 0.769 45 80 0.661 0.667 -0.499
6 0.770 0.910 18 106 0.838 0.859 -0.536
そして est_via_blocking メソッドを用いて、これらの層の推定値を1つの全体的なATEに結合します。
cm.est_via_blocking()
print(cm.estimates)
Treatment Effect Estimates: Blocking
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE -0.559 0.033 -16.722 0.000 -0.625 -0.494
ATC -0.571 0.039 -14.797 0.000 -0.646 -0.495
ATT -0.548 0.038 -14.304 0.000 -0.624 -0.473
これもよく機能します。
傾向スコアによってデータをグループに層別化することは、何が「類似した」単位を構成するかについて予備知識がない場合に有用ですが、それだけが唯一の方法ではありません。標本の異なるグループが類似した方法で介入の影響を受ける可能性が高いことを事前に知っている場合は、ATEを推定する前にサンプルをこれらのグループに分割し、結果をプールしてグローバルなATEを取得することは理にかなっています。
どの手法を使うか？
これで、観測データからの因果推論のための一般的な手法のほとんどをカバーしました。残りの問題は「どのように、どの方法を使用するかを決定するか？」です。これは簡単な質問ではありません。この論文のような自動化された手法もありますが、私はそれらを試してみたことがありません。
最終的には、あなたのテクニックを選択するためには、あなたがどのように反事実を構築するかについていくつかの仮定をする必要があります。共変量空間での重なりが良好であるとデータを信頼しているならば、マッチングは良いアプローチです。これがそうでない場合は、未踏の領域にうまく外挿するために信頼できるモデルを使用するか、傾向スコアのようなものが無視可能性を仮定するのに十分な情報を捉えるという仮定を立てる必要があります。
これらの方法がすべて失敗する可能性があることを強調するために、もう一つ例を挙げます。前の例とは異なり、共変量は1つ以上ある。これまでのすべてのデータ・ジェネレーターと同様に、このデータ・ジェネレーターもまた、設計上以下の仮定に従います。
$Y_{1}, Y_{0} \mathrel{\unicode{x2AEB}} X , | ,Z$
これまで議論してきた方法を盲目的に試して、どうなるか見てみましょう。
data_gen = dg.generate_exercise_dataset_2
ds = data_gen()
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(ds)))
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(data_gen)))
zs = [c for c in ds.columns if c.startswith("z")]
cm = CausalModel(
Y=ds.y.values,
D=ds.x.values,
X=ds[zs].values)
cm.est_via_ols()
cm.est_via_matching()
cm.est_propensity_s()
cm.est_via_weighting()
cm.stratify_s()
cm.est_via_blocking()
print(cm.estimates)
Observed ATE: -0.423 (0.711)
Real ATE: 4.570 (0.184)
Treatment Effect Estimates: OLS
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE -0.165 0.363 -0.455 0.649 -0.878 0.547
ATC 0.438 0.402 1.088 0.276 -0.350 1.226
ATT -0.468 0.384 -1.218 0.223 -1.220 0.285
Treatment Effect Estimates: Matching
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE 1.241 0.581 2.137 0.033 0.103 2.379
ATC 1.314 0.695 1.890 0.059 -0.049 2.677
ATT 1.204 0.667 1.804 0.071 -0.104 2.512
Treatment Effect Estimates: Weighting
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE 3.190 0.478 6.667 0.000 2.252 4.128
Treatment Effect Estimates: Blocking
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE -0.003 0.378 -0.009 0.993 -0.745 0.738
ATC -0.003 0.378 -0.009 0.993 -0.745 0.738
ATT -0.003 0.378 -0.009 0.993 -0.745 0.738
y = []
yerr = []
x_label = []
for method, result in dict(cm.estimates).items():
y.append(result["ate"])
yerr.append(result["ate_se"])
x_label.append(method)
x = np.arange(len(y))
plt.errorbar(x=x, y=y, yerr=yerr, linestyle="none", capsize=5, marker="o")
plt.xticks(x, x_label)
plt.title("Estimated Effect Size", fontsize=18)
plt.hlines(4.6, -0.5, 3.5, linestyles="dashed")
plt.xlim(-0.5,3.5);
横の破線は、このデータセットの真のATEを示しています。
それぞれの手法の結果が異なるだけでなく、すべての手法が真の値を見逃しています。
これは、この種の手法の限界についての警告であるべきです。データセットのどのような性質が、これらの手法が真の値を外れる原因となっているのか、読者が試してみるのは、非常に興味深い練習になるかもしれません。
因果推論の構造
願わくば、この時点までに、あなたは無視可能性仮定の重要性に気付いているでしょう。
$Y_{1}, Y_{0} \mathrel{\unicode{x2AEB}} X , | ,Z$
これまで話してこなかったのは、これが真実であるように、どのようにして$Z$を選択するかということです。最終的には、研究されているシステムに関する領域の知識から来る必要があります。因果関係モデルと呼ばれる強力なツールのセットがあります。これを使えば、研究しているシステムに関する知識をシステムのグラフィカルモデルとしてエンコードしたり、上記のような条件付き独立性の仮定について推論したりすることができます。
この投稿が提起するもう一つの疑問は、因果推論を行う唯一の方法は交絡変数を調整することであるかどうかということです。そうではありません - 後の投稿では 他にも使えるテクニックを見てみます。
コード
この投稿のノートはgithub こちらにあります。