【因果推断python】26_双重稳健估计1

2024-06-11 06:04

本文主要是介绍【因果推断python】26_双重稳健估计1,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

目录

不要把所有的鸡蛋放在一个篮子里

双重稳健估计

关键思想


不要把所有的鸡蛋放在一个篮子里

我们已经学会了如何使用线性回归和倾向得分加权来估计 E[Y|T=1]-E[Y|T=0]|X。但是我们应该在什么时候使用哪一个呢?在不明确的情况下,请同时使用两者!双重稳健估计是一种将倾向得分和线性回归相结合的方法,您不必依赖它们中的任何一种。

为了了解这是如何工作的,让我们考虑一下心态实验。这是一项在美国公立高中进行的随机研究,旨在发现成长心态的影响。它的工作方式是学校邀请学生参加一个研讨会,向他们灌输一种成长的心态。然后,他们跟踪学生在大学期间的表现,并衡量他们在学业上的表现。这个衡量结果被编译为标准化的成就分数。为了保护学生的隐私,这项研究的真实数据没有公开。但是,我们有一个与 Athey 和 Wager 提供的统计属性相同的模拟数据集,因此我们将改为使用这个数据来进行分析。

import warnings
warnings.filterwarnings('ignore')import pandas as pd
import numpy as np
from matplotlib import style
from matplotlib import pyplot as plt
import seaborn as sns%matplotlib inlinestyle.use("fivethirtyeight")
pd.set_option("display.max_columns", 6)
data = pd.read_csv("./data/learning_mindset.csv")
data.sample(5, random_state=5)

虽然这项研究做了随机化处理,但这些数据并不是没有出现混淆的情况。其中一个可能的原因是,干预变量是通过学生是否收到研讨会邀请来衡量的。因此,尽管被邀请参与的机会是随机的,但是否真的参与却不是。我们在这里处理一个不服从(non-compliance)的情况。这方面的一个证据是学生对成功的期望是如何与是否参加研讨会相关联的。自我报告中期望较高的学生更有可能参加成长心态研讨会。

data.groupby("success_expect")["intervention"].mean()

正如已经学习到的,我们可以通过使用线性回归或者逻辑回归估计倾向得分模型的方法来调整不服从的情况。在做回归之前,我们需要将份类变量转化为虚拟变量。

categ = ["ethnicity", "gender", "school_urbanicity"]
cont = ["school_mindset", "school_achievement", "school_ethnic_minority", "school_poverty", "school_size"]data_with_categ = pd.concat([data.drop(columns=categ), # dataset without the categorical featurespd.get_dummies(data[categ], columns=categ, drop_first=False) # categorical features converted to dummies
], axis=1)print(data_with_categ.shape)
(10391, 32)

我们现在已经准备好了解双重稳健估计的工作原理。

双重稳健估计\hat{\mu_0}(x)

我不会推导出估算器,而是首先向您展示它,然后才告诉您为什么它很棒。

A\hat{T}E=\frac1N\sum\left(\frac{T_i(Y_i-\hat{\mu_1}(X_i))}{\hat{P}(X_i)}+\hat{\mu_1}(X_i)\right)-\frac1N\sum\left(\frac{(1-T_i)(Y_i-\hat{\mu_0}(X_i))}{1-\hat{P}(X_i)}+\hat{\mu_0}(X_i)\right)

其中 \hat{P}(x) 是对倾向得分的估计(例如,使用逻辑回归),\hat{\mu_1}(x) 是对倾向得分的估计E[Y|X,T=1](例如使用线性回归),而 \mu_0(x) 是对 E[Y|X,T=0]。正如您可能已经猜到的那样,双重稳健估计器的第一部分估计 E[Y_{1}],第二部分估计 E[Y_{0}]。让我们检查第一部分,因为所有直觉也将通过类比适用于第二部分。

因为我一开始就知道这个公式很吓人(但别担心,你会看到它超级简单),我将首先展示如何编写这个估计器。我觉得有些人对代码的恐惧不如对公式的恐惧。让我们看看这个估计器在实践中是如何工作的,好吗?

