因果推断笔记——因果图建模之微软开源的EconML

Posted 2021-12-11 悟乙己

文章目录

• 1 EconML介绍
• • 1.1 EconML介绍
  • 1.2 一些理论解答
  • 1.3 常规CATE的估计器
  • 1.4 IV工具变量 + CATE的估计器
  • 1.5 动态处理效应的估计器
• 2 智能营销案例一：（econml）不同优惠折扣下的用户反应
• • 2.1 背景
  • 2.2 数据准备与理解
  • 2.3 生成假的Groud Truth的估计效应
  • 2.4 X~Y 分析：线性估计器LinearDML
  • 2.5 X~Y 分析：非线性估计器：因果树CausalForestDML
  • 2.6 X|T~Y分析：SingleTreeCateInterpreter哪些用户比较积极/消极
  • 2.7 T~X分析：SingleTreePolicyInterpreter 什么收入的人该打折
  • 2.8 干预后的结果对比（该打折的人打折后）
• 3 智能营销案例二：（econml+dowhy）不同优惠折扣下的用户反应
• • 3.1 背景
  • 3.2 数据准备与理解
  • 3.3 生成假的Groud Truth的估计效应
  • 3.4 利用dowhy创建因果图 + EconML创建线性估计器LinearDML
  • 3.5 EconML因果树估计器CausalForestDML + dowhy构建因果图
  • 3.6 dowhy基于因果树的估计稳定性检验：反驳
  • 3.7 T~X分析：SingleTreePolicyInterpreter 什么收入的人该打折
  • 3.8 干预后的结果对比（该打折的人打折后）
• 4 案例：EconML + SHAP丰富可解释性
• • 4.1 单干预 - 单输出
  • 4.2 多干预 - 多输出

---

1 EconML介绍

微软EconML简介：基于机器学习的Heterogeneous Treatment Effects估计

GITHUB ： microsoft/EconML
econml DOC：user guide

1.1 EconML介绍

机器学习最大的promise之一是在许多领域实现决策的自动化。
许多数据驱动的决策场景的核心问题是对heterogeneous treatment effects的估计，也即：对于具有特定特征集的样本，干预对输出结果的causal effect是什么？
简言之，这个Python工具包旨在：

• 测量某些干预变量T对结果变量Y的causal effect
• 控制一组特征X和W，来衡量causal effect如何随X的变化而变化

这个Python工具包旨在测量某些干预变量T对结果变量Y的causal effect，控制一组特征X和W，来衡量causal effect如何随X的变化而变化。

这个包里实现的方法适用于观察的数据集（非实验或历史）。为了使估计结果具有因果解释，一些算法假设数据中没有未观察到的混杂因子（即同时对T和Y产生影响的因子均被包含在X，W中），而其他一些算法则假设可以使用工具变量Z（即观测变量Z对干预T有影响，但对结果变量Y没有直接影响）。
并且包中的大多数算法都可以提供置信区间和推断结果。

我们将概述最近将机器学习与因果推理结合起来的方法，以及机器学习给因果推理估计方法带来的重要统计性能。
相比DoWhy，EconML借助一些更复杂的机器学习算法来进行因果推断。在EconML中可以使用的因果推断方法有：

• 双机器学习(Double Machine Learning)。
• 双重鲁棒学习(Doubly Robust Learner)。
• 树型学习器(Forest Learners)。
• 元学习器(Meta Learners)。
• 深度工具变量法(Deep IV).
• 正交随机树(Orthogonal Random Forest)
• 加权双机器学习(Weighted Double Machine Learning)

我们还将概述置信区间构建(例如自举、小袋自举、去偏lasso)、可解释性(形状值、树解释器)和策略学习(双鲁棒策略学习)的方法。
EconML希望做到以下几点：

• 实现同时涉及计量经济学和机器学习的最新方法
• 保持对effect heterogeneity建模的灵活性（通过随机森林、boosting、lasso和神经网络等技术），同时- 保留所学模型的因果解释，并同时提供有效的置信区间
• 使用统一的API
• 构建标准的Python包以用于机器学习和数据分析

---

1.2 一些理论解答

参考：
因果推断笔记—— 相关理论：Rubin Potential、Pearl、倾向性得分、与机器学习异同（二）的【3.0.2 章节】

