因果推断杂记——因果推断与线性回归、SHAP值理论的关系（十九）_causalforestdml

---

1 因果推断与线性回归的关系

第一个问题也是从知乎的这个问题开始：
因果推断（causal inference）是回归（regression）问题的一种特例吗？

其中经济学大佬慧航提到过，回归只是工具，因果推断可以用，其他研究方向也可以用。

在此给出我的看法，

因果推断，是需要考虑干预得（Y|X,T），其中干预效应是主要的差异点；
而一般的多元，只是(Y|X)，并没有考量到干预T的影响

1.1 DML的启发

所以，之前在做DML的时候，可以看到整个异质性HTE的求解经过：
因果推断笔记——DML ：Double Machine Learning案例学习（十六）

我们首先基于X使用ML获得T的残差和Y的残差，之后使用lr拟合残差，不同的是，这次我们把X和T的交互项加进来，即

$Y_i-M_y(X_i)=\tau (X_i)\cdot (T_i-M_t(X_i))+{ϵ}_i$

$\widetilde{Y_i}=\alpha +{\beta }_1\widetilde{T_i}+{\beta \text{β}\beta }_2{X}_i\text{Xi}{X}_i\widetilde{T_i}+{ϵ}_i$

然后我们就可以计算CATE的值了：

$\hat{\mu }(\partial Sale{s}_i,X_i)=M(Price=1,X_i)-M(Price=0,X_i)$

其中，M即最后的lr模型。

从以上DML求解无偏异质性CATE的过程看到，如果要得到无偏解，是需要经过一些求解步骤的；
关于残差正交化可得到无偏差因果效应的数学原理：https://zhuanlan.zhihu.com/p/41993542

1.2 特殊的离散回归 = 因果？

当然，这里感觉有个特例，$（Y|X,T）$ 中 如果不考虑任何协变量的影响，只有$（Y|T）$
那么此时，因果关系的ATE，应该就是等于$（Y|T）$ 离散回归的系数

---

2 因果推断中的ITE 与SHAP值理论的思考

本问题是由 多篇顶会看个体因果推断（ITE）的前世今生
和
机器学习模型可解释性进行到底 —— SHAP值理论（一）
引发的思考。

2.1 一些奇思妙想

本节可不参考，比较胡言乱语…

ITE代表的是无偏个体效应

再来看一下SHAP值中，可以“量化”不同特征，对个体的影响值，那么这个值，可以认为是RM的ITE吗？
虽然，SHAP值肯定是有偏的，但是也想沿着这个问题来看，SHAP值理论中的SHAP代表的怎么样的 “ITE”？在有偏的结论下，该如何解读？
之后简称sITE （此处应该需要公式推导，笔者水平就解读有限了）
个人理解:
$sITE=Predict(Y|X)-Predict(Y|X)的均值$

那么这里的实验组 - 对照组中的对照组就是，模型预测情况下，所有个体的“平均水平”
如果其中有一个特征是，是否有优惠券，

• 特征SHAP值>0，就代表，优惠券对其的刺激，比大家反应要强烈一些，更能刺激购买
• 特征SHAP值<0，代表，优惠券对其的刺激，要弱于常人，不利购买，不建议推送优惠券

沿着这个解读，给一个当下 “不负责任” 的结论： 值有偏，正负方向无偏：

• SHAP值量化出来的sITE 是有偏的，具体的值不具有参考意义；
• 但方向（正负号）代表整体趋势，还是可以借用的。

所以，不知道看到这里的看客，
有木有人，想用SHAP值来直接做“个性化推荐”的？

2.2 因果推断 -> shap值

经过一夜思考，有了一些Update，如下参考：案例：EconML + SHAP丰富可解释性

也就是，如果模型训练集中存在合理的实验组与对照组，且按照某些规范进行组合，那么就形同uplift model中的s-learner的形态（参考：思考：差分响应模型升级版（One-Model Approach）），此时shap值即为无偏，
此时eITE 含义应该就是 -> 某个样本 与 总体均值的处理效应ITE

2.3 ecoml中的shap值的案例

类似于如何用SHAP解释黑箱预测机器学习模型，我们也可以解释黑箱效应异质性模型。这种方法解释了为什么异质因果效应模型对特定人群产生了较大或较小的效应值。是哪些特征导致了这种差异?

当模型被简洁地描述时，这个问题很容易解决，例如线性异质性模型的情况，人们可以简单地研究模型的系数。
然而，当人们开始使用更具表现力的模型(如随机森林和因果森林)来建模效果异质性时，就会变得困难。
SHAP值对于理解模型从训练数据中提取的影响异质性的主导因素有很大的帮助。
我们的软件包提供了与SHAP库的无缝集成。每个CateEstimator都有一个方法shap_values，它返回每个处理和结果对的估计器输出的SHAP值解释。

然后，可以使用SHAP库提供的大量可视化功能对这些值进行可视化。
此外，只要有可能，我们的库就会为每种最终模型类型从SHAP库中调用快速专用算法，这可以大大减少计算时间。

2.3.1 单干预 - 单输出

## Ignore warnings
from econml.dml import CausalForestDML, LinearDML, NonParamDML
from econml.dr import DRLearner
from econml.metalearners import DomainAdaptationLearner, XLearner
from econml.iv.dr import LinearIntentToTreatDRIV
import numpy as np
import scipy.special
import matplotlib.pyplot as plt
import shap
from sklearn.ensemble import RandomForestRegressor, RandomForestClassifier
from sklearn.linear_model import Lasso