from sklearn.linear_model import LogisticRegression, LinearRegressiondef doubly_robust(df, X, T, Y):ps = LogisticRegression(C=1e6).fit(df[X], df[T]).predict_proba(df[X])[:, 1]mu0 = LinearRegression().fit(df.query(f"{T}==0")[X], df.query(f"{T}==0")[Y]).predict(df[X])mu1 = LinearRegression().fit(df.query(f"{T}==1")[X], df.query(f"{T}==1")[Y]).predict(df[X])return (np.mean(df[T]*(df[Y] - mu1)/ps + mu1) -np.mean((1-df[T])*(df[Y] - mu0)/(1-ps) + mu0))
T = 'intervention'
Y = 'achievement_score'
X = data_with_categ.columns.drop(['schoolid', T, Y])doubly_robust(data_with_categ, X, T, Y)
0.3882222817222756

双重稳健估计者说,就成就而言,我们应该期望参加心态研讨会的个人比未经治疗的同伴高 0.388 个标准差。 再一次,我们可以使用 bootstrap 来构建置信区间。

from joblib import Parallel, delayed # for parallel processingnp.random.seed(88)
# run 1000 bootstrap samples
bootstrap_sample = 1000
ates = Parallel(n_jobs=4)(delayed(doubly_robust)(data_with_categ.sample(frac=1, replace=True), X, T, Y)for _ in range(bootstrap_sample))
ates = np.array(ates)
print(f"ATE 95% CI:", (np.percentile(ates, 2.5), np.percentile(ates, 97.5)))
sns.distplot(ates, kde=False)
plt.vlines(np.percentile(ates, 2.5), 0, 20, linestyles="dotted")
plt.vlines(np.percentile(ates, 97.5), 0, 20, linestyles="dotted", label="95% CI")
plt.title("ATE Bootstrap Distribution")
plt.legend();

现在我们已经了解了双重稳健估计器,让我们来看看为什么它如此出色。首先,它被称为双重鲁棒,因为它只需要模型之一,\hat{P}(x)\hat{\mu}(x)是正确的指定的。要了解这一点,请查看估计 E[Y_{1}]的第一部分并仔细查看它。

\hat{E}[Y_1]=\frac{1}{N}\sum\left(\frac{T_i(Y_i-\hat{\mu_1}(X_i))}{\hat{P}(X_i)}+\hat{\mu_1}(X_i)\right)

假设 \hat{\mu_1}(x) 是正确的。如果倾向得分模型是错误的,我们不需要担心。因为如果\hat{\mu_1}(x) 是正确的,那么 E[T_i(Y_i-\hat{\mu_1}(X_i))]=0。那是因为 T_{i}的乘法只选择了被处理的,并且\hat{\mu_1} 在被处理的残差上,根据定义,均值为零。这导致整个公式等于 \hat{\mu_1}(X_i),这是通过假设正确估计的 E[Y_{1}]。所以,你看,通过正确,\hat{\mu_1}(X_i) 消除了倾向得分模型的相关性。我们可以应用相同的推理来理解 E[Y_0] 的估计量。

但不要相信我的话。让代码告诉你方向!在下面的估计器中,我用一个从 0.1 到 0.9 的随机统一变量替换了估计倾向得分的逻辑回归(我不希望非常小的权重破坏我的倾向得分方差)。由于这是随机的,所以它不可能是一个好的倾向得分模型,但我们将看到双重稳健估计器仍然设法产生一个非常接近于使用逻辑回归估计倾向得分时的估计值。

from sklearn.linear_model import LogisticRegression, LinearRegressiondef doubly_robust_wrong_ps(df, X, T, Y): # wrong PS model np.random.seed(654) ps = np.random.uniform(0.1, 0.9, df.shape[0]) mu0 = LinearRegression().fit(df.query(f"{T}==0")[X], df.query(f"{T}==0")[Y]).predict(df[X]) mu1 = LinearRegression().fit(df.query(f"{T}==1")[X], df.query(f"{T}==1")[Y]).predict(df[X]) return ( np.mean(df[T](df[Y] - mu1)/ps + mu1) - np.mean((1-df[T])(df[Y] - mu0)/(1-ps) + mu0) )doubly_robust_wrong_ps(data_with_categ, X, T, Y)

如果我们使用自助采样法,我们可以看到,相比基于逻辑回归的倾向得分,方差会稍高一点。