异质处理效应 （Heterogenous treatment effects,HTE）
简单理解就是按某些规律进行分组之后的ATE，也就是“条件平均处理效应”(Conditional Average Treatment Effect，CATE)
每个子组的平均处理效应被称为“条件平均处理效应”(Conditional Average Treatment Effect，CATE) ，也就是说，每个子组的平均处理效应被称为条件平均处理效应，以子组内的分类方式为条件。

1.3 常规CATE的估计器

在econML中，“条件平均处理效应”(Conditional Average Treatment Effect，CATE) 的四种估计方式：

• 双机器学习(Double Machine Learning)。
• 双重鲁棒学习(Doubly Robust Learner)。
• 元学习器(Meta Learners)。
• 正交随机树(Orthogonal Random Forest)

DML的几种方法包括：

• econml.dml.DML(*, model_y, model_t, model_final) # 基础双重ML的模型
• econml.dml.LinearDML(*[, model_y, model_t, …]) # 线性估计器
• econml.dml.SparseLinearDML(*[, model_y, …]) # 稀疏线性估计器
• econml.dml.CausalForestDML(*[, model_y, …]) # 因果树
• econml.dml.NonParamDML(*, model_y, model_t, …) # 非参数ML估计器， that can have arbitrary final ML models of the CATE.

双重鲁棒学习(Doubly Robust Learner)几种方法：

• econml.dr.DRLearner(*[, model_propensity, …])
• econml.dr.LinearDRLearner(*[, …])
• econml.dr.SparseLinearDRLearner(*[, …])
• econml.dr.ForestDRLearner(*[, …])

CATE estimator that uses doubly-robust correction techniques to account for covariate shift (selection bias) between the treatment arms.

元学习器(Meta Learners)四种方法：

• econml.metalearners.XLearner(*, models[, …])
• econml.metalearners.TLearner(*, models[, …])
• econml.metalearners.SLearner(*, overall_model)
• econml.metalearners.DomainAdaptationLearner(*, …) 元算法，使用领域适应技术来解释

正交随机树(Orthogonal Random Forest) 的两种方法：

• econml.orf.DMLOrthoForest(*[, n_trees, …])
• econml.orf.DROrthoForest(*[, n_trees, …])

1.4 IV工具变量 + CATE的估计器

Double Machine Learning (DML) IV的几种方法：

• econml.iv.dml.OrthoIV(*[, model_y_xw, …])
  用IV实现正交/双ml方法进行CATE估计
• econml.iv.dml.DMLIV(*[, model_y_xw, …])
  参数DMLIV估计器的基类来估计一个CATE。
• econml.iv.dml.NonParamDMLIV(*[, model_y_xw, …])
  非参数DMLIV的基类，允许在DMLIV算法的最后阶段采用任意平方损耗的ML方法。

Doubly Robust (DR) IV 稳健+IV几种方法：

• econml.iv.dr.DRIV(*[, model_y_xw, …])
• econml.iv.dr.LinearDRIV(*[, model_y_xw, …])
• econml.iv.dr.SparseLinearDRIV(*[, …]) # 稀疏
• econml.iv.dr.ForestDRIV(*[, model_y_xw, …]) # 因果树
• econml.iv.dr.IntentToTreatDRIV(*[, …])
• econml.iv.dr.LinearIntentToTreatDRIV(*[, …]) # 为线性意图处理A/B测试设置实现DRIV算法

DeepIV方法：

• econml.iv.nnet.DeepIV(*, n_components, m, h, …)

Sieve Methods 的几种方法：

• econml.iv.sieve.SieveTSLS(*, t_featurizer, …)
• econml.iv.sieve.HermiteFeatures(degree[, …])
• econml.iv.sieve.DPolynomialFeatures([…])

1.5 动态处理效应的估计器

• econml.dynamic.dml.DynamicDML(*[, model_y, …])
  用于动态处理效果估计的CATE估计器。

2 智能营销案例一：（econml）不同优惠折扣下的用户反应

参考：Case Study - Customer Segmentation at An Online Media Company.ipynb

2.1 背景

