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

Posted 2021-09-10 悟乙己

文章目录

• 1 EconML介绍
• • 1.1 EconML介绍
  • 1.2 一些理论解答
• 2 智能营销案例：不同优惠折扣下的用户反应
• • 2.1 背景
  • 2.2 数据准备与理解
  • 2.3 生成假的Groud Truth的估计效应
  • 2.4 线性估计器LinearDML
• 相关问题
• • 案例：EconML + SHAP丰富可解释性

---

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) ，也就是说，每个子组的平均处理效应被称为条件平均处理效应，以子组内的分类方式为条件。

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

参考：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 线性估计器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")

相关问题

案例：EconML + SHAP丰富可解释性

Interpretability with SHAP

类似于如何用SHAP解释黑箱预测机器学习模型，我们也可以解释黑箱效应异质性模型。这种方法解释了为什么异质因果效应模型对特定人群产生了较大或较小的效应值。是哪些特征导致了这种差异?当模型被简洁地描述时，这个问题很容易解决，例如线性异质性模型的情况，人们可以简单地研究模型的系数。然而，当人们开始使用更具表现力的模型(如随机森林和因果森林)来建模效果异质性时，就会变得困难。SHAP值对于理解模型从训练数据中提取的影响异质性的主导因素有很大的帮助。
我们的软件包提供了与SHAP库的无缝集成。每个CateEstimator都有一个方法shap_values，它返回每个处理和结果对的估计器输出的SHAP值解释。然后，可以使用SHAP库提供的大量可视化功能对这些值进行可视化。此外，只要有可能，我们的库就会为每种最终模型类型从SHAP库中调用快速专用算法，这可以大大减少计算时间。