np.random.seed(123)
n_samples = 5000
n_features = 10
true_te = lambda X: (X[:, 0]>0) * X[:, 0]
X = np.random.normal(0, 1, size=(n_samples, n_features))
W = np.random.normal(0, 1, size=(n_samples, n_features))
T = np.random.binomial(1, scipy.special.expit(X[:, 0]))
y = true_te(X) * T + 5.0 * X[:, 0] + np.random.normal(0, .1, size=(n_samples,))
X_test = X[:min(100, n_samples)].copy()
X_test[:, 0] = np.linspace(np.percentile(X[:, 0], 1), np.percentile(X[:, 0], 99), min(100, n_samples))


# 因果树估计器
est = CausalForestDML(random_state=123)
est.fit(y, T, X=X, W=W)

# 线性估计器
est = LinearDML(random_state=123)
est.fit(y, T, X=X, W=W)

# 非参数
est = NonParamDML(model_y=RandomForestRegressor(min_samples_leaf=20, random_state=123),
                  model_t=RandomForestRegressor(min_samples_leaf=20, random_state=123),
                  model_final=RandomForestRegressor(min_samples_leaf=20, random_state=123),
                  random_state=123)
est.fit(y.ravel(), T.ravel(), X=X, W=W)

# 双重鲁棒学习
est = DRLearner(model_regression=RandomForestRegressor(min_samples_leaf=20, random_state=123),
                model_propensity=RandomForestClassifier(min_samples_leaf=20, random_state=123),
                model_final=RandomForestRegressor(min_samples_leaf=20, random_state=123),
                random_state=123)
est.fit(y.ravel(), T.ravel(), X=X, W=W)

# 元学习 DomainAdaptationLearner
est = DomainAdaptationLearner(models=RandomForestRegressor(min_samples_leaf=20, random_state=123),
                              final_models=RandomForestRegressor(min_samples_leaf=20, random_state=123),
                              propensity_model=RandomForestClassifier(min_samples_leaf=20, random_state=123))
est.fit(y.ravel(), T.ravel(), X=X)

# Xlearner 元学习
# Xlearner.shap_values uses a slow shap exact explainer, as there is no well defined final model
# for the XLearner method.
est = XLearner(models=RandomForestRegressor(min_samples_leaf=20, random_state=123),
               cate_models=RandomForestRegressor(min_samples_leaf=20, random_state=123),
               propensity_model=RandomForestClassifier(min_samples_leaf=20, random_state=123))
est.fit(y.ravel(), T.ravel(), X=X)

# 工具变量 
est = LinearIntentToTreatDRIV(model_y_xw=RandomForestRegressor(min_samples_leaf=20, random_state=123),
                              model_t_xwz=RandomForestClassifier(min_samples_leaf=20, random_state=123),
                              flexible_model_effect=RandomForestRegressor(min_samples_leaf=20, random_state=123),
                              random_state=123)
est.fit(y.ravel(), T.ravel(), Z=T.ravel(), X=X, W=W)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67

输出一个图来看：

est = CausalForestDML(random_state=123)
est.fit(y, T, X=X, W=W)
shap_values = est.shap_values(X[:20])
shap.plots.beeswarm(shap_values['Y0']['T0'])

1
2
3
4
5

详细的SHAP可参考：机器学习模型可解释性进行到底 —— SHAP值理论（一）

特征密度散点图：beeswarm

下图中每一行代表一个特征，横坐标为Shap值。特征的排序是按照shap 的平均绝对值，对模型来说的最重要特征。宽的地方表示有大量的样本聚集。

一个点代表一个样本，颜色越红说明特征本身数值越大，颜色越蓝说明特征本身数值越小。

可以看做一种特征重要性的排列图，从这里可以看到X0重要性高（排位）；
同时，X0值越大，对模型的效果影响越好（SHAP值为正）

2.3.2 多干预 - 多输出

# 数据生成
np.random.seed(123)
n_samples = 5000
n_features = 10
n_treatments = 2
n_outputs = 3
true_te = lambda X: np.hstack([(X[:, [0]]>0) * X[:, [0]],
                               np.ones((X.shape[0], n_treatments - 1))*np.arange(1, n_treatments).reshape(1, -1)])
X = np.random.normal(0, 1, size=(n_samples, n_features))
W = np.random.normal(0, 1, size=(n_samples, n_features))
T = np.random.normal(0, 1, size=(n_samples, n_treatments))
for t in range(n_treatments):
    T[:, t] = np.random.binomial(1, scipy.special.expit(X[:, 0]))
y = np.sum(true_te(X) * T, axis=1, keepdims=True) + 5.0 * X[:, [0]] + np.random.normal(0, .1, size=(n_samples, 1))
y = np.tile(y, (1, n_outputs))
for j in range(n_outputs):
    y[:, j] = (j + 1) * y[:, j]
X_test = X[:min(100, n_samples)].copy()
X_test[:, 0] = np.linspace(np.percentile(X[:, 0], 1), np.percentile(X[:, 0], 99), min(100, n_samples))

# 
est = CausalForestDML(n_estimators=400, random_state=123)
est.fit(y, T, X=X, W=W)

#
shap_values = est.shap_values(X[:200])

# 
plt.figure(figsize=(25, 15))
for j in range(n_treatments):
    for i in range(n_outputs):
        plt.subplot(n_treatments, n_outputs, i + j * n_outputs + 1)
        plt.title("Y{}, T{}".format(i, j))
        shap.plots.beeswarm(shap_values['Y' + str(i)]['T' + str(j)], plot_size=None, show=False)
