【因果推断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: 多模块(.py)中全局变量的导入

文章目录 global关键字可变类型和不可变类型数据的内存地址单模块(单个py文件)的全局变量示例总结 多模块(多个py文件)的全局变量from x import x导入全局变量示例 import x导入全局变量示例 总结 global关键字 global 的作用范围是模块(.py)级别: 当你在一个模块(文件)中使用 global 声明变量时,这个变量只在该模块的全局命名空

【Python编程】Linux创建虚拟环境并配置与notebook相连接

1.创建 使用 venv 创建虚拟环境。例如,在当前目录下创建一个名为 myenv 的虚拟环境: python3 -m venv myenv 2.激活 激活虚拟环境使其成为当前终端会话的活动环境。运行: source myenv/bin/activate 3.与notebook连接 在虚拟环境中,使用 pip 安装 Jupyter 和 ipykernel: pip instal

【机器学习】高斯过程的基本概念和应用领域以及在python中的实例

引言 高斯过程(Gaussian Process,简称GP)是一种概率模型,用于描述一组随机变量的联合概率分布,其中任何一个有限维度的子集都具有高斯分布 文章目录 引言一、高斯过程1.1 基本定义1.1.1 随机过程1.1.2 高斯分布 1.2 高斯过程的特性1.2.1 联合高斯性1.2.2 均值函数1.2.3 协方差函数(或核函数) 1.3 核函数1.4 高斯过程回归(Gauss

【学习笔记】 陈强-机器学习-Python-Ch15 人工神经网络(1)sklearn

系列文章目录 监督学习:参数方法 【学习笔记】 陈强-机器学习-Python-Ch4 线性回归 【学习笔记】 陈强-机器学习-Python-Ch5 逻辑回归 【课后题练习】 陈强-机器学习-Python-Ch5 逻辑回归(SAheart.csv) 【学习笔记】 陈强-机器学习-Python-Ch6 多项逻辑回归 【学习笔记 及 课后题练习】 陈强-机器学习-Python-Ch7 判别分析 【学

nudepy,一个有趣的 Python 库!

更多资料获取 📚 个人网站:ipengtao.com 大家好,今天为大家分享一个有趣的 Python 库 - nudepy。 Github地址:https://github.com/hhatto/nude.py 在图像处理和计算机视觉应用中,检测图像中的不适当内容(例如裸露图像)是一个重要的任务。nudepy 是一个基于 Python 的库,专门用于检测图像中的不适当内容。该

国产游戏崛起:技术革新与文化自信的双重推动

近年来,国产游戏行业发展迅猛,技术水平和作品质量均得到了显著提升。特别是以《黑神话:悟空》为代表的一系列优秀作品,成功打破了过去中国游戏市场以手游和网游为主的局限,向全球玩家展示了中国在单机游戏领域的实力与潜力。随着中国开发者在画面渲染、物理引擎、AI 技术和服务器架构等方面取得了显著进展,国产游戏正逐步赢得国际市场的认可。然而,面对全球游戏行业的激烈竞争,国产游戏技术依然面临诸多挑战,未来的

pip-tools:打造可重复、可控的 Python 开发环境,解决依赖关系,让代码更稳定

在 Python 开发中,管理依赖关系是一项繁琐且容易出错的任务。手动更新依赖版本、处理冲突、确保一致性等等,都可能让开发者感到头疼。而 pip-tools 为开发者提供了一套稳定可靠的解决方案。 什么是 pip-tools? pip-tools 是一组命令行工具,旨在简化 Python 依赖关系的管理,确保项目环境的稳定性和可重复性。它主要包含两个核心工具:pip-compile 和 pip

HTML提交表单给python

python 代码 from flask import Flask, request, render_template, redirect, url_forapp = Flask(__name__)@app.route('/')def form():# 渲染表单页面return render_template('./index.html')@app.route('/submit_form',

Python QT实现A-star寻路算法

目录 1、界面使用方法 2、注意事项 3、补充说明 用Qt5搭建一个图形化测试寻路算法的测试环境。 1、界面使用方法 设定起点: 鼠标左键双击,设定红色的起点。左键双击设定起点,用红色标记。 设定终点: 鼠标右键双击,设定蓝色的终点。右键双击设定终点,用蓝色标记。 设置障碍点: 鼠标左键或者右键按着不放,拖动可以设置黑色的障碍点。按住左键或右键并拖动,设置一系列黑色障碍点

Python:豆瓣电影商业数据分析-爬取全数据【附带爬虫豆瓣,数据处理过程,数据分析,可视化,以及完整PPT报告】

**爬取豆瓣电影信息,分析近年电影行业的发展情况** 本文是完整的数据分析展现,代码有完整版,包含豆瓣电影爬取的具体方式【附带爬虫豆瓣,数据处理过程,数据分析,可视化,以及完整PPT报告】   最近MBA在学习《商业数据分析》,大实训作业给了数据要进行数据分析,所以先拿豆瓣电影练练手,网络上爬取豆瓣电影TOP250较多,但对于豆瓣电影全数据的爬取教程很少,所以我自己做一版。 目