解释因果推断中的“do”操作符

引言

在互联网上或书籍中，很难找到完全解释的、端到端示例的因果推断，其中包含实际的、可工作的源代码，正如我在探索这项新兴技术的工作原理以及其重要性的过程中所发现的那样。

但是，如果你坚持下去，这绝对是值得的，因为它能够解决其他机器学习技术无法解决的问题类型。

传统的机器学习模型可以预测如果未来与过去类似，可能会发生什么，但它们无法告诉你应该采取哪些不同的行动来实现期望的结果。

例如，分类算法可以预测银行贷款客户是否可能违约，但它无法回答诸如“如果我们改变贷款的还款期限，是否会有更多客户避免违约？”这样的问题。

以下是因果推断可以回答的一些传统预测模型无法回答的问题类型的示例：

• 一个对系统的建议性变更是否改善了人们的结果？

• 导致系统结果变化的原因是什么？

• 对系统的哪些变更可能会改善人们的结果？

有许多在线文章详细介绍了因果推断涉及的数学细节，但很少有提供完整解释和所有源代码的实际示例。

如果你继续阅读本文直到结束，我承诺将为你提供完整的解释和所有源代码，使你能够使用其他机器学习技术无法实现的方式做一些真正了不起的事情。

数据

我们首先需要一些数据。我创建了一个纯粹的合成数据集，受到了著名的LaLonde数据的启发，该数据观察和记录了20世纪70年代就业技能培训计划对收入的影响。

由于LaLonde数据和研究提供了灵感，所以在文章末尾的参考部分有一个引用。

import pandas as pd

df_training=pd.read_excel("data/training.xlsx")
df_training["age_group"] = df_training["age_group"].astype("category")
df_training

值得花点时间了解合成数据集的关键方面

• received_training如果个体参加了旨在提供就业技能并增加收入潜力的虚构培训计划，则为1。在合成数据集中，有640人参加了培训计划，1360人没有参加。

• age是年龄，以年为单位。

• education_years持有学校教育年限。

• received_benefits如果个体曾经领取过失业救济，则为1。

• university_degree如果个体在大学学习并获得学位，则为1。

• single如果个体单身（即未婚或未注册的伴侣关系），则为1。

• top_earner如果个体属于收入最高的四分之一人群，则为1。

• age_group是年龄的分类版本。

• earnings是个体在虚构就业技能培训计划完成后3年内的收入金额，是“目标”或感兴趣的特征。

培训计划产生了什么影响？

现在让我们来看看数据，看看培训计划对参与者收入产生了什么影响…

received_training_filter = df_training["received_training"] == 1

impact_of_training = df_training[received_training_filter]["earnings"].mean() - df_training[~received_training_filter]["earnings"].mean()

top_earner_received_trainingment = df_training[received_training_filter]["top_earner"].value_counts(normalize=True)
top_earner_no_received_trainingment = df_training[~received_training_filter]["top_earner"].value_counts(normalize=True)

p_top_earner_received_trainingment = top_earner_received_trainingment[1]
p_top_earner_no_received_trainingment = top_earner_no_received_trainingment[1]

print(f"The average impact of participation in the employment skills training program on earnings is ${impact_of_training:+0,.2f}\n")
print(f"P(top earner=1 | received_trainingment=1) = {p_top_earner_received_trainingment}")
print(f"P(top earner=1 | received_trainingment=0) = {p_top_earner_no_received_trainingment}")

display(df_training.groupby("received_training")["earnings"].agg(["median","mean"]))

根据分析，参加培训计划的影响是负面的 -

• 参加培训计划的明显影响是年收入减少了1065.29美元。

• 参加培训的人中，成为高收入者的概率为0.19，而不参加培训的人中为0.28。

• 接受培训的人的收入中位数为3739美元，未接受培训的人为4893美元。

• 接受培训的人的平均收入为6067美元，未接受培训的人为7132美元。

你会给出什么建议？

基于分析，明确的建议是停止培训计划，因为使用四种不同的测量方法可以证明对收入的影响始终是负面的。

然而，从直觉上看，这个结论似乎不正确。即使培训非常糟糕，参加培训也不应该会使参与者的就业能力降低并损害他们未来的收入。