plt.tight_layout()
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

这里可以看到有生成了三个outcome，
两个干预T，直接使用CausalForestDML因果树估计器，
所以就有6种情况可以得到：

---

3 [翻译]SHAP官方案例：回归与因果的差异

Be careful when interpreting predictive models in search of causal insights

3.1 笔者总结

这篇文献讲的是真啰嗦。。

本文主旨是在说，预测类模型，即使回归系数显著了，也仅代表关联，并不能表征因果关系，会有这样、那样的情况:

• 情况一：多个协变量中，可能互相存在相关性，从而混淆因果关系；比如本案例中 sales_Calls和interactions是高度依赖的，两者应该共同作用影响复订率，但是交互因子却比电话销售高很多
• 情况二：回归系数显著并不代表直接因果，广告支出由其他两个换算而成，其并不是直接影响复订率，是通过其他指标间接达到的，所以广告归因并非直接原因
• 情况三：树模型中的L1正则化会剔除低贡献特征，所以探究因果时候，模型的剔除可能会丢失一些因果线索

当你能够穷举混杂因子的时候，DML无疑是比较好的选择。

在混杂因子无法穷举的情况下，以下几种方法也是有使用场景的：

• IV：在我们不能随机分配治疗方案的情况下，工具变量技术可以用来识别因果影响，但我们可以随机引导一些客户接受治疗，比如发送电子邮件鼓励他们探索新产品功能。
• DID：当新治疗方法在不同群体间错开引入时，差异中的差异方法可能会有所帮助。
• 断点回归：当治疗模式显示出明显的界限时(例如，根据特定的、可衡量的特征(如每月收入超过5000美元)，断点回归是一个很好的选择。

3.2 产品重复订阅项目背景

像XGBoost这样的预测机器学习模型与SHAP这样的可解释性工具配对时变得更加强大。这些工具确定了输入特征和预测结果之间最有信息的关系，这对于解释模型正在做什么、让涉众参与进来以及诊断潜在的问题非常有用。

我们很容易进一步进行分析，假设解释工具也可以确定决策者在未来想要改变结果时应该操纵哪些特征。
在本文中，我们将讨论预测模型容易产生策略选择、智能决策上的误导。

原因与相关性与因果性的根本区别有关。SHAP可以更好捕捉预测ML模型获得的相关性。

但是，让相关性透明并不意味着它们就有因果关系!
所有预测模型都隐含地假设每个人在未来都会保持同样的行为方式，因此相关模式将保持不变。

为了理解如果某人的行为开始不同会发生什么，我们需要建立因果模型，这需要做出假设并使用因果分析工具。

假设我们的任务是构建一个模型，预测客户是否会复订他们的产品订阅。让我们假设，经过一些挖掘，
我们能够获得8个对预测客户流失很重要的功能:客户折扣、广告支出、客户每月使用情况、最后一次升级、客户报告的漏洞、与客户的互动、
与客户的销售电话以及宏观经济活动。

然后，我们使用这些功能来训练一个基本的XGBoost模型来预测客户是否会在订阅到期后续订:

# This cell defines the functions we use to generate the data in our scenario

import numpy as np
import pandas as pd
import scipy.stats
import sklearn
import xgboost

class FixableDataFrame(pd.DataFrame):
    """ Helper class for manipulating generative models.
    """
    def __init__(self, *args, fixed={}, **kwargs):
        self.__dict__["__fixed_var_dictionary"] = fixed
        super(FixableDataFrame, self).__init__(*args, **kwargs)
    def __setitem__(self, key, value):
        out = super(FixableDataFrame, self).__setitem__(key, value)
        if isinstance(key, str) and key in self.__dict__["__fixed_var_dictionary"]:
            out = super(FixableDataFrame, self).__setitem__(key, self.__dict__["__fixed_var_dictionary"][key])
        return out

# generate the data
def generator(n, fixed={}, seed=0):
    """ The generative model for our subscriber retention example.
    """
    if seed is not None:
        np.random.seed(seed)
    X = FixableDataFrame(fixed=fixed)

    # the number of sales calls made to this customer
    X["Sales calls"] = np.random.uniform(0, 4, size=(n,)).round()

    # the number of sales calls made to this customer
    X["Interactions"] = X["Sales calls"] + np.random.poisson(0.2, size=(n,))

    # the health of the regional economy this customer is a part of
    X["Economy"] = np.random.uniform(0, 1, size=(n,))

    # the time since the last product upgrade when this customer came up for renewal
    X["Last upgrade"] = np.random.uniform(0, 20, size=(n,))

    # how much the user perceives that they need the product
    X["Product need"] = (X["Sales calls"] * 0.1 + np.random.normal(0, 1, size=(n,)))

    # the fractional discount offered to this customer upon renewal
    X["Discount"] = ((1-scipy.special.expit(X["Product need"])) * 0.5 + 0.5 * np.random.uniform(0, 1, size=(n,))) / 2

    # What percent of the days in the last period was the user actively using the product
    X["Monthly usage"] = scipy.special.expit(X["Product need"] * 0.3 + np.random.normal(0, 1, size=(n,)))

    # how much ad money we spent per user targeted at this user (or a group this user is in)
    X["Ad spend"] = X["Monthly usage"] * np.random.uniform(0.99, 0.9, size=(n,)) + (X["Last upgrade"] < 1) + (X["Last upgrade"] < 2)

    # how many bugs did this user encounter in the since their last renewal
    X["Bugs faced"] = np.array([np.random.poisson(v*2) for v in X["Monthly usage"]])

    # how many bugs did the user report?
    X["Bugs reported"] = (X["Bugs faced"] * scipy.special.expit(X["Product need"])).round()

    # did the user renew?
    X["Did renew"] = scipy.special.expit(7 * (
          0.18 * X["Product need"] \
        + 0.08 * X["Monthly usage"] \
        + 0.1 * X["Economy"] \
        + 0.05 * X["Discount"] \
        + 0.05 * np.random.normal(0, 1, size=(n,)) \
        + 0.05 * (1 - X['Bugs faced'] / 20) \
        + 0.005 * X["Sales calls"] \
        + 0.015 * X["Interactions"] \
        + 0.1 / (X["Last upgrade"]/4 + 0.25)
        + X["Ad spend"] * 0.0 - 0.45
    ))

    # in real life we would make a random draw to get either 0 or 1 for if the
    # customer did or did not renew. but here we leave the label as the probability
    # so that we can get less noise in our plots. Uncomment this line to get
    # noiser causal effect lines but the same basic results
    X["Did renew"] = scipy.stats.bernoulli.rvs(X["Did renew"])

    return X

def user_retention_dataset():
    """ The observed data for model training.
    """
    n = 10000
    X_full = generator(n)
    y = X_full["Did renew"]
    X = X_full.drop(["Did renew", "Product need", "Bugs faced"], axis=1)
    return X, y

def fit_xgboost(X, y):
    """ Train an XGBoost model with early stopping.
    """
    X_train,X_test,y_train,y_test = sklearn.model_selection.train_test_split(X, y)
    dtrain = xgboost.DMatrix(X_train, label=y_train)
    dtest = xgboost.DMatrix(X_test, label=y_test)
    model = xgboost.train(
        { "eta": 0.001, "subsample": 0.5, "max_depth": 2, "objective": "reg:logistic"}, dtrain, num_boost_round=200000,
        evals=((dtest, "test"),), early_stopping_rounds=20, verbose_eval=False
    )
    return model

X, y = user_retention_dataset()
model = fit_xgboost(X, y)


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105

以上是数据加工过程，需要明确，这里所有特征加工，都是自己创造，所以知道某些特征是由谁产生

import shap

explainer = shap.Explainer(model)
shap_values = explainer(X)

# 层次聚类 + SHAP值
clust = shap.utils.hclust(X, y, linkage="single")
shap.plots.bar(shap_values, clustering=clust, clustering_cutoff=1)
1
2
3
4
5
6
7
8
9

这个条形图显示，提供的折扣、广告支出和报告的漏洞数量是驱动模型预测客户留存率的三大因素。
这很有趣，乍一看也很合理。条形图还包括我们稍后将使用的特征冗余聚类。
把层次聚类与SHAP结合,是一个亮点;那么这个组合图就可以表示：
- 特征重要性
- 特征交互关系/紧密程度（越底层的合并代表两个变量的相关）

然而，当我们深入挖掘并观察改变每个特征的值如何影响模型的预测时，我们发现了一些非直觉的模式。

SHAP散点图显示改变特征值如何影响模型对复订概率的预测。

如果蓝点遵循递增的模式，这意味着特征越大，模型预测的复订概率越高。

散布图显示了一些令人惊讶的发现:

• 报告更多bug的用户更有可能复订
• 折扣越大，用户续订的可能性就越小!

我们反复检查代码和数据管道，以排除漏洞，然后与一些业务合作伙伴交谈，他们给出了直观的解释:

• 高使用率、重视产品的用户更有可能报告漏洞并复订订阅。
• 销售人员往往会给他们认为不太可能对产品感兴趣的客户很高的折扣，这些客户的流失率更高。

模型中这些最初的反直觉关系是一个问题吗?那要看我们的目标是什么了!

我们对这个模型的最初目标是预测客户保留率，这对于财务规划中估算未来收入等项目非常有用。

因为报告更多bug的用户实际上更有可能复订，所以在模型中捕获这种关系有助于预测。
只要我们的模型具有良好的样本外拟合，我们就应该能够为金融提供良好的预测，因此不应该担心模型中这种关系。

# 不同特征SHAP值分布可视化
shap.plots.scatter(shap_values, ylabel="SHAP value\n(higher means more likely to renew)")

#shap.plots.scatter(shap_values,color=shap_values[:,"DIS"] ,ylabel="SHAP value\n(higher means more likely to renew)")

1
2
3
4
5

这是一类任务的一个例子，称为预测任务。在预测任务中，目标是预测给定一组特征X的结果Y(例如renew)。预测练习的一个关键组成部分是，我们只关心预测模型(X)在类似于我们的训练集的数据分布中接近Y。X和Y之间的简单关联可能有助于这些类型的预测。

但是，假设另一个团队使用我们的预测模型，其新目标是确定我们公司可以采取什么行动来留住更多的客户。这个团队非常关心每个X特性与Y之间的关系，这不仅仅是在我们的训练分布中，而且是在世界变化时产生的反事实场景中。在那个用例中，识别变量之间的稳定相关性不再足够;这个团队想知道操纵特性X是否会导致y的改变。当你告诉首席工程师你想让他引入新的bug来增加客户复订时，想象一下他的表情吧!

这是一类任务的一个例子叫做因果任务。在一个因果任务中，我们想知道改变世界的一个方面X(例如报告的bug)如何影响结果Y(复订)。在这种情况下，关键是要知道改变X是否会导致Y的增加，或者数据中的关系是否仅仅是相关的。

3.3 所有特征加工过程

理解因果关系的一个有用工具是写下我们感兴趣的数据生成过程的因果图。
我们例子中的因果图说明了为什么XGBoost客户留存率模型获得的稳健预测关系不同于希望计划干预措施以增加留存率的团队感兴趣的因果关系。

这张图只是对真实数据生成机制(上面定义的)的一个总结。
实心椭圆表示我们观察到的特征，而虚线椭圆表示我们没有测量到的隐藏特征。
每个特征都是带有箭头的所有特征的函数，再加上一些随机效果。

在我们的例子中，我们知道因果图，因为我们模拟了数据。
在实践中，真正的因果图是未知的，但是我们可以使用上下文特定的领域知识来推断哪些关系可以存在，哪些关系不可以存在。

import graphviz

import os
os.environ["PATH"] += os.pathsep + 'C:\\Users\\Desktop\\Graphviz\\bin\\'
import pydotplus 


names = [
         "Bugs reported", "Monthly usage", "Sales calls", "Economy",
         "Discount", "Last upgrade", "Ad spend", "Interactions"
]
g = graphviz.Digraph()
for name in names:
    g.node(name, fontsize="10")
g.node("Product need", style="dashed", fontsize="10")
g.node("Bugs faced", style="dashed", fontsize="10")
g.node("Did renew", style="filled", fontsize="10")

g.edge("Product need", "Did renew")
g.edge("Product need", "Discount")
g.edge("Product need", "Bugs reported")
g.edge("Product need", "Monthly usage")
g.edge("Discount", "Did renew")
g.edge("Monthly usage", "Bugs faced")
g.edge("Monthly usage", "Did renew")
g.edge("Monthly usage", "Ad spend")
g.edge("Economy", "Did renew")
g.edge("Sales calls", "Did renew")
g.edge("Sales calls", "Product need")
g.edge("Sales calls", "Interactions")
g.edge("Interactions", "Did renew")
g.edge("Bugs faced", "Did renew")
g.edge("Bugs faced", "Bugs reported")
g.edge("Last upgrade", "Did renew")
g.edge("Last upgrade", "Ad spend")
g



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

在这张图中有很多关系，但首先要关注的是，我们可以度量的一些特性会受到未度量的混杂特性(如产品需求和所面临的bug)的影响。

例如，报告更多错误的用户会遇到更多错误，因为他们使用产品的次数更多，而且他们也更有可能报告这些错误，因为他们更需要产品。

产品需求对更新有其直接的因果效应。因为我们不能直接测量产品需求，所以我们最终在预测模型中捕捉到的bug报告和更新之间的相关性结合了bug所带来的小的直接负面影响和产品需求带来的大的正面混杂效应。

下图将我们示例中的SHAP值与每个特征的真正因果效应进行了对比(由于我们生成了数据，因此在本示例中已知)。

'''
marginal_effects
'''
def marginal_effects(generative_model, num_samples=100, columns=None,\
                     max_points=20, logit=True, seed=0):
    """ Helper function to compute the true marginal causal effects.
    """
    X = generative_model(num_samples)
    if columns is None:
        columns = X.columns
    ys = [[] for _ in columns]
    xs = [X[c].values for c in columns]
    xs = np.sort(xs, axis=1)
    xs = [xs[i] for i in range(len(xs))]
    for i,c in enumerate(columns):
        xs[i] = np.unique([np.nanpercentile(xs[i], v, interpolation='nearest') for v in np.linspace(0, 100, max_points)])
        for x in xs[i]:
            Xnew = generative_model(num_samples, fixed={c: x}, seed=seed)
            val = Xnew["Did renew"].mean()
            if logit:
                val = scipy.special.logit(val)
            ys[i].append(val)
        ys[i] = np.array(ys[i])
    ys = [ys[i] - ys[i].mean() for i in range(len(ys))]
    return list(zip(xs, ys))

shap.plots.scatter(shap_values, ylabel="SHAP value\n(higher means more likely to renew)", overlay={
    "True causal effects": marginal_effects(generator, 10000, X.columns)
})

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

预测模型捕获了报告的bug对留存率的总体积极影响(如SHAP所示)，即使报告bug的因果影响为零，而遇到bug的影响为正。（这里原文是说负向的，但是从结果来看是正面的）

我们在折扣中也看到了类似的问题，折扣也是由未观察到的客户对产品的需求驱动的。我们的预测模型发现，折扣和留存率之间存在负相关关系，这是由这种未观察到的特性，即产品需求驱动的，尽管实际上在续订上存在一个小的正向因果效应!换句话说，如果两个客户有相同的产品需求，在其他方面是相似的，那么折扣较大的客户更有可能续订。

当我们开始把预测模型当作因果关系来解释时，这个情节还揭示了第二个更狡猾的问题。

注意，广告支出也有类似的问题——它对用户留存率没有因果关系(黑线是平的)，但预测模型却获得了积极的影响!

在这种情况下，广告支出只受上次升级和每月使用量的驱动，所以我们不存在不可观察的混淆问题，而是存在 可观察的混淆问题。
广告支出和影响广告支出的功能之间存在统计冗余。
当我们从几个特征中获得相同的信息时，预测模型可以使用这些特征中的任何一个进行预测，尽管它们并不都是因果关系。

虽然广告支出对复订本身没有因果关系，但它与推动更新的若干功能有着密切关系。我们的预测类模型将广告支出定义为一个有用的预测器，因为它涵盖了多个因果驱动因素为一体(从而形成一个稀疏的模型)，但如果我们开始将其解释为因果效应，这就会产生严重的误导。

现在，我们将依次处理示例中的每个部分，以说明何时预测模型可以准确地测量因果效应，何时不能。
我们还将介绍一些因果工具，这些工具有时可以在预测模型失败的情况下估计因果效应。

以广告支出作为举例，来说明，预测类模型，广告支出与复订率有很大关系，但并非一定是因果关系，
这里会给出一些案例证明

3.4 因果推断：DML

When predictive models can answer causal questions

注意，我们的预测模型在捕捉经济特征的真正因果效应方面做得很好(更好的经济对留存率有积极的影响)。那么，什么时候预测模型才能捕捉真正的因果效应呢?

让XGBoost能够为经济现象获得良好因果效应估计的重要因素是该特性的强大独立组件(在这个模拟中);
它对留存率的预测能力与任何其他可测量的功能或任何未测量的混杂因素都没有明显的冗余。因此，它不会受到来自未测量混杂物或特征冗余的偏见。

由于我们在SHAP柱状图的右侧添加了聚类，所以我们可以将数据特征的结构看作一个树状图。
当特征在树状图的底部(左侧)合并在一起时，这意味着这些特征包含的关于结果(复订)的信息是非常冗余的（多重共线），模型可以使用其中的任何一个特征，比如这里的ad_spend和monthly_usage。当特征在树状图的顶部(右边)合并在一起时，这意味着它们所包含的关于结果的信息是相互独立的。

通过注意到Economy直到集群树状图的最顶端才与任何其他特征合并，其独立于所有其他测量的特征。
这告诉我们，经济不会受到观察混淆。
但要相信经济效应是有因果关系的，我们还需要检查未观察到混淆。检查未测量混杂因素比较困难，需要使用领域知识(由我们上面示例中的业务合作伙伴提供)。

对于经典的ML预测模型来说，要传递因果结果，特征不仅需要独立于模型中的其他特征，而且也需要独立于未观察混杂因子。兴趣驱动的例子自然地表现出这种程度的独立性并不常见，但当我们的数据包含一些实验时，我们经常可以找到独立特征的例子。

在大多数真实世界的数据集中，特征不是独立的和不混淆的，所以标准的预测模型不会了解真正的因果效应。
因此，用SHAP解释它们并不能揭示因果效应。
但这并非一无是处，有时我们可以使用观察因果推理工具来解决或至少将这个问题最小化。

因果推理可以帮助的第一个场景被观察到的混淆因子。
当有另一个特征同时影响原始特征和我们所预测的结果时，该特征就被“混淆”了。如果我们能测量其他特征，它就被称为可观测混杂因子observed confounder。

在我们的场景中，这方面的一个例子就是ad_spend广告支出。
尽管广告支出对留存率没有直接的因果关系，但它与“最后升级”和“每月使用”功能相关，这些功能确实会推动留存率。

我们的预测模型认为广告支出是留存率最好的单一预测因素之一，因为它通过相关性捕捉了许多真正的因果驱动因素。

正则对特征解释的弊端 ： 自动忽略了某些特征

XGBoost引入了正则化，这是一种花哨的说法，它试图选择预测良好+简单的概率模型。
如果它使用一个特征而不是三个特征能够同样准确地预测，它将倾向于这样做以避免过拟合。（如果是L1就会直接删除）

但这意味着，如果广告支出与最后升级和每月使用量高度相关，XGBoost可能只会使用广告支出这一个特征

XGBoost(或任何其他具有正则化的机器学习模型)的这种特性对于生成未来留存率的稳健预测非常有用，但对于理解如果我们想要提高留存率，我们应该操纵哪些特征却不是很好。

这突出了为每个问题匹配正确的建模工具的重要性。

与漏洞报告的例子不同，增加广告支出可以提高留存率的结论并没有什么直观的错误。如果没有适当地关注我们的预测模型是什么或不是什么，我们很可能会继续进行这一发现，并只是在增加广告支出后才认识到自己的错误，而未能获得我们预期的更新结果。

# Economy is independent of other measured features.
shap.plots.bar(shap_values, clustering=clust, clustering_cutoff=1)
1
2

3.5 Observational Causal Inference

对于广告支出来说，好消息是我们能够衡量所有可能混淆它的特征。因此，这是观察到的混淆的一个例子，我们应该能够只使用我们已经收集的数据来解开相关模式;我们只需要使用观察因果推理的正确工具。这些工具让我们能够明确哪些功能会干扰广告支出，然后根据这些功能进行调整，从而准确估计广告支出对产品更新的因果影响。

一个特别灵活的观察因果推理工具是双重/去偏机器学习。
它使用任何你想要的机器学习模型，首先解出感兴趣的特征(即广告花费)，然后估计改变该特征的平均因果效应(即因果效应的平均斜率)。

DML的工作原理如下（讲得真是绕来绕去）:

• 使用一组可能的混淆因素(即任何不是由Ad Spend引起的功能)训练模型来预测用户感兴趣的功能(不受ad_spend影响)。
• 训练一个模型来预测结果(DID Renew)，使用相同的一组可能的混杂因素。
• 使用感兴趣的因果特征的残差变异训练模型预测结果的残差变异(减去我们的预测后剩下的变异)。

我们的直觉是，如果Ad Spend能够引起复订，那么其他混淆功能所不能预测的Ad Spend部分应该与其他混淆功能所不能预测的复订部分相关联。

说另一种方式,DML假定有一个独立的(未被注意的)噪声特征,影响广告支出(因为广告支出并不完全取决于其他特性),所以我们可以归罪于这个独立的噪声特性的值,然后训练模型这个独立特性来预测输出。

虽然我们可以手动执行所有的DML步骤，但使用像econML或CausalML这样的因果推断包更容易。这里我们使用econML的LinearDML模型。

这将返回一个表示治疗是否具有非零因果效应的p值，并在我们的场景中很好地运行，正确地识别出没有证据表明在更新上的广告支出具有因果效应(p值= 0.85):

绕来绕去就是在说DML的作用以及使用场景

'''
联谊econml
'''
from econml.dml import LinearDML
from sklearn.base import BaseEstimator, clone
import matplotlib.pyplot as plt

class RegressionWrapper(BaseEstimator):
    """ Turns a classifier into a 'regressor'.

    We use the regression formulation of double ML, so we need to approximate the classifer
    as a regression model. This treats the probabilities as just quantitative value targets
    for least squares regression, but it turns out to be a reasonable approximation.
    """
    def __init__(self, clf):
        self.clf = clf

    def fit(self, X, y, **kwargs):
        self.clf_ = clone(self.clf)
        self.clf_.fit(X, y, **kwargs)
        return self

    def predict(self, X):
        return self.clf_.predict_proba(X)[:, 1]

# Run Double ML, controlling for all the other features
def double_ml(y, causal_feature, control_features):
    """ Use doubleML from econML to estimate the slope of the causal effect of a feature.
    """
    xgb_model = xgboost.XGBClassifier(objective="binary:logistic", random_state=42)
    est = LinearDML(model_y=RegressionWrapper(xgb_model))
    est.fit(y, causal_feature, W=control_features)
    return est.effect_inference()

def plot_effect(effect, xs, true_ys, ylim=None):
    """ Plot a double ML effect estimate from econML as a line.

    Note that the effect estimate from double ML is an average effect *slope* not a full
    function. So we arbitrarily draw the slope of the line as passing through the origin.
    """
    plt.figure(figsize=(5, 3))

    pred_xs = [xs.min(), xs.max()]
    mid = (xs.min() + xs.max())/2
    pred_ys = [effect.pred[0]*(xs.min() - mid), effect.pred[0]*(xs.max() - mid)]

    plt.plot(xs, true_ys - true_ys[0], label='True causal effect', color="black", linewidth=3)
    point_pred = effect.point_estimate * pred_xs
    pred_stderr = effect.stderr * np.abs(pred_xs)
    plt.plot(pred_xs, point_pred - point_pred[0], label='Double ML slope', color=shap.plots.colors.blue_rgb, linewidth=3)
    # 99.9% CI
    plt.fill_between(pred_xs, point_pred - point_pred[0] - 3.291 * pred_stderr,
                     point_pred - point_pred[0] + 3.291 * pred_stderr, alpha=.2, color=shap.plots.colors.blue_rgb)
    plt.legend()
    plt.xlabel("Ad spend", fontsize=13)
    plt.ylabel("Zero centered effect")
    if ylim is not None:
        plt.ylim(*ylim)
    plt.gca().xaxis.set_ticks_position('bottom')
    plt.gca().yaxis.set_ticks_position('left')
    plt.gca().spines['right'].set_visible(False)
    plt.gca().spines['top'].set_visible(False)
    plt.show()

# estimate the causal effect of Ad spend controlling for all the other features
causal_feature = "Ad spend"
control_features = [
    "Sales calls", "Interactions", "Economy", "Last upgrade", "Discount",
    "Monthly usage", "Bugs reported"
]
effect = double_ml(y, X[causal_feature], X.loc[:,control_features])

# plot the estimated slope against the true effect
xs, true_ys = marginal_effects(generator, 10000, X[["Ad spend"]], logit=False)[0]
plot_effect(effect, xs, true_ys, ylim=(-0.2, 0.2))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

记住，DML(或任何其他观察因果推理方法)只有在您能够测量和识别所有可能的混杂因素时才有效。

在这里我们知道因果图，可以看到每月使用率和最后升级是我们需要控制的两个直接混淆因素。

但如果我们不知道因果关系图，我们仍然可以查看SHAP柱状图中的冗余，并看到每月使用和最后升级是最冗余的功能，也是控制的好候选功能(折扣和bug报告)。

3.6 DML中处理混杂因子

非混淆冗余

因果推理可以帮助的第二个场景是非混淆冗余。
这种情况发生在我们想要因果效应的特性，或者是由模型中包含的另一个特性所驱动的特性，但这个特性并不是我们感兴趣的特性的混杂物。

绕来绕去，就是模型不仅有混杂因子W ，也有协变量X

这方面的一个例子就是Sales_Calls。销售电话直接影响留存率，但也通过互动间接影响留存率。

当我们在模型中同时包含Interactions和Sales_Calls时，这两个特征共享的因果效应将被迫在它们之间展开。

我们可以在上面的SHAP散点图中看到这一点，这表明XGBoost如何低估了Sales Calls的真正因果效应，因为大多数因果效应都被放到了interaction特性上。

原则上，协变量可以通过从模型中去除冗余变量来固定(见下文)。例如，如果我们从模型中删除了交互，那么我们将捕捉到进行销售电话更新概率的全部效果。

这种删除对于DML也很重要，因为如果控制由感兴趣的特征引起的下游特征，DML将无法捕获间接因果效应。

在这种情况下，DML将只测量不通过其他特征的“直接”效应。然而，DML对于控制上游非混杂冗余(冗余特征导致感兴趣的特征)来说是稳健的，尽管这将降低检测真实效果的统计能力。

不幸的是，我们通常不知道真正的因果图，所以很难知道另一个特征何时与我们感兴趣的特征是冗余的，因为观察到的混杂冗余与非混杂冗余。

如果这是因为混淆，那么我们应该使用像DML这样的方法来控制该特征，而如果它是一个下游结果，那么如果我们想要完整的因果效应而不仅仅是直接效应，那么我们应该将该特征从我们的模型中删除。

Controlling for a feature we shouldn’t tends to hide or split up causal effects, while failing to control for a feature we should have controlled for tends to infer causal effects that do not exist.

控制一个我们不应该倾向于隐藏或分割因果效应的功能，而未能控制一个我们应该控制的功能则倾向于推断根本不存在的因果效应。当你不确定的时候，这通常使得控制功能成为更安全的选择。

# Interactions and sales calls are very redundant with one another.
shap.plots.bar(shap_values, clustering=clust, clustering_cutoff=1)
1
2

# Fit, explain, and plot a univariate model with just Sales calls
# Note how this model does not have to split of credit between Sales calls and
# Interactions, so we get a better agreement with the true causal effect.
sales_calls_model = fit_xgboost(X[["Sales calls"]], y)
sales_calls_shap_values = shap.Explainer(sales_calls_model)(X[["Sales calls"]])
shap.plots.scatter(sales_calls_shap_values, overlay={
    "True causal effects": marginal_effects(generator, 10000, ["Sales calls"])
})
1
2
3
4
5
6
7
8

3.7 混杂因子无法穷尽的时候

(或任何其他假定无混淆的因果推理方法)只有在您能够测量和识别您想要估计因果影响的特征的所有可能混淆时才有效。
如果你不能测量所有的混杂因素，那么你将处于最困难的情况:未被观察的混杂因素。

折扣和bug报告功能都受到未观察到的混淆，因为不是所有重要的变量(例如，Product Need和Bugs Faced)都在数据中测量

。尽管这两个特征相对独立于模型中的所有其他特征，但仍有一些重要的驱动因素是无法测量的。
在这种情况下，需要观察混杂物的预测模型和因果模型，如DML，都将失败。这就是为什么在控制了所有其他观察特征的情况下，DML估计了Discount特征的巨大负因果效应:

def double_ml(y, causal_feature, control_features):
    """ Use doubleML from econML to estimate the slope of the causal effect of a feature.
    """
    xgb_model = xgboost.XGBClassifier(objective="binary:logistic", random_state=42)
    est = LinearDML(model_y=RegressionWrapper(xgb_model))
    est.fit(y, causal_feature, W=control_features)
    return est.effect_inference()


# estimate the causal effect of Ad spend controlling for all the other features
causal_feature = "Discount"
control_features = [
    "Sales calls", "Interactions", "Economy", "Last upgrade",
    "Monthly usage", "Ad spend", "Bugs reported"
]
effect = double_ml(y, X[causal_feature], X.loc[:,control_features])

# plot the estimated slope against the true effect
xs, true_ys = marginal_effects(generator, 10000, X[[causal_feature]], logit=False)[0]
plot_effect(effect, xs, true_ys, ylim=(-0.5, 0.2))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

除非能够测量之前未测量的特征(或与之相关的特征)，否则很难在未观察到的混杂存在的情况下找到因果效应。
在这些情况下，确定能够为政策提供信息的因果效应的唯一方法是创建或利用一些打破兴趣特征与未测量混杂物之间相关性的随机化。
在这种情况下，随机实验仍然是发现因果效应的黄金标准。

基于工具变量、DID或断点回归 的专门因果工具有时可以利用部分随机化，即使在不可能进行完整实验的情况下也是如此。

例如，

• IV：在我们不能随机分配治疗方案的情况下，工具变量技术可以用来识别因果影响，但我们可以随机引导一些客户接受治疗，比如发送电子邮件鼓励他们探索新产品功能。
• DID：当新治疗方法在不同群体间错开引入时，差异中的差异方法可能会有所帮助。
• 最后，当治疗模式显示出明显的界限时(例如，根据特定的、可衡量的特征(如每月收入超过5000美元)，断点回归是一个很好的选择。

3.8 小结

灵活的预测模型，如XGBoost或LightGBM是解决预测问题的强大工具。

然而，它们并不是内在的因果模型，所以用SHAP解释它们在许多常见情况下无法准确回答因果问题。除非模型中的特征是实验变化的结果，否则将SHAP应用于预测模型而不考虑混淆通常不是衡量用于告知政策的因果影响的适当工具。

SHAP和其他解释性工具对于因果推理可能很有用，而且SHAP被集成到许多因果推理包中，但是这些用例在本质上是明确的因果性的。

为此，使用我们为预测问题收集的相同数据，并使用像double ML这样的因果推理方法，这些方法是特别设计来返回因果结果的，通常是告知策略的好方法。在其他情况下，只有实验或其他随机化的来源才能真正回答“如果”的问题。

因果推理总是需要我们做出重要的假设。本文的主要观点是，我们通过将正常预测模型解释为因果关系而做出的假设通常是不现实的。