如今，商业决策者依靠评估干预的因果效应来回答有关策略转变的假设问题，比如打折促销特定产品、向网站添加新功能或增加销售团队的投资。
然而，人们越来越感兴趣的是了解不同用户对这两种选择的不同反应，而不是学习是否要针对所有用户的特定干预采取行动。 识别对干预反应最强烈的用户特征，有助于制定规则，将未来用户划分为不同的群体。 这可以帮助优化政策，使用最少的资源，获得最大的利润。

在本案例研究中，我们将使用一个个性化定价示例来解释EconML库如何适应这个问题，并提供健壮可靠的因果解决方案。

目的是了解不同收入水平人群的需求的异质性价格弹性，了解哪个用户对小折扣的反应最强烈。
此外，他们的最终目标是确保，尽管对一些消费者降低了价格，但需求有足够的提高，以提高整体收入。

EconML的基于“DML”的估计器可用于获取现有数据中的折扣变化，以及一组丰富的用户特性，以估计随多个客户特性而变化的异构价格敏感性。

然后,EconML的可解释性的两个函数：

• SingleTreeCateInterpreter
  提供了一个presentation-ready总结的关键特性,解释折扣最大的响应能力的差异,
• SingleTreePolicyInterpreter
• 建议谁应该获得折扣政策为了增加收入(不仅需求),这可能有助于该公司在未来为这些用户设定一个最优价格。

2.2 数据准备与理解

# Some imports to get us started
# Utilities
import os
import urllib.request
import numpy as np
import pandas as pd

# Generic ML imports
from sklearn.preprocessing import PolynomialFeatures
from sklearn.ensemble import GradientBoostingRegressor

# EconML imports
from econml.dml import LinearDML, CausalForestDML
from econml.cate_interpreter import SingleTreeCateInterpreter, SingleTreePolicyInterpreter

import matplotlib.pyplot as plt

%matplotlib inline

# Import the sample pricing data
file_url = "https://msalicedatapublic.blob.core.windows.net/datasets/Pricing/pricing_sample.csv"
train_data = pd.read_csv(file_url)

# Define estimator inputs
Y = train_data["demand"]  # outcome of interest
T = train_data["price"]  # intervention, or treatment
X = train_data[["income"]]  # features
W = train_data.drop(columns=["demand", "price", "income"])  # confounders
# Get test data
X_test = np.linspace(0, 5, 100).reshape(-1, 1)
X_test_data = pd.DataFrame(X_test, columns=["income"])

其中train_data是10000*11完整的数据集，而X_test_data是新生成的一列自变量X，
这里可以看到，test数据集 其实不用捎上混杂因子的W变量们，也不用额外生成干预变量T

数据格式为：

Feature Name	Type	Details
account_age	W	user’s account age
age	W	user’s age
avg_hours	W	the average hours user was online per week in the past
days_visited	W	the average number of days user visited the website per week in the past
friend_count	W	number of friends user connected in the account
has_membership	W	whether the user had membership
is_US	W	whether the user accesses the website from the US
songs_purchased	W	the average songs user purchased per week in the past
income	X	user’s income
price	T	the price user was exposed during the discount season (baseline price * small discount)
demand	Y	songs user purchased during the discount season

数据集*有~ 10000个观察，包括9个连续和分类变量，代表用户的特征和在线行为历史，如年龄，日志收入，以前的购买，每周以前的在线时间等。
为了保护公司的隐私，我们这里以模拟数据为例。数据是综合生成的，特征分布与真实分布不一致。
然而，特征名称保留了它们的名称和含义.

这里来盘一下具体的因果逻辑：

• 其他变量：account_age ~ songs_purchased - W - 混杂因子
• income - X - 考察变量 - 用户收入
• demand - Y - outcome - 销量
• Price - T - 干预，折扣，取值为[1,0.9,0.8]，根据下面的公式的来