np.random.seed(88)
parallel_fn = delayed(doubly_robust_wrong_ps)
wrong_ps = Parallel(n_jobs=4)(parallel_fn(data_with_categ.sample(frac=1, replace=True), X, T, Y)for _ in range(bootstrap_sample))
wrong_ps = np.array(wrong_ps)
print(f"ATE 95% CI:", (np.percentile(ates, 2.5), np.percentile(ates, 97.5)))

这涵盖了倾向模型错误但结果模型正确的情况。其他情况呢?让我们再好好看看估计器的第一部分,但让我们重新排列一些术语

\hat{E}[Y_1]=\frac1N\sum\left(\frac{T_i(Y_i-\hat{\mu_1}(X_i))}{\hat{P}(X_i)}+\hat{\mu_1}(X_i)\right)

\hat{E}[Y_1]=\frac1N\sum\left(\frac{T_iY_i}{\hat{P}(X_i)}-\frac{T_i\hat{\mu_1}(X_i)}{\hat{P}(X_i)}+\hat{\mu_1}(X_i)\right)

\hat{E}[Y_1]=\frac1N\sum\left(\frac{T_iY_i}{\hat{P}(X_i)}-\left(\frac{T_i}{\hat{P}(X_i)}-1\right)\hat{\mu_1}(X_i)\right)

\hat{E}[Y_1]=\frac1N\sum\left(\frac{T_iY_i}{\hat{P}(X_i)}-\left(\frac{T_i-\hat{P}(X_i)}{\hat{P}(X_i)}\right)\hat{\mu_1}(X_i)\right)

现在,假设正确指定了倾向得分 \hat{P}(X_i)。在这种情况下,E[T_i-\hat{P}(X_i)]=0,它消除了依赖于 \hat{\mu_1}(X_i) 的部分。这使得双重鲁棒估计器减少为倾向得分加权估计器 \frac{T_iY_i}{\hat{P}(X_i)},假设是正确的。因此,即使 \hat{\mu_1}(X_i)是错误的,只要正确指定了倾向得分,估计器仍然是正确的。

再一次,如果你更相信代码而不是公式,这里就是实际验证。在下面的代码中,我用随机正态变量替换了两个回归模型。毫无疑问 \hat{\mu}(X_i)是错误的。尽管如此,我们仍将看到双重稳健估计仍设法恢复我们之前看到的大约 0.38 的相同 \hat{ATE}

from sklearn.linear_model import LogisticRegression, LinearRegressiondef doubly_robust_wrong_model(df, X, T, Y):np.random.seed(654)ps = LogisticRegression(C=1e6).fit(df[X], df[T]).predict_proba(df[X])[:, 1]# wrong mu(x) modelmu0 = np.random.normal(0, 1, df.shape[0])mu1 = np.random.normal(0, 1, df.shape[0])return (np.mean(df[T]*(df[Y] - mu1)/ps + mu1) -np.mean((1-df[T])*(df[Y] - mu0)/(1-ps) + mu0))
doubly_robust_wrong_model(data_with_categ, X, T, Y)

同样的,我们可以通过使用自助采样法看到方差还是相对高一点。

np.random.seed(88)
parallel_fn = delayed(doubly_robust_wrong_model)
wrong_mux = Parallel(n_jobs=4)(parallel_fn(data_with_categ.sample(frac=1, replace=True), X, T, Y)for _ in range(bootstrap_sample))
wrong_mux = np.array(wrong_mux)
print(f"ATE 95% CI:", (np.percentile(ates, 2.5), np.percentile(ates, 97.5)))

我希望我已经让你相信双重稳健估计的力量。它之所以神奇,是因为在因果推理中,有两种方法可以从我们的因果估计中消除偏见:您可以对干预机制或结果机制进行建模如果这些模型中的任何一个都是正确的,那么您就可以开始了。

一个需要警惕的地方是,在实践中,很难对其中任何一个进行精确建模更常见的情况是,倾向得分和结果模型都不是 100% 正确的。他们都错了,但方式不同。发生这种情况时, 是使用单一模型最好还是使用双重稳健估计更佳,目前还没有一个定论[1] [2] [3]。至于我,我仍然喜欢对两种方法都考虑一下,因为至少给了我两种正确的可能性。

关键思想