那么，分析出了什么问题呢？

此时，传统方法，包括概率和预测模型，在这方面已经无法再提供更多帮助。任何这些技术的应用都会得出应取消培训的结论。

为了突破这些限制，真正理解发生了什么，我们需要建立一个因果模型，并应用神奇的“do”操作符。

如果你想了解培训计划的真实影响，以及为什么使用传统概率和预测模型可能会导致错误甚至危险的结果，请继续阅读…

使用“do”操作符的因果推断解决方案

让我们通过更详细地查看数据集中的一些特征来开始更准确的评估之旅…

import matplotlib.pyplot as plt

def plot_comparison(feature : str, normalize : bool = True):

    fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(15, 6))

    df_training[~received_training_filter][feature].value_counts(normalize=normalize).sort_index().plot(ax=axes[0], kind="bar", title="No Training", xlabel=feature, ylim=(0,1) if normalize else False)
    df_training[received_training_filter][feature].value_counts(normalize=normalize).sort_index().plot(ax=axes[1], kind="bar", title="Received Training", xlabel=feature, ylim=(0,1) if normalize else False)

    plt.show()

plot_comparison("education_years", normalize=False)

显然，在参加培训和未参加培训的人之间，在教育方面有明显不同的模式。

在一个现实世界的项目中，我们会与领域专家合作，以了解这些模式，但即使没有领域专业知识，合理地得出结论是，教育可能对参与培训的人员和他们的收入能力产生因果影响。

plot_comparison("age_group", normalize=False)

年龄也有类似的情况。对于接受培训的人，年龄呈线性模式，年轻人更多，老年人较少。对于没有接受培训的人，30至40岁年龄组有一个峰值。一个假设可能是，许多30至40岁的人已经赚得不错，不需要任何与工作技能相关的培训。

同样，这表明年龄影响了个体是否有可能参加培训以及他们的收入潜力。

因果推断和“混杂变量”

这个简单的附加分析揭示了“混杂变量”的存在。这个术语在许多现有的示例中被随意使用，通常伴随着微积分和公式，但通常没有清楚的解释。

简单地说，诸如年龄和教育之类的特征的影响与感兴趣的主要影响（即培训对收入的影响）混在一起，当我们应用传统方法时，这种单独的影响无法分离出来。

构建“有向无环图”（DAG）

单独从数据中无法发现因果关系。数据需要通过“有向无环图”（DAG）进行补充，该图通过利用领域专业知识和其他技术来“发现”因果关系。

有关更详细的探讨，请参阅我的有关发现因果关系的文章 -

https://towardsdatascience.com/causal-discovery-does-the-cockerel-crowing-cause-the-sun-to-rise-f4308453ecfa

下一步将使用我的DirectedAcyclicGraph类。为了使文章更简洁，我没有在文章中包含源代码，但是以下是完整源代码的链接，以防你想要运行代码 - https://gist.github.com/grahamharrison68/9733b0cd4db8e3e049d5be7fc17b7602。

如果你决定使用它，如果你喜欢它，为什么不考虑选择给我买杯咖啡呢？

这是我对数据中因果关系的建议 -

from dag_tools import DirectedAcyclicGraph

training_model_edges : list = [("age", "received_training"), ("age", "earnings"),
                               ("education_years", "received_training"), ("education_years", "earnings"),
                               ("received_benefits", "received_training"), ("received_benefits", "earnings"),
                               ("single", "received_training"), ("single", "earnings"),
                               ("university_degree", "received_training"), ("university_degree", "earnings"),
                               ("received_training", "earnings")]

training_model_pos : dict = {"received_training": [1,1], "age": [3, 5], "education_years": [5, 5], "received_benefits": [7, 5], "single": [9, 5], "university_degree": [11, 5], "earnings": [13, 1]}

training_model = DirectedAcyclicGraph(edges=training_model_edges)    
training_model.display_pgm_model(pos=training_model_pos)

DAG的解释如下 -

• received_training（即参加培训计划）对收入（即未来收入）产生因果影响。

• 所有其他特征都会对个体是否可能参加培训计划产生因果影响。

• 所有其他特征也会对未来收入产生因果影响。

