【因果推断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 FastAPI+Celery+RabbitMQ实现分布式图片水印处理系统

《PythonFastAPI+Celery+RabbitMQ实现分布式图片水印处理系统》这篇文章主要为大家详细介绍了PythonFastAPI如何结合Celery以及RabbitMQ实现简单的分布式... 实现思路FastAPI 服务器Celery 任务队列RabbitMQ 作为消息代理定时任务处理完整

Python Websockets库的使用指南

《PythonWebsockets库的使用指南》pythonwebsockets库是一个用于创建WebSocket服务器和客户端的Python库,它提供了一种简单的方式来实现实时通信,支持异步和同步... 目录一、WebSocket 简介二、python 的 websockets 库安装三、完整代码示例1.

揭秘Python Socket网络编程的7种硬核用法

《揭秘PythonSocket网络编程的7种硬核用法》Socket不仅能做聊天室,还能干一大堆硬核操作,这篇文章就带大家看看Python网络编程的7种超实用玩法,感兴趣的小伙伴可以跟随小编一起... 目录1.端口扫描器:探测开放端口2.简易 HTTP 服务器:10 秒搭个网页3.局域网游戏:多人联机对战4.

使用Python实现快速搭建本地HTTP服务器

《使用Python实现快速搭建本地HTTP服务器》:本文主要介绍如何使用Python快速搭建本地HTTP服务器,轻松实现一键HTTP文件共享,同时结合二维码技术,让访问更简单,感兴趣的小伙伴可以了... 目录1. 概述2. 快速搭建 HTTP 文件共享服务2.1 核心思路2.2 代码实现2.3 代码解读3.

Python使用自带的base64库进行base64编码和解码

《Python使用自带的base64库进行base64编码和解码》在Python中,处理数据的编码和解码是数据传输和存储中非常普遍的需求,其中,Base64是一种常用的编码方案,本文我将详细介绍如何使... 目录引言使用python的base64库进行编码和解码编码函数解码函数Base64编码的应用场景注意

Python基于wxPython和FFmpeg开发一个视频标签工具

《Python基于wxPython和FFmpeg开发一个视频标签工具》在当今数字媒体时代,视频内容的管理和标记变得越来越重要,无论是研究人员需要对实验视频进行时间点标记,还是个人用户希望对家庭视频进行... 目录引言1. 应用概述2. 技术栈分析2.1 核心库和模块2.2 wxpython作为GUI选择的优

Python如何使用__slots__实现节省内存和性能优化

《Python如何使用__slots__实现节省内存和性能优化》你有想过,一个小小的__slots__能让你的Python类内存消耗直接减半吗,没错,今天咱们要聊的就是这个让人眼前一亮的技巧,感兴趣的... 目录背景:内存吃得满满的类__slots__:你的内存管理小助手举个大概的例子:看看效果如何?1.

Python+PyQt5实现多屏幕协同播放功能

《Python+PyQt5实现多屏幕协同播放功能》在现代会议展示、数字广告、展览展示等场景中,多屏幕协同播放已成为刚需,下面我们就来看看如何利用Python和PyQt5开发一套功能强大的跨屏播控系统吧... 目录一、项目概述:突破传统播放限制二、核心技术解析2.1 多屏管理机制2.2 播放引擎设计2.3 专

Python中随机休眠技术原理与应用详解

《Python中随机休眠技术原理与应用详解》在编程中,让程序暂停执行特定时间是常见需求,当需要引入不确定性时,随机休眠就成为关键技巧,下面我们就来看看Python中随机休眠技术的具体实现与应用吧... 目录引言一、实现原理与基础方法1.1 核心函数解析1.2 基础实现模板1.3 整数版实现二、典型应用场景2

Python实现无痛修改第三方库源码的方法详解

《Python实现无痛修改第三方库源码的方法详解》很多时候,我们下载的第三方库是不会有需求不满足的情况,但也有极少的情况,第三方库没有兼顾到需求,本文将介绍几个修改源码的操作,大家可以根据需求进行选择... 目录需求不符合模拟示例 1. 修改源文件2. 继承修改3. 猴子补丁4. 追踪局部变量需求不符合很