在这里,我们看到了一种将线性回归与倾向得分相结合的简单方法,以产生双重稳健的估计量。这个估计器之所以有这个名字,是因为它只需要一个模型是正确的。如果倾向得分模型是正确的,即使结果模型是错误的,我们也能够识别因果效应。另一方面,如果结果模型是正确的,即使倾向评分模型是错误的,我们也能够识别因果效应。

这篇关于【因果推断python】26_双重稳健估计1的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/1050363

相关文章

如何使用 Python 读取 Excel 数据

《如何使用Python读取Excel数据》:本文主要介绍使用Python读取Excel数据的详细教程,通过pandas和openpyxl,你可以轻松读取Excel文件,并进行各种数据处理操... 目录使用 python 读取 Excel 数据的详细教程1. 安装必要的依赖2. 读取 Excel 文件3. 读

Python的time模块一些常用功能(各种与时间相关的函数)

《Python的time模块一些常用功能(各种与时间相关的函数)》Python的time模块提供了各种与时间相关的函数,包括获取当前时间、处理时间间隔、执行时间测量等,:本文主要介绍Python的... 目录1. 获取当前时间2. 时间格式化3. 延时执行4. 时间戳运算5. 计算代码执行时间6. 转换为指

利用Python调试串口的示例代码

《利用Python调试串口的示例代码》在嵌入式开发、物联网设备调试过程中,串口通信是最基础的调试手段本文将带你用Python+ttkbootstrap打造一款高颜值、多功能的串口调试助手,需要的可以了... 目录概述:为什么需要专业的串口调试工具项目架构设计1.1 技术栈选型1.2 关键类说明1.3 线程模

Python ZIP文件操作技巧详解

《PythonZIP文件操作技巧详解》在数据处理和系统开发中,ZIP文件操作是开发者必须掌握的核心技能,Python标准库提供的zipfile模块以简洁的API和跨平台特性,成为处理ZIP文件的首选... 目录一、ZIP文件操作基础三板斧1.1 创建压缩包1.2 解压操作1.3 文件遍历与信息获取二、进阶技

Python Transformers库(NLP处理库)案例代码讲解

《PythonTransformers库(NLP处理库)案例代码讲解》本文介绍transformers库的全面讲解,包含基础知识、高级用法、案例代码及学习路径,内容经过组织,适合不同阶段的学习者,对... 目录一、基础知识1. Transformers 库简介2. 安装与环境配置3. 快速上手示例二、核心模

Python正则表达式语法及re模块中的常用函数详解

《Python正则表达式语法及re模块中的常用函数详解》这篇文章主要给大家介绍了关于Python正则表达式语法及re模块中常用函数的相关资料,正则表达式是一种强大的字符串处理工具,可以用于匹配、切分、... 目录概念、作用和步骤语法re模块中的常用函数总结 概念、作用和步骤概念: 本身也是一个字符串,其中

Python使用getopt处理命令行参数示例解析(最佳实践)

《Python使用getopt处理命令行参数示例解析(最佳实践)》getopt模块是Python标准库中一个简单但强大的命令行参数处理工具,它特别适合那些需要快速实现基本命令行参数解析的场景,或者需要... 目录为什么需要处理命令行参数?getopt模块基础实际应用示例与其他参数处理方式的比较常见问http

python实现svg图片转换为png和gif

《python实现svg图片转换为png和gif》这篇文章主要为大家详细介绍了python如何实现将svg图片格式转换为png和gif,文中的示例代码讲解详细,感兴趣的小伙伴可以跟随小编一起学习一下... 目录python实现svg图片转换为png和gifpython实现图片格式之间的相互转换延展:基于Py

Python中的getopt模块用法小结

《Python中的getopt模块用法小结》getopt.getopt()函数是Python中用于解析命令行参数的标准库函数,该函数可以从命令行中提取选项和参数,并对它们进行处理,本文详细介绍了Pyt... 目录getopt模块介绍getopt.getopt函数的介绍getopt模块的常用用法getopt模

Python利用ElementTree实现快速解析XML文件

《Python利用ElementTree实现快速解析XML文件》ElementTree是Python标准库的一部分,而且是Python标准库中用于解析和操作XML数据的模块,下面小编就来和大家详细讲讲... 目录一、XML文件解析到底有多重要二、ElementTree快速入门1. 加载XML的两种方式2.