例如，一个人的年龄“导致”他们是否参加培训，可能是因为更多年轻人想接受培训，年龄还“导致”收入，可能是因为更有经验的老年人可以赚更多钱。

这种模式非常常见。当统计学家进行随机对照试验（RCT）时，他们可能会对与主要影响混合的变量进行条件控制。

这意味着对于年龄，他们可以将观察结果分成年龄组，查看每组中处理与年龄之间的关系，然后对每组进行比例平均以估计真实的整体效果。

然而，这种方法存在一些问题。例如，你如何定义组的界限？如果关键影响是16至18岁的人，但边界已设置为16至30岁怎么办？如果40至45岁的人没有观察结果怎么办？

另一种方法不是观察，而是干预。我们可以简单地强制每个人都参加培训，然后我们将看到真正的影响。但是，如果观察结果是历史性的（如LaLonde数据），并且干预已经太晚了呢？或者如果这是有关吸烟或肥胖的研究呢？

为了证明我们的理论，主体不能被强制吸烟或变得肥胖！

这就是“do”操作符的作用。它听起来像是魔法，但实际上可以构建一个因果推断模型，可以在不必在现实世界中进行干预的情况下准确地模拟这些干预。

这将节省大量的时间和金钱，消除了需要为大量变量进行条件控制的随机对照试验，并使像吸烟和肥胖这样的道德和伦理问题在现实世界的研究中也能够进行研究。

“Do”操作符的魔力

让我们想象一下，如果不是观察一组参加培训和未参加培训的人，而是可以回到过去，干预而不是观察，并让他们都参加培训。

在这种情况下，DAG将如下所示 -

这实际上是魔法“do”操作符正在做的事情。如果你干预并执行𝑑𝑜(𝑡𝑟𝑒𝑎𝑡=1)，则实际上是在“擦除”所有因果关系的输入线，因为无论年龄、教育和其他特征如何影响参加培训的概率，都会发生。

DoWhy库在幕后模拟了这种干预。它可以通过使用do微积分的规则将𝑝(earnings|𝑑𝑜(received_training=1))转换为一组可以从数据中计算出来的观察规则来实现这一点，但我们不进行物理干预。

我特意在本文中略去了数学细节。有很多文章展示了数学内容，但很少有文章显示了带有Python代码的实际示例，因此这是本文的重点。

注意：如果你想运行代码，你将需要我的DirectedAcyclicGraph类，如果你尚未下载它，请转到https://gist.github.com/grahamharrison68/9733b0cd4db8e3e049d5be7fc17b7602，并不要忘记如果你喜欢它的话，考虑给我买杯咖啡！

以下是在数据上执行“do”操作符的完整源代码…

import numpy as np
import dowhy.api

variable_types = {'received_training': 'd',  'age': 'c', 'education_years': 'c', 'received_benefits': 'd', 'single': 'd', 'university_degree': 'd', 'earnings': 'c'}

np.random.seed(1) 

df_do = df_training.causal.do(x={"received_training": 1},
                              outcome="earnings",
                              dot_graph=training_model.gml_graph, 
                              variable_types=variable_types,
                              proceed_when_unidentifiable=True)

display(df_do.groupby("received_training")["earnings"].agg(["median","mean"]))
display(df_do)

代码的含义是什么？

在进入真正令人惊叹的结果之前，逐行浏览代码可能会很有用。

首先，导入dowhy.api会神奇地扩展pandas DataFrame，以便该类获得一个新的causal.do方法。

接下来，在numpy中设置随机种子确保do方法的结果是可复制的。DoWhy文档没有提到任何关于设置随机种子的信息，这是通过纯粹的试错发现的。还要注意的是，随机种子需要在每次调用causal.do之前的上一条语句中设置，而不仅仅是在第一次调用之前。

causal.do的下一个神秘之处是variable_types参数。DoWhy文档不完整且不一致。尝试了许多不同的方法后，得出了以下结论 -

与文档所说的不同，只有两种类型是重要的 - “d”表示离散和“c”表示连续。

• 在统计学中，整数是离散的，但如果将整数声明为“d”，DoWhy会产生一些非常奇怪的结果。

• 根据对DoWhy文档和示例的深入研究，我的结论是整数需要声明为“c”表示连续。

