计算因果效应的基本示例#
这是对DoWhy因果推断库的快速介绍。我们将加载一个样本数据集，并估计一个（预先指定的）处理变量对一个（预先指定的）结果变量的因果效应。
首先，让我们加载所有需要的包。
[1]:
import numpy as np
from dowhy import CausalModel
import dowhy.datasets
现在，让我们加载一个数据集。为了简单起见，我们模拟了一个数据集，其中常见原因与治疗之间以及常见原因与结果之间存在线性关系。
Beta 是真实的因果效应。
[2]:
data = dowhy.datasets.linear_dataset(beta=10,
num_common_causes=5,
num_instruments = 2,
num_effect_modifiers=1,
num_samples=5000,
treatment_is_binary=True,
stddev_treatment_noise=10,
num_discrete_common_causes=1)
df = data["df"]
[3]:
df.head()
[3]:
| X0 | Z0 | Z1 | W0 | W1 | W2 | W3 | W4 | v0 | y | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0.229215 | 1.0 | 0.803310 | 0.969484 | 0.708693 | -0.472480 | 1.134963 | 0 | True | 11.467510 |
| 1 | -1.424967 | 0.0 | 0.163824 | -0.161973 | -0.495899 | -0.971206 | -0.177170 | 0 | False | -2.703782 |
| 2 | 0.462122 | 1.0 | 0.287209 | 2.114670 | -1.484899 | -0.760696 | 1.512796 | 1 | True | 13.401639 |
| 3 | -2.460213 | 0.0 | 0.188397 | 0.524076 | 0.090692 | -1.059184 | -1.477797 | 1 | False | 0.378057 |
| 4 | -2.152985 | 1.0 | 0.276386 | 1.981324 | 0.648104 | -1.882990 | -0.061003 | 0 | True | 9.043652 |
请注意，我们正在使用pandas dataframe来加载数据。目前，DoWhy仅支持pandas dataframe作为输入。
接口 1 (推荐): 输入因果图#
我们现在以GML图形格式（推荐）输入因果图。你也可以使用DOT格式。
要为您的数据集创建因果图，您可以使用像DAGitty这样的工具，它提供了一个图形用户界面来构建图。您可以导出它生成的图字符串。图字符串非常接近DOT格式：只需将dag重命名为digraph，删除换行符并在每行后添加分号，即可将其转换为DOT格式并输入到DoWhy中。
[4]:
# With graph
model=CausalModel(
data = df,
treatment=data["treatment_name"],
outcome=data["outcome_name"],
graph=data["gml_graph"]
)
[5]:
model.view_model()
[6]:
from IPython.display import Image, display
display(Image(filename="causal_model.png"))
上述因果图展示了因果模型中编码的假设。我们现在可以使用这个图首先识别因果效应（从因果估计量到概率表达式），然后估计因果效应。
DoWhy 哲学：保持识别和估计分离#
识别可以在不访问数据的情况下实现，仅访问图。这导致需要计算一个表达式。然后可以在估计步骤中使用可用数据评估此表达式。重要的是要理解这些是正交的步骤。
识别#
[7]:
identified_estimand = model.identify_effect(proceed_when_unidentifiable=True)
print(identified_estimand)
Estimand type: EstimandType.NONPARAMETRIC_ATE
### Estimand : 1
Estimand name: backdoor
Estimand expression:
d
─────(E[y|W3,W4,W1,W0,W2])
d[v₀]
Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W3,W4,W1,W0,W2,U) = P(y|v0,W3,W4,W1,W0,W2)
### Estimand : 2
Estimand name: iv
Estimand expression:
⎡ -1⎤
⎢ d ⎛ d ⎞ ⎥
E⎢─────────(y)⋅⎜─────────([v₀])⎟ ⎥
⎣d[Z₁ Z₀] ⎝d[Z₁ Z₀] ⎠ ⎦
Estimand assumption 1, As-if-random: If U→→y then ¬(U →→{Z1,Z0})
Estimand assumption 2, Exclusion: If we remove {Z1,Z0}→{v0}, then ¬({Z1,Z0}→y)
### Estimand : 3
Estimand name: frontdoor
No such variable(s) found!
注意参数标志 proceed_when_unidentifiable。它需要设置为 True 以传达我们忽略任何未观察到的混杂因素的假设。默认行为是提示用户再次确认未观察到的混杂因素可以被忽略。
估计#
[8]:
causal_estimate = model.estimate_effect(identified_estimand,
method_name="backdoor.propensity_score_stratification")
print(causal_estimate)
*** Causal Estimate ***
## Identified estimand
Estimand type: EstimandType.NONPARAMETRIC_ATE
### Estimand : 1
Estimand name: backdoor
Estimand expression:
d
─────(E[y|W3,W4,W1,W0,W2])
d[v₀]
Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W3,W4,W1,W0,W2,U) = P(y|v0,W3,W4,W1,W0,W2)
## Realized estimand
b: y~v0+W3+W4+W1+W0+W2
Target units: ate
## Estimate
Mean value: 9.666348684992126
你可以向estimate_effect方法输入额外的参数。例如，要估计对任何单位子集的影响，你可以指定“target_units”参数，它可以是一个字符串（“ate”，“att”或“atc”），过滤数据框行的lambda函数，或计算效果的新数据框。你也可以指定“效果修饰符”来估计这些变量之间的异质效果。参见help(CausalModel.estimate_effect)。
[9]:
# Causal effect on the control group (ATC)
causal_estimate_att = model.estimate_effect(identified_estimand,
method_name="backdoor.propensity_score_stratification",
target_units = "atc")
print(causal_estimate_att)
print("Causal Estimate is " + str(causal_estimate_att.value))
*** Causal Estimate ***
## Identified estimand
Estimand type: EstimandType.NONPARAMETRIC_ATE
### Estimand : 1
Estimand name: backdoor
Estimand expression:
d
─────(E[y|W3,W4,W1,W0,W2])
d[v₀]
Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W3,W4,W1,W0,W2,U) = P(y|v0,W3,W4,W1,W0,W2)
## Realized estimand
b: y~v0+W3+W4+W1+W0+W2
Target units: atc
## Estimate
Mean value: 9.670293023783838
Causal Estimate is 9.670293023783838
接口2：指定常见原因和工具#
[10]:
# Without graph
model= CausalModel(
data=df,
treatment=data["treatment_name"],
outcome=data["outcome_name"],
common_causes=data["common_causes_names"],
effect_modifiers=data["effect_modifier_names"])
[11]:
model.view_model()
[12]:
from IPython.display import Image, display
display(Image(filename="causal_model.png"))
我们得到了相同的因果图。现在像以前一样进行识别和估计。
[13]:
identified_estimand = model.identify_effect(proceed_when_unidentifiable=True)
[14]:
estimate = model.estimate_effect(identified_estimand,
method_name="backdoor.propensity_score_stratification")
print(estimate)
print("Causal Estimate is " + str(estimate.value))
*** Causal Estimate ***
## Identified estimand
Estimand type: EstimandType.NONPARAMETRIC_ATE
### Estimand : 1
Estimand name: backdoor
Estimand expression:
d
─────(E[y|W3,W4,W1,W0,W2])
d[v₀]
Estimand assumption 1, Unconfoundedness: If U→{v0} and U→y then P(y|v0,W3,W4,W1,W0,W2,U) = P(y|v0,W3,W4,W1,W0,W2)
## Realized estimand
b: y~v0+W3+W4+W1+W0+W2
Target units: ate
## Estimate
Mean value: 9.666348684992126
Causal Estimate is 9.666348684992126
反驳估计#
现在让我们看看如何反驳所获得的估计。反驳方法提供了每个正确估计器都应通过的测试。因此，如果估计器未能通过反驳测试（p值<0.05），则意味着估计器存在一些问题。
请注意，我们无法验证估计是否正确，但如果它违反了某些预期行为，我们可以拒绝它（这类似于可以被证伪但不能被证明为真的科学理论）。以下的反驳测试基于以下两种方法之一：1）不变变换：数据的变化不应改变估计结果。任何在原始数据和修改数据之间结果显著变化的估计器都将无法通过测试；
随机共同原因
数据子集
无效化转换：数据变化后，因果真实估计为零。任何在新数据上结果显著偏离零的估计器都未能通过测试。
安慰剂治疗
添加一个随机共同原因变量#
[15]:
res_random=model.refute_estimate(identified_estimand, estimate, method_name="random_common_cause", show_progress_bar=True)
print(res_random)
Refute: Add a random common cause
Estimated effect:9.666348684992126
New effect:9.66634868499213
p value:1.0
用随机（安慰剂）变量替换治疗#
[16]:
res_placebo=model.refute_estimate(identified_estimand, estimate,
method_name="placebo_treatment_refuter", show_progress_bar=True, placebo_type="permute")
print(res_placebo)
Refute: Use a Placebo Treatment
Estimated effect:9.666348684992126
New effect:0.029698445089780416
p value:0.94
移除数据的随机子集#
[17]:
res_subset=model.refute_estimate(identified_estimand, estimate,
method_name="data_subset_refuter", show_progress_bar=True, subset_fraction=0.9)
print(res_subset)
Refute: Use a subset of data
Estimated effect:9.666348684992126
New effect:9.702493576161388
p value:0.5800000000000001
如你所见，倾向评分分层估计器对反驳具有合理的鲁棒性。
可重复性：为了确保可重复性，您可以向任何反驳方法添加一个参数“random_seed”，如下所示。
并行化：您还可以使用内置的并行化功能来加速反驳过程。只需将n_jobs设置为大于1的值，以将工作负载分配到多个CPU上，或者将n_jobs=-1设置为使用所有CPU。目前，这仅适用于random_common_cause、placebo_treatment_refuter和data_subset_refuter。
[18]:
res_subset=model.refute_estimate(identified_estimand, estimate,
method_name="data_subset_refuter", show_progress_bar=True, subset_fraction=0.9, random_seed = 1, n_jobs=-1, verbose=10)
print(res_subset)
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 5 tasks | elapsed: 3.1s
[Parallel(n_jobs=-1)]: Done 10 tasks | elapsed: 3.6s
[Parallel(n_jobs=-1)]: Done 17 tasks | elapsed: 4.4s
[Parallel(n_jobs=-1)]: Done 24 tasks | elapsed: 4.9s
[Parallel(n_jobs=-1)]: Done 33 tasks | elapsed: 6.0s
[Parallel(n_jobs=-1)]: Done 42 tasks | elapsed: 6.9s
[Parallel(n_jobs=-1)]: Done 53 tasks | elapsed: 8.1s
[Parallel(n_jobs=-1)]: Done 64 tasks | elapsed: 9.1s
[Parallel(n_jobs=-1)]: Done 77 tasks | elapsed: 10.6s
[Parallel(n_jobs=-1)]: Done 90 tasks | elapsed: 12.0s
Refute: Use a subset of data
Estimated effect:9.666348684992126
New effect:9.701936861633671
p value:0.6000000000000001
[Parallel(n_jobs=-1)]: Done 100 out of 100 | elapsed: 12.9s finished
添加一个未观察到的共同原因变量#
这种反驳不返回p值。相反，它提供了一个敏感性测试，用于评估如果识别假设（在identify_effect中使用）无效时，估计值变化的速度。具体来说，它检查对违反后门假设的敏感性：即所有共同原因都被观察到。
为此，它创建了一个新的数据集，其中在治疗和结果之间增加了一个共同的因果因素。为了捕捉共同因果因素的影响，该方法将共同因果因素对治疗和结果的影响强度作为输入。基于这些关于共同因果因素影响的输入，它改变了治疗和结果的值，然后重新运行估计器。希望新的估计值在未观察到的共同因果因素影响较小的情况下不会发生剧烈变化，这表明对任何未观察到的混杂因素具有鲁棒性。
解释此过程的另一种等效方式是假设输入数据中已经存在未观察到的混杂因素。处理和结果值的变化消除了原始数据中存在的任何未观察到的共同原因的影响。然后在此修改后的数据上重新运行估计器，提供了正确的识别估计值，我们希望对于未观察到的共同原因影响的某个有界值，新估计值与原始估计值之间的差异不会太大。
领域知识的重要性：此测试需要领域知识来设置未观察到的混杂效应的合理输入值。我们首先展示了混杂因素对治疗和结果影响的单一值的结果。
[19]:
res_unobserved=model.refute_estimate(identified_estimand, estimate, method_name="add_unobserved_common_cause",
confounders_effect_on_treatment="binary_flip", confounders_effect_on_outcome="linear",
effect_strength_on_treatment=0.01, effect_strength_on_outcome=0.02)
print(res_unobserved)
Refute: Add an Unobserved Common Cause
Estimated effect:9.666348684992126
New effect:8.664434813571749
随着未观察到的混杂效应增加，检查趋势通常更有用。为此，我们可以提供一系列假设的混杂效应。输出是在不同未观察到的混杂情况下估计效应的(最小, 最大)范围。
[20]:
res_unobserved_range=model.refute_estimate(identified_estimand, estimate, method_name="add_unobserved_common_cause",
confounders_effect_on_treatment="binary_flip", confounders_effect_on_outcome="linear",
effect_strength_on_treatment=np.array([0.001, 0.005, 0.01, 0.02]), effect_strength_on_outcome=0.01)
print(res_unobserved_range)
Refute: Add an Unobserved Common Cause
Estimated effect:9.666348684992126
New effect:(7.184834387565013, 9.395788216420343)
上述图表显示了随着假设的治疗混杂因素增加，估计值如何减少。根据领域知识，我们可能知道治疗的最大合理混杂效应。由于我们看到效应没有超过零，我们可以安全地得出结论，治疗v0的因果效应是正的。
我们还可以改变对治疗和结果的混淆效应。我们获得了一个热图。
[21]:
res_unobserved_range=model.refute_estimate(identified_estimand, estimate, method_name="add_unobserved_common_cause",
confounders_effect_on_treatment="binary_flip", confounders_effect_on_outcome="linear",
effect_strength_on_treatment=[0.001, 0.005, 0.01, 0.02],
effect_strength_on_outcome=[0.001, 0.005, 0.01,0.02])
print(res_unobserved_range)
Refute: Add an Unobserved Common Cause
Estimated effect:9.666348684992126
New effect:(3.2311645505216724, 9.63495028468222)
自动推断效应强度参数。 最后，DoWhy 支持自动选择效应强度参数。这是基于一个假设，即未观察到的混杂因素对治疗或结果的影响不能比任何观察到的混杂因素更强。也就是说，我们至少已经收集了最相关的混杂因素的数据。如果是这种情况，那么我们可以通过观察到的混杂因素的效应强度来限制 effect_strength_on_treatment 和 effect_strength_on_outcome 的范围。还有一个额外的可选参数，表示未观察到的混杂因素的效应强度是否应该与观察到的最高效应强度相同，或者是其一部分。您可以使用可选的 effect_fraction_on_treatment 和 effect_fraction_on_outcome 参数来设置它。默认情况下，这两个参数为 1。
[22]:
res_unobserved_auto = model.refute_estimate(identified_estimand, estimate, method_name="add_unobserved_common_cause",
confounders_effect_on_treatment="binary_flip", confounders_effect_on_outcome="linear")
print(res_unobserved_auto)
/github/home/.cache/pypoetry/virtualenvs/dowhy-oN2hW5jr-py3.8/lib/python3.8/site-packages/sklearn/utils/validation.py:1183: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().
y = column_or_1d(y, warn=True)
Refute: Add an Unobserved Common Cause
Estimated effect:9.666348684992126
New effect:(0.7690576672487185, 9.398648715588278)
结论: 假设未观察到的混杂因素对处理或结果的影响不超过任何已观察到的混杂因素，则可以得出因果效应为正的结论。