因果推断反事实归因：从PN/PS公式到Python代码的3个实战案例

因果推断反事实归因从PN/PS公式到Python代码的3个实战案例在数据科学领域我们常常需要回答如果...会怎样这类问题。比如如果用户没有点击这个广告他们的购买概率会降低多少如果患者服用了另一种药物康复时间会缩短多少这类问题本质上都是在探讨因果关系而反事实推理正是解决这类问题的利器。反事实推理的核心思想是构建一个与现实相反的假设场景通过对比实际结果与假设结果来推断因果关系。这种思维方式不仅存在于学术研究中也渗透在我们的日常决策里——当你想着如果昨天早点出门就不会迟到时就已经在无意识地进行反事实推理了。1. 反事实推理基础与核心公式1.1 必要性概率(PN)与充分性概率(PS)在反事实推理中两个核心概念是必要性概率(PN)和充分性概率(PS)。让我们通过一个简单的例子来理解它们假设X表示戴口罩(X1)或不戴口罩(X0)Y表示感染(Y0)或未感染(Y1)。对于一位戴口罩且未感染的人(A先生)PN回答的问题是如果A先生没有戴口罩他被感染的概率是多少数学表达式为PN P(Y_{X0}0 | X1, Y1)而对于一位不戴口罩且被感染的人(B先生)PS回答的问题是如果B先生戴了口罩他不会被感染的概率是多少数学表达式为PS P(Y_{X1}1 | X0, Y0)1.2 单调性假设下的PN计算在实际应用中PN的计算常引入单调性假设即干预的结果不会比不干预更差(Y_X1 ≥ Y_X0)。在单调性假设下PN可以表示为def calculate_pn(p_y1, p_y1_given_do_x0, p_x1y1): 计算单调性假设下的必要性概率(PN) 参数: p_y1: P(Y1) p_y1_given_do_x0: P(Y1|do(X0)) p_x1y1: P(X1,Y1) 返回: PN值 err (p_y1 - p_y1_given_do_x0) / p_x1y1 return err这个公式包含两个关键部分过剩风险率(ERR)反映干预与不干预的风险差异混杂因子(CF)校正观测数据与实验数据的偏差1.3 非单调情况下的边界估计当单调性假设不成立时我们可以估计PN的上下界def pn_bounds(p_y1, p_y1_given_do_x0, p_x0y0, p_x1y1, p_y0_given_x1): 计算PN的上下界 参数: p_y1: P(Y1) p_y1_given_do_x0: P(Y1|do(X0)) p_x0y0: P(X0,Y0) p_x1y1: P(X1,Y1) p_y0_given_x1: P(Y0|X1) 返回: (下界, 上界) # 计算ERR和CF err (p_y1 - p_y1_given_do_x0) / p_x1y1 cf (p_y1_given_do_x0 - p_x0y0) / p_x1y1 # 计算q值 q p_y0_given_x1 / (1 - p_y0_given_x1) lower_bound err cf upper_bound err q cf return max(0, lower_bound), min(1, upper_bound)2. 案例一戴口罩与感染率的因果分析2.1 数据生成与问题设定让我们首先模拟一个关于戴口罩与感染率关系的数据集。在现实中戴口罩行为可能与其他健康习惯相关混杂因素我们需要在模拟中体现这一点。import numpy as np import pandas as pd def generate_mask_data(sample_size10000): 生成戴口罩与感染率的模拟数据 np.random.seed(42) # 生成混杂因素健康意识水平影响戴口罩和感染概率 health_conscious np.random.normal(0, 1, sample_size) # 生成戴口罩行为受健康意识影响 mask_prob 1 / (1 np.exp(-health_conscious)) mask np.random.binomial(1, mask_prob) # 生成感染结果受戴口罩和健康意识影响 infection_prob 1 / (1 np.exp(-(-2*mask - health_conscious np.random.normal(0, 0.5, sample_size)))) infection np.random.binomial(1, infection_prob) return pd.DataFrame({ health_conscious: health_conscious, mask: mask, infection: infection }) mask_data generate_mask_data() print(mask_data.head())2.2 观测数据与实验数据对比在现实世界中我们通常只能获得观测数据但PN计算需要干预数据(P(Y|do(X)))。我们可以通过两种方式获得干预数据随机对照试验数据理想情况通过调整混杂因素估计干预效果from sklearn.linear_model import LogisticRegression def estimate_do_effect(data): 通过调整混杂因素估计P(Y1|do(X0)) # 拟合感染概率模型 model LogisticRegression() model.fit(data[[mask, health_conscious]], data[infection]) # 创建反事实场景所有人都没戴口罩 cf_data data.copy() cf_data[mask] 0 # 预测反事实感染概率 cf_infection_prob model.predict_proba(cf_data[[mask, health_conscious]])[:,1] return cf_infection_prob.mean() # 计算各项概率 p_y1 1 - mask_data[infection].mean() p_x1y1 ((mask_data[mask]1) (mask_data[infection]0)).mean() p_y1_given_do_x0 estimate_do_effect(mask_data) # 计算PN pn calculate_pn(p_y1, p_y1_given_do_x0, p_x1y1) print(f计算得到的PN值为: {pn:.3f})2.3 结果解读与应用假设我们得到PN0.65这意味着对于戴口罩且未感染的人群中有65%的人之所以未感染是因为戴了口罩。这个结果可以用于公共卫生政策制定量化戴口罩的防护效果个人行为决策帮助个人理解防护措施的重要性进一步研究设计识别需要控制的混杂因素注意在实际应用中确保所有相关混杂因素都被纳入模型至关重要。遗漏重要混杂因素会导致估计偏差。3. 案例二药品责任诉讼中的因果归因3.1 案例背景与数据准备考虑一起药品责任诉讼原告声称某种止痛药导致患者死亡而被告药厂提供的临床试验数据显示风险增加很小。我们需要通过反事实分析评估药品的实际责任。def generate_drug_lawsuit_data(): 生成药品诉讼案例的模拟数据 np.random.seed(123) # 试验数据随机分配无混杂 trial_size 10000 trial_drug np.random.binomial(1, 0.5, trial_size) trial_death np.random.binomial(1, 0.014 0.002*trial_drug, trial_size) # 观察数据存在自我选择偏差 obs_size 5000 health_status np.random.normal(0, 1, obs_size) # 未观测的健康状况 obs_drug np.random.binomial(1, 1/(1np.exp(-(-1 health_status))), obs_size) obs_death np.random.binomial(1, 0.01 0.02*obs_drug - 0.03*health_status, obs_size) trial_data pd.DataFrame({drug: trial_drug, death: trial_death, source: trial}) obs_data pd.DataFrame({drug: obs_drug, death: obs_death, source: obs}) return pd.concat([trial_data, obs_data]) drug_data generate_drug_lawsuit_data()3.2 传统风险率与PN对比传统方法计算过剩风险率def traditional_risk_ratio(data): 计算传统过剩风险率 p_yx data[data[drug]1][death].mean() p_yx_prime data[data[drug]0][death].mean() return (p_yx - p_yx_prime) / p_yx # 分别计算试验数据和观察数据的风险率 trial_rr traditional_risk_ratio(drug_data[drug_data[source]trial]) obs_rr traditional_risk_ratio(drug_data[drug_data[source]obs]) print(f试验数据计算的过剩风险率: {trial_rr:.3f}) print(f观察数据计算的过剩风险率: {obs_rr:.3f})计算PN值考虑单调性假设# 使用试验数据作为do(X0)的估计 p_y1_trial 1 - drug_data[drug_data[source]trial][death].mean() p_y1_given_do_x0 1 - drug_data[(drug_data[source]trial) (drug_data[drug]0)][death].mean() p_x1y1_obs ((drug_data[source]obs) (drug_data[drug]1) (drug_data[death]0)).mean() pn_drug calculate_pn(p_y1_trial, p_y1_given_do_x0, p_x1y1_obs) print(f药品案例的PN值: {pn_drug:.3f})3.3 法律与商业决策启示在这个案例中传统方法可能低估药品风险而PN分析揭示了更完整的因果图景。这对法律和商业决策有多方面启示法律责任评估PN值接近1表明药品很可能是死亡原因药品安全监测需要更全面的数据收集和分析方法临床试验设计考虑如何将观察数据纳入风险评估4. 案例三广告点击的归因分析4.1 多触点归因问题在数字营销中用户可能接触多个广告触点后转化。反事实推理可以帮助我们评估每个触点的贡献。def generate_ad_attribution_data(num_users10000): 生成广告归因模拟数据 np.random.seed(456) # 用户属性 user_activity np.random.exponential(1, num_users) # 用户活跃度 # 广告曝光受用户活跃度影响 ad1_prob 1 - 1/(1np.exp(-(user_activity-1))) ad2_prob 1 - 1/(1np.exp(-(user_activity-0.5))) ad1 np.random.binomial(1, ad1_prob) ad2 np.random.binomial(1, ad2_prob) # 转化受广告和用户活跃度影响 conv_prob 1/(1np.exp(-(-3 1.5*ad1 1.2*ad2 0.8*user_activity np.random.normal(0, 0.3, num_users)))) conversion np.random.binomial(1, conv_prob) return pd.DataFrame({ user_activity: user_activity, ad1: ad1, ad2: ad2, conversion: conversion }) ad_data generate_ad_attribution_data()4.2 反事实归因模型实现我们可以为每个广告触点计算其必要性概率def attribution_pn(data, ad_var): 计算单个广告触点的PN值 # 拟合转化模型 model LogisticRegression() model.fit(data[[ad1, ad2, user_activity]], data[conversion]) # 创建反事实场景该广告未展示 cf_data data.copy() cf_data[ad_var] 0 # 预测反事实转化概率 cf_conv_prob model.predict_proba(cf_data[[ad1, ad2, user_activity]])[:,1] # 计算各项概率 p_y1 data[conversion].mean() p_y1_given_do_x0 cf_conv_prob.mean() p_x1y1 (data[ad_var]1 data[conversion]1).mean() return calculate_pn(p_y1, p_y1_given_do_x0, p_x1y1) # 计算各广告的PN pn_ad1 attribution_pn(ad_data, ad1) pn_ad2 attribution_pn(ad_data, ad2) print(f广告1的PN值: {pn_ad1:.3f}) print(f广告2的PN值: {pn_ad2:.3f})4.3 营销预算分配建议基于PN值的归因分析可以指导营销决策预算优化将更多资源分配给PN值高的渠道跨渠道协同分析不同触点间的交互效应用户分群针对不同用户群体优化触点策略我们可以进一步将PN值转换为归因权重total_pn pn_ad1 pn_ad2 ad1_weight pn_ad1 / total_pn ad2_weight pn_ad2 / total_pn print(f广告1的归因权重: {ad1_weight:.2f}) print(f广告2的归因权重: {ad2_weight:.2f})5. 混杂因子处理与敏感性分析5.1 常见混杂因子识别在因果推断中混杂因子是同时影响处理变量和结果的变量。在前面的案例中戴口罩案例健康意识水平药品案例患者健康状况广告案例用户活跃度识别混杂因子的关键在于领域知识。以下是一个通用的混杂因子检查流程def check_confounders(data, treatment, outcome, potential_confs): 检查潜在混杂因子 参数: data: 数据集 treatment: 处理变量名 outcome: 结果变量名 potential_confs: 潜在混杂因子列表 返回: 与处理和结果都相关的混杂因子列表 confounders [] for var in potential_confs: # 检查与处理的关系 corr_treat data[[treatment, var]].corr().iloc[0,1] # 检查与结果的关系控制处理后 model LinearRegression() model.fit(data[[treatment, var]], data[outcome]) coef model.coef_[1] if abs(corr_treat) 0.1 and abs(coef) 0.1: confounders.append(var) return confounders5.2 基于倾向得分的校正倾向得分匹配是处理混杂的常用方法from sklearn.linear_model import LogisticRegression from sklearn.neighbors import NearestNeighbors def propensity_score_matching(data, treatment, outcome, confounders): 倾向得分匹配实现 # 计算倾向得分 ps_model LogisticRegression() ps_model.fit(data[confounders], data[treatment]) data[ps] ps_model.predict_proba(data[confounders])[:,1] # 对每个处理组个体寻找对照组匹配 treated data[data[treatment]1] control data[data[treatment]0] nbrs NearestNeighbors(n_neighbors1).fit(control[[ps]]) distances, indices nbrs.kneighbors(treated[[ps]]) # 创建匹配数据集 matched_control control.iloc[indices.flatten()].copy() matched_data pd.concat([treated, matched_control]) return matched_data5.3 敏感性分析实现评估结果对未观测混杂的敏感性def sensitivity_analysis(ate, se, confounder_strength): 因果效应估计的敏感性分析 参数: ate: 平均处理效应估计 se: 标准误 confounder_strength: 未观测混杂的影响强度 返回: 调整后的效应估计和p值 adjusted_ate ate - confounder_strength z_score adjusted_ate / se p_value 2 * (1 - stats.norm.cdf(abs(z_score))) return adjusted_ate, p_value6. 工程实践建议与常见陷阱6.1 数据需求与准备可靠的因果推断需要处理变量变异足够的变化以估计效应混杂因子测量尽可能测量所有相关混杂样本量充足特别是使用匹配方法时数据检查清单def causal_data_check(data, treatment, outcome, confounders): 因果分析数据质量检查 checks {} # 处理变量变异检查 treat_balance data[treatment].value_counts(normalizeTrue) checks[treatment_balance] min(treat_balance) 0.1 # 混杂因子与处理的关系 conf_treat_corr data[confounders].corrwith(data[treatment]) checks[conf_treat_corr] any(abs(conf_treat_corr) 0.1) # 混杂因子与结果的关系 model LinearRegression() model.fit(data[confounders [treatment]], data[outcome]) checks[conf_outcome_coef] any(abs(model.coef_[:-1]) 0.1) return checks6.2 模型选择考量不同场景下的模型选择建议场景特征推荐方法优点局限性线性关系少量混杂回归调整简单直观易于实现对模型误设敏感非线性关系倾向得分匹配/加权不假设线性关系需要重叠假设高维混杂双重机器学习处理高维数据计算复杂动态处理结构方程模型处理复杂因果结构需要强假设6.3 结果解释注意事项因果推断结果的合理解释需要明确反事实框架清楚说明比较的对象标注假设条件特别是识别假设量化不确定性提供置信区间或敏感性分析避免过度解读区分统计显著与实际显著结果报告模板def report_causal_result(estimate, ci_lower, ci_upper, assumptions): 因果效应结果报告模板 report f 因果效应分析结果 ------------------ 点估计值: {estimate:.3f} 95%置信区间: [{ci_lower:.3f}, {ci_upper:.3f}] 关键假设 1. 无未观测混杂: {assumptions[unobserved_confounding]} 2. 正值假设: {assumptions[positivity]} 3. 一致性: {assumptions[consistency]} 敏感性分析建议 - 未观测混杂需要解释{abs(estimate)*0.5:.2f}的效应量才能推翻结论 - 倾向得分重叠情况: {assumptions[overlap]} return report7. 扩展应用与前沿方向7.1 推荐系统中的反事实推理在推荐系统中反事实推理可以解决以下问题流行度偏差用户点击可能仅因物品流行曝光偏差用户只能看到被推荐的物品位置偏差用户倾向于点击靠前物品解决方案框架class CounterfactualRecommender: def __init__(self, base_model, propensity_model): self.base_model base_model # 基础推荐模型 self.propensity_model propensity_model # 倾向得分模型 def train(self, X, y, exposure): # 训练基础模型 self.base_model.fit(X, y) # 训练倾向模型曝光概率 self.propensity_model.fit(X, exposure) def predict(self, X, counterfactualFalse): # 基础预测 pred self.base_model.predict_proba(X)[:,1] if counterfactual: # 反事实调整 propensities self.propensity_model.predict_proba(X)[:,1] pred pred / np.clip(propensities, 1e-3, 1) return pred7.2 医疗决策中的个性化治疗反事实推理可以预测个体对不同治疗的反应def personalized_treatment_effect(data, features, treatment, outcome): 估计个体化处理效应 from sklearn.ensemble import RandomForestRegressor # 分开处理组和对照组 treated data[data[treatment]1] control data[data[treatment]0] # 分别在两组训练结果模型 model_treated RandomForestRegressor() model_treated.fit(treated[features], treated[outcome]) model_control RandomForestRegressor() model_control.fit(control[features], control[outcome]) # 为每个人预测反事实结果 data[y_treated] model_treated.predict(data[features]) data[y_control] model_control.predict(data[features]) # 计算个体处理效应 data[ite] data[y_treated] - data[y_control] return data7.3 与大语言模型的结合大语言模型(LLM)可以增强因果推断混杂因子识别利用领域知识建议潜在混杂因果图构建帮助构建合理的因果结构结果解释用自然语言解释复杂结果LLM辅助的混杂识别示例def llm_confounder_suggestion(variables, domain): 使用LLM API获取混杂因子建议 prompt f 在{domain}领域研究{variables[treatment]}对{variables[outcome]}的影响时 除了{variables[observed]}还有哪些重要的混杂变量需要考虑 请列出5个最关键的变量并简要说明理由。 # 这里应该是调用LLM API的代码 # 实际实现需要接入具体的LLM服务 response 1. 社会经济地位: 影响治疗选择和结果\n2. 家族病史: 可能影响结果\n3. 环境因素: 如居住地区\n4. 健康行为: 如吸烟饮酒\n5. 医疗可及性: 影响治疗质量 return response