• 在DoWhy源代码中，有一个称为infer_variable_types的方法，但它被存根化，没有代码，所以我编写了自己的实现，可以在DirectedAcyclicGraph.infer_variable_types()中找到这个静态方法。

这是极其重要的causal.do方法参数的含义 -

• x={"received_training": 1} 表示我们要“执行”的操作。在这种情况下，我们想知道如果每个人都被迫接受培训会发生什么，这在数据中通过received_training=1表示。

• outcome="earnings" - 这是我们正在寻找的结果或影响，即“执行”received_training=1对个体收入的影响是什么？

• dot_graph=training_model.gml_graph 通知do操作符在数据中存在的因果关系。training_model是我DirectedAcyclicGraph类的一个实例，我给它添加了一个以gml格式呈现结构的属性

• do方法需要传递common_causes或dot_graph来描述因果关系。

• dot_graph参数将接受dot或gml格式的结构，但在文档中没有提到；在我看来，gml更好，因为在DoWhy的其他地方都使用了gml。

• 指定图要比将common_causes设置为图好得多，因为图可以捕捉任何类型的结构，而common_causes则限制更多。同样，在DoWhy文档中没有提到这一点。

• variable_types参数已经解释过了。

• proceed_when_unidentifiable=True避免了干扰计算的讨厌用户提示。

它是如何工作的？

causal.do方法返回一个新的DataFrame，实际上模拟了强制干预，并提供了在每个人都接受培训的情况下可能收集到的数据 -

df_do["received_training"].value_counts()

与其他大多数Python因果关系库不同，DoWhy在这方面是不同的，因为大多数其他库只返回一个数字而不是DataFrame。

返回DataFrame最初可能会有点困惑，但稍微深入一点，它是一种强大、灵活和信息丰富的方法。

在不知道内部实现细节的情况下，我的结论是DoWhy通过根据需要用于“去混杂”本文中前面描述的混合效应的组进行数据采样来模拟随机对照试验（RCT）。

例如，比较原始观察数据和新的干预数据中的以下特征 -

def plot_do_comparison(feature : str, normalize : bool = True):

    fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(15, 6))

    df_training[received_training_filter][feature].value_counts(normalize=normalize).sort_index().plot(ax=axes[0], kind="bar", title="Received Training - Observation", xlabel=feature, ylim=(0,1) if normalize else False)
    df_do[feature].value_counts(normalize=normalize).sort_index().plot(ax=axes[1], kind="bar", title="Received Training - Intervention", xlabel=feature, ylim=(0,1) if normalize else False)

    plt.show()

plot_do_comparison("received_benefits", normalize=True)

显然，DoWhy已经以一种非常不同的方式重新对干预数据进行了采样。

培训对收入的真正“去混杂”影响是什么？

现在，我们只需通过查看df_do DataFrame来解释培训对收入的真正影响 -

print("Median and mean earnings of those receiving training from the interventional data")
display(df_do.groupby("received_training")["earnings"].agg(["median","mean"]))

print("Median and mean earnings of those receiving training from the original observational data")
display(df_training[received_training_filter].groupby("received_training")["earnings"].agg(["median","mean"]))

结论

使用概率在观察数据上的传统方法表明，那些参加培训的人实际上会比没有参加培训的人工资更低。

从观察数据中培训的平均工资为7,392。

在应用因果推断方法后，不建议取消培训计划，而是建议扩大培训计划，因为它为需要帮助提高长期收入的群体提供了更公平的机会。

我在文章末尾承诺了一些令人惊奇的事情，如果这个结果对你产生的启发性影响与对我产生的影响相同，我希望它能够兑现这个承诺。

每当数据中存在因果效应时，传统的预测方法可能会导致错误的结论和建议，这使得因果推断成为所有数据科学家必备的工具。

✄-----------------------------------------------

看到这里，说明你喜欢这篇文章，请点击「在看」或顺手「转发」「点赞」。

欢迎微信搜索「panchuangxx」，添加小编磐小小仙微信，每日朋友圈更新一篇高质量推文（无广告），为您提供更多精彩内容。

▼  ▼   扫描二维码添加小编  ▼  ▼