ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

因果推断实战指南:从RCT到工具变量与断点回归

因果推断实战指南:从RCT到工具变量与断点回归 简介因果推断是连接相关性与因果性的关键工具在数据科学、公共政策、医药与社会科学中应用广泛。这份资源以同名指南为核心面向希望系统掌握因果推断原理与实操方法的研究者、数据分析师及高年级学生内容覆盖随机对照试验、潜在结果框架、混杂变量控制以及倾向得分匹配、工具变量、差分法等主流技术同时讨论数据缺失、测量误差等现实挑战与人工智能时代的发展方向。资源共189个文件压缩包大小22.19MB以134个PNG示意图、23个IPython Notebook程序、21个CSV实验数据为主另含少量图片、Python脚本与说明文档。Notebook代码与配套数据可直接运行帮助读者边学边练理解如何区分相关性与因果效应、构建反事实模型并评估结果。目前已有62人学习浏览适合作为课程辅助、论文方法复现或项目落地前的速查手册。通过学习读者可获得从理论推导到代码实现的一套完整资源显著提升因果推断的应用信心。1. 因果推断不是统计学选修课而是数据产品的底线大多数数据分析师的第一课是相关矩阵第一工具是线性回归但这两样东西回答不了最值钱的问题如果我把价格提高10%销量会掉多少如果用户看了新的推荐算法停留时长是否真的变长相关性只能说“有关”因果推断要回答“如果”。哪怕是A/B测试也常因样本偏倚、违反稳定单元处理价值假设而失效。手头这份CSV数据集组合恰好覆盖了从随机对照试验到观察性研究的完整链路能让你用一两天时间把倾向得分匹配、工具变量、差分法、断点回归逐一跑通。它适合已经被回归搞到麻木、想认真区分相关与因果的工程师和数据科学家。2. 从潜在结果模型出发先搞清楚你在估什么效应2.1 潜在结果与ATE/ATT/CATE的定义因果推断的现代基础是Rubin的潜在结果框架。对每一个个体假设存在两个潜在结果处理状态下的结果Y(1)和控制状态下的结果Y(0)。个体因果效应定义为其差值Y(1) - Y(0)。然而一个个体在同一时间只能观察到其中一种结果这是因果推断的根本难题。研究者转而估计平均处理效应ATE即E[Y(1) - Y(0)]或者处理组平均处理效应ATT即E[Y(1) - Y(0) | D1]。随机化实验之所以被奉为金标准是因为它同时满足两个关键假定的强版本可忽略性处理分配独立于潜在结果也称无混杂和重叠假设每个个体都有一定概率进入处理组。当你拿到一份RCT数据时直接对比两组的均值差异就是ATE的无偏估计。但现实中数据分析师手里的数据大多是观察性的处理组和对照组在协变量分布上本身就不同这时候需要更复杂的识别策略。2.2 用medicine_impact_recovery.csv做RCT分析medicine_impact_recovery.csv模拟的是某个药物对患者康复周期的影响。我先假设它的列结构为patient_id、treatment0/1、age、severity病情严重度0-100、recovery_days。当数据来自完全随机化实验时可以直接用t检验或者线性回归。import pandas as pd import statsmodels.api as sm df pd.read_csv(medicine_impact_recovery.csv) # 先做协变量平衡性检查确认随机化是否成功 for col in [age, severity]: treated df.loc[df[treatment] 1, col] control df.loc[df[treatment] 0, col] print(col, diff:, treated.mean() - control.mean()) # 直接回归Y recovery_days, D treatment X df[[treatment, age, severity]] X sm.add_constant(X) model sm.OLS(df[recovery_days], X).fit() print(model.summary())treatment的系数就是调整了年龄和严重程度后的处理效应。如果随机化成功加入协变量前后系数变化很小如果变化很大说明随机化可能被破坏或者样本量太小导致分组不平衡。这里有两个参数需要注意add_constant会添加截距项否则回归会强制通过原点severity作为连续变量进入模型默认假设对恢复天数是线性影响这在医学数据里往往过于乐观常见做法是加入二次项或样条项。2.3 用learning_mindset.csv理解相关性与因果的差距learning_mindset.csv对应的场景是教育干预一组学生接受了“成长心态”辅导对照组没有观测变量包括pretest基线成绩、school_id、intervention、posttest。这类数据经常不是纯随机的因为学校可能自己选择是否参与干预。如果直接把干预组和对照组对比会混入学校层面的选择性偏差。# 查看不同学校的干预比例识别潜在聚类效应 import pandas as pd df pd.read_csv(learning_mindset.csv) print(df.groupby(school_id)[intervention].mean()) # 简单对比未调整的均值差 print(df.groupby(intervention)[posttest].mean())当你看到不同学校的干预比例不均匀时就要警惕了。一个常见做法是在回归中加入学校固定效应用C(school_id)来吸收学校层面的不变特征这就相当于在校内做比较。这个时候实际估计的效应是校内ATE受到该校样本量权重的影响。如果想要总体效应应该用学校的规模做加权。对于这类场景我会先用线性概率模型检查干预对posttest的影响然后对比不加固定效应和加固定效应的结果差异越大说明选择偏差越严重。场景识别策略关键假设主要风险完全随机RCT均值差/回归随机分配随机化失败干预非随机学校固定效应组内可比忽略动态混杂观察性数据PSM/DID无未观测混杂维度灾难3. 观察性数据里的主力工具倾向得分匹配与差分法3.1 为什么不能直接回归混杂与选择偏差当观测数据不能随机分配时直接回归暗含了一个强假设处理组和对照组在协变量完全可比的条件下才成立。比如customer_transactions.csv里高消费用户可能更倾向参与会员活动他们本身的消费趋势就跟普通用户不同。即使回归里控制了历史消费金额也无法排除某个未被观测到的变量同时在影响会员活动和当前消费。这时需要更结构化的识别策略。一个自然的思路是“造一个随机化”给定协变量X处理分配D独立于潜在结果。倾向得分e(X) P(D1|X)就是这个随机化概率的估计值。用逻辑回归估计倾向得分后可以做匹配、加权或分层。匹配的直观方式是把倾向得分相近的处理组和控制组个体配成对从而在协变量分布上实现平衡。3.2 PSM实战customer_transactions.csv配置customer_transactions.csv我处理过类似的结构customer_id、member_days入会时长、total_spend历史消费、promo_flag是否收到促销、future_spend未来消费。这里我们关心的是促销活动是否提升了未来消费但促销往往发给活跃用户活跃本身就会带来更多消费。import pandas as pd from sklearn.linear_model import LogisticRegression import numpy as np df pd.read_csv(customer_transactions.csv) # 特征工程对数化层数多的变量 df[log_spend] np.log1p(df[total_spend]) X df[[member_days, log_spend, prev_orders]] y df[promo_flag] # 训练倾向得分模型 lr LogisticRegression(max_iter1000) lr.fit(X, y) df[ps] lr.predict_proba(X)[:, 1] # 最近邻匹配1:1无放回 treated df[df[promo_flag] 1].sort_values(ps) control_pool df[df[promo_flag] 0] used set() matched_pairs [] for idx, row in treated.iterrows(): candidates control_pool.index.difference(used) if len(candidates) 0: break c control_pool.loc[candidates] diff (c[ps] - row[ps]).abs().argmin() matched_control_idx c.index[diff] used.add(matched_control_idx) matched_pairs.append((idx, matched_control_idx))这段代码的关键在于argmin会找距离最近的对照组个体但需要注意卡尺问题如果最小距离仍然超过0.01应该将该处理组样本剔除。我没有在代码里加卡尺判断实际中建议改用sklearn.neighbors.NearestNeighbors并设置radius或者用pymatch库。匹配之后必须做平衡性验证检查匹配后各协变量的标准化差异是否小于0.1否则匹配不充分。3.3 DID实战ice_cream_sales.csv与政策冲击差分法DID适用于面板数据或重复截面数据某个政策、事件冲击在特定时点发生我们比较处理组在冲击前后的变化量与对照组同期变化量的差值。ice_cream_sales.csv记录了两家店在夏季某周突然调价后的销量一家调价处理组一家未调价对照组。周度数据天然适合DID。import pandas as pd import statsmodels.api as sm df pd.read_csv(ice_cream_sales.csv) # 假设列: store_id, week, price_change(1/0), post(1/0), sales df[did] df[price_change] * df[post] X df[[price_change, post, did]] X sm.add_constant(X) model sm.OLS(df[sales], X).fit() print(model.summary()) # did系数即价格调降对销量的因果效应DID的核心假设是平行趋势如果调价没有发生处理组和对照组的销量变化趋势一致。这个假设无法直接检验但可以通过画图验证冲击前几个周期的趋势是否平行。还有一个细节是标准误需要聚类到store_id层面因为同一个店的跨期误差存在自相关。上面代码没有聚类实践中应使用model.get_robustcov_results(cov_typecluster, groupsdf[store_id])。4. 工具变量与断点回归随机化失效时的备选方案4.1 工具变量的适用条件与诊断当存在未观测混杂且没有自然实验时工具变量IV是最后武器之一。有效的工具变量Z要满足两条相关性Z影响处理D和排他性Z只通过D影响结果Y。Z的常见来源包括制度规则、距离远近、时间窗口等。ak91.csv是经典的Angrist-Krueger教育回报数据用出生季度作为受教育年限的工具变量。理由是晚出生的孩子入学年龄偏晚在义务教育法允许退学的年龄时已经接受更多教育而出生季度本身不应该直接影响工资除非通过教育。用工具变量做估计常用两阶段最小二乘2SLS。第一阶段用Z预测D第二阶段用D的预测值预测Y。在Python里可以直接使用linearmodels.IV2SLS。4.2 用ak91.csv做2SLS估计import pandas as pd from linearmodels.iv import IV2SLS df pd.read_csv(ak91.csv) # 假设列: wage, education, quarter_of_birth, age, married df[educ_int] df[education].astype(int) exog df[[age, married]] exog sm.add_constant(exog) endog df[wage] # 或log_wage treatment df[educ_int] instruments pd.get_dummies(df[quarter_of_birth], prefixq, drop_firstTrue) iv_model IV2SLS(endog, exog, treatment, instruments).fit() print(iv_model.summary)IV2SLS的第一个参数是结果变量第二个是外生控制变量第三个是内生处理变量第四个是工具变量。检查第一阶段F统计量一般要大于10否则存在弱工具变量问题。这里的quarter_of_birth被拆成了三个虚拟变量相当于用出生季度的非线性形式做工具。排他性检验只能靠逻辑论证无法用统计方法直接验证。如果结果不稳健可以尝试只用出生一季度作为单一工具看系数变化范围。4.3 断点回归设计在enem_scores.csv上的应用断点回归RDD利用某个强制规则造成的阈值比如成绩超过某个分数就能获得奖学金比较阈值附近学生的表现。enem_scores.csv存储了巴西高考ENEM成绩、是否获得助学金、后续大学表现。如果助学金由分数是否高于600分决定就能在600分附近构造一个准实验。import pandas as pd import numpy as np df pd.read_csv(enem_scores.csv) # 假设列: score, scholarship(600), college_grade df[above] (df[score] 600).astype(int) # 取阈值附近±10分的带宽 bw 10 rdd_df df[(df[score] 600 - bw) (df[score] 600 bw)] # 局部线性回归左侧与右侧分别拟合 for side in [0, 1]: sub rdd_df[rdd_df[above] side] X sub[score] - 600 X sm.add_constant(X) model sm.OLS(sub[college_grade], X).fit() print(fSide {side}: intercept {model.params[0]:.3f})在600分左侧和右侧分别拟合直线截距之差就是局部平均处理效应。带宽选择是关键带宽太窄样本少、方差大带宽太宽则引入非线性偏差。主流做法是使用MSE最优带宽如rdrobust包计算。RDD只对阈值附近个体有效外部效度有限但内部效度非常高。方法数据类型需要假设输出效应2SLS截面排他性、相关性LATERDD截面/面板连续性、无操纵阈值附近LATE5. 数据中的系统偏差缺失值、样本选择与稳健性检查5.1 从invest_email_biased.csv看选择偏差invest_email_biased.csv对应的是营销邮件的A/B测试数据但回拒信的用户被排除在分析之外导致样本不再是随机样本。假设实验发出投资建议邮件处理组和对照组各5000人但只有50%的用户打开了邮件并记录了结果。如果我们只分析这些“打开并响应”的用户就会犯内生分层错误打开行为跟用户本身的投资倾向相关这个倾向也让用户更容易获得高回报。这种偏差无法通过重新加权修复除非知道未打开用户的反事实结果。5.2 缺失数据的处理与敏感性分析collections_email.csv可能记录了催收邮件的完整发送状态和后续还款行为。处理缺失数据的常见套路是三种完整记录分析列表删除、多重插补、逆概率加权。但因果推断场景对插补要格外谨慎如果缺失机制与处理或结果相关插补会引入偏倚更安全的做法是用缺失指示器加多重插补做敏感性分析。import pandas as pd from sklearn.impute import IterativeImputer df pd.read_csv(collections_email.csv) # 假设列: email_sent(1/0), opened, due_amount, payment_after missing_cols [opened, due_amount] for col in missing_cols: df[col _missing] df[col].isna().astype(int) # 用多重插补填充并对比两种方式下的效应估计 imp IterativeImputer(max_iter10, random_state42) df[[opened, due_amount]] imp.fit_transform(df[[opened, due_amount]])参数max_iter10表示多重插补循环次数random_state固定随机种子来保证结果可复现。多重插补产生多个数据集最后用Rubin法则合并估计值。但在因果推断中我更倾向把它当作稳健性检验而非主结果。如果完整记录分析和插补后的结论一致说明结论对缺失机制不敏感如果不一致必须明确报告缺失机制造成了多大影响。5.3 最后的检查清单当你用这些数据集跑完一个因果效应估计至少要回答以下几个问题识别假设是什么随机化、可忽略性、平行趋势还是排他性处理组和对照组是否在协变量上平衡标准化差异小于0.1吗有没有做placebo测试用虚假处理或虚假时间点重跑模型系数应当不显著。标准误是否聚类聚类层级是否与处理分配或抽样层级一致最后的技巧把多份CSV文件看作同一套因果推断流程的不同阶段而不是孤立的数据集。先用medicine_impact_recovery.csv验证你对RCT的理解再用ice_cream_sales.csv练习DID最后用ak91.csv挑战工具变量。每一个文件都配一个“如果假设不成立会看到什么迹象”的检查项这比背一堆统计测试更有用。本文还有配套的精品资源点击获取
返回列表