$$T=\{\begin{array}{llc}1& \text{with }p=0.2,& \\ 0.9& \text{with }p=0.3,& \text{if income}<1\\ 0.8& \text{with }p=0.5,& \\ & & \\ 1& \text{with }p=0.7,& \\ 0.9& \text{with }p=0.2,& \text{if income}\ge 1\\ 0.8& \text{with }p=0.1,& \end{array}$$

2.3 生成假的Groud Truth的估计效应

估计效应的对比，此时没有一个标准答案，所以这里就搞了一个手算版本，当作真实的估计量，来进行对比。
{ γ ( X ) = − 3 − 14 ⋅ { income < 1 } β ( X , W ) = 20 + 0.5 ⋅ avg-hours + 5 ⋅ { day-svisited > 4 } \\begin{cases} \\gamma(X) = -3 - 14 \\cdot \\{\\text{income}<1\\} \\\\ \\beta(X,W) = 20 + 0.5 \\cdot \\text{avg-hours} + 5 \\cdot \\{\\text{day-svisited}>4\\} \\\\ \\end{cases} {γ(X)=−3−14⋅{income<1}β(X,W)=20+0.5⋅avg-hours+5⋅{day-svisited>4}​
最终要算的是价格弹性elasticity
$$\begin{gathered} Y=\gamma (X)\cdot T+\beta (X,W) \\ elasticity(手算)=(\gamma (X)\cdot T)/Y \end{gathered}$$
当然，在实际计算中有取对数的操作，公式会变成：
$$\begin{gathered} log(Y)=\theta (X)\cdot log(T)+f(X,W)+ϵ \\ log(T)=g(X,W)+\eta \end{gathered}$$

整体生成手算的估计效应：

# Define underlying treatment effect function given DGP
def gamma_fn(X):
    return -3 - 14 * (X["income"] < 1)

def beta_fn(X):
    return 20 + 0.5 * (X["avg_hours"]) + 5 * (X["days_visited"] > 4)

def demand_fn(data, T):
    Y = gamma_fn(data) * T + beta_fn(data)
    return Y

def true_te(x, n, stats):
    if x < 1:
        subdata = train_data[train_data["income"] < 1].sample(n=n, replace=True)
    else:
        subdata = train_data[train_data["income"] >= 1].sample(n=n, replace=True)
    te_array = subdata["price"] * gamma_fn(subdata) / (subdata["demand"])
    if stats == "mean":
        return np.mean(te_array)
    elif stats == "median":
        return np.median(te_array)
    elif isinstance(stats, int):
        return np.percentile(te_array, stats)

# Get the estimate and range of true treatment effect
truth_te_estimate = np.apply_along_axis(true_te, 1, X_test, 1000, "mean")  # estimate
truth_te_upper = np.apply_along_axis(true_te, 1, X_test, 1000, 95)  # upper level
truth_te_lower = np.apply_along_axis(true_te, 1, X_test, 1000, 5)  # lower level

这里的truth_te_estimate就是最终估计数字了，公式即为:$$elasticity(手算)=(\gamma (X)\cdot T)/Y$$

  te_array = subdata["price"] * gamma_fn(subdata) / (subdata["demand"])

2.4 X~Y 分析：线性估计器LinearDML

# Get log_T and log_Y
log_T = np.log(T)
log_Y = np.log(Y)

# Train EconML model
est = LinearDML(
    model_y=GradientBoostingRegressor(),
    model_t=GradientBoostingRegressor(),
    featurizer=PolynomialFeatures(degree=2, include_bias=False),
)
est.fit(log_Y, log_T, X=X, W=W, inference="statsmodels")
# Get treatment effect and its confidence interval    得到治疗效果及其置信区间
te_pred = est.effect(X_test)
# 因为是配对实验，PSM，那就是 ps(y=1) - ps(y=0)，每个人都有，
# 后面求ATE的时候，会进行平均
te_pred_interval = est.effect_interval(X_test)   # 置信区间 上限与下限



# Compare the estimate and the truth
plt.figure(figsize=(10, 6))
plt.plot(X_test.flatten(), te_pred, label="Sales Elasticity Prediction")
plt.plot(X_test.flatten(), truth_te_estimate, "--", label="True Elasticity")
plt.fill_between(
    X_test.flatten(),
    te_pred_interval[0],
    te_pred_interval[1],
    alpha=0.2,
    label="95% Confidence Interval",
)
plt.fill_between(
    X_test.flatten(),
    truth_te_lower,
    truth_te_upper,
    alpha=0.2,
    label="True Elasticity Range",
)
plt.xlabel("Income")
plt.ylabel("Songs Sales Elasticity")
plt.title("Songs Sales Elasticity vs Income")
plt.legend(loc="lower right")

先来解读一下步骤，LinearDML初始化之后，就直接fit建模，
这里初始化了model_y + model_t两个模型，也就是在估计器里，当出现了X/T/W/Y拟合了两个模型：
$$\begin{gathered} log(Y)=\theta (X)\cdot log(T)+f(X,W)+ϵ \\ log(T)=g(X,W)+\eta \end{gathered}$$
此时可以说，X|T~Y的异质性处理效应就是弹性
（笔者喜欢在这里把X叫做异质变量，异质处理效应在这里就是需求弹性 ），
直接计算处理效应：est.effect(X_test)，这里可以看到估计效应只需要X(income)一列就可以了，
同时给到处理效应的区间估计effect_interval，
之后就是把，手算的真效应估计量，pred预测估计量进行对比：

从手算的弹性~income的关系来看，是一条非线性曲线，当收入小于1时，弹性在-1.75左右，当收入大于1时，有一个较小的负值。
模型拟合的是二次项的函数，拟合不是很好。
但它仍然抓住了总体趋势:弹性是负的，如果人们有更高的收入，他们对价格变化的敏感度就会降低。

# Get the final coefficient and intercept summary
est.summary()

输出结果：

“线性dml”估计器也可以返回最终模型的系数和截距的总结，包括点估计、p值和置信区间。
由上表可知，$income$具有正效应，${income}^2$具有负效应，两者在统计学上都是显著的。

这里的两个结果，coef 和 CATE分别代表了 X -> Y 和 T -> Y的影响
CATE的代表有优惠券则会降低收入水平。

2.5 X~Y 分析：非线性估计器：因果树CausalForestDML

# Train EconML model
est = CausalForestDML(
    model_y=GradientBoostingRegressor(), model_t=GradientBoostingRegressor()
)
est.fit(log_Y, log_T, X=X, W=W, inference="blb")
# Get treatment effect and its confidence interval
te_pred = est.effect(X_test)
te_pred_interval = est.effect_interval(X_test)

# Compare the estimate and the truth
plt.figure(figsize=(10, 6))
plt.plot(X_test.flatten(), te_pred, label="Sales Elasticity Prediction")
plt.plot(X_test.flatten(), truth_te_estimate, "--", label="True Elasticity")
plt.fill_between(
    X_test.flatten(),
    te_pred_interval[0],
    te_pred_interval[1],
    alpha=0.2,
    label="95% Confidence Interval",
)
plt.fill_between(
    X_test.flatten(),
    truth_te_lower,
    truth_te_upper,
    alpha=0.2,
    label="True Elasticity Range",
)
plt.xlabel("Income")
plt.ylabel("Songs Sales Elasticity")
plt.title("Songs Sales Elasticity vs Income")
plt.legend(loc="lower right")

我们注意到，这个模型比“线性dml”更适合，95%置信区间正确地覆盖了真实的处理效果估计，并捕捉了收入约为1时的变化。
总体而言，该模型表明，低收入人群比高收入人群对价格变化更加敏感。

2.6 X|T~Y分析：SingleTreeCateInterpreter哪些用户比较积极/消极

EconML包括可解释性工具，以更好地理解治疗效果。

治疗效果可能很复杂，但我们通常感兴趣的是一些简单的规则，这些规则可以区分哪些用户对提议的变化做出积极回应，哪些用户保持中立，哪些用户做出消极回应。

EconML SingleTreeCateInterpreter通过训练由任何EconML估计器输出的处理效果的单一决策树来提供可解释性。

intrp = SingleTreeCateInterpreter(include_model_uncerta
以上是关于因果推断笔记——因果图建模之微软开源的EconML的主要内容，如果未能解决你的问题，请参考以下文章
 
 因果推断笔记——因果图建模之微软开源的dowhy
 因果推断笔记——因果图建模之微软开源的dowhy
 因果推断笔记——工具变量内生性以及DeepIV
 因果推断笔记——工具变量内生性以及DeepIV
 因果推断笔记——因果图建模之Uber开源的CausalML
 因果推断笔记——因果图建模之Uber开源的CausalML