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

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

在数据科学领域，我们常常需要回答"如果...会怎样"这类问题。比如：如果用户没有点击这个广告，他们的购买概率会降低多少？如果患者服用了另一种药物，康复时间会缩短多少？这类问题本质上都是在探讨因果关系，而反事实推理正是解决这类问题的利器。

反事实推理的核心思想是构建一个与现实相反的假设场景，通过对比实际结果与假设结果来推断因果关系。这种思维方式不仅存在于学术研究中，也渗透在我们的日常决策里——当你想着"如果昨天早点出门就不会迟到"时，就已经在无意识地进行反事实推理了。

1. 反事实推理基础与核心公式

1.1 必要性概率(PN)与充分性概率(PS)

在反事实推理中，两个核心概念是必要性概率(PN)和充分性概率(PS)。让我们通过一个简单的例子来理解它们：

假设X表示"戴口罩"(X=1)或"不戴口罩"(X=0)，Y表示"感染"(Y=0)或"未感染"(Y=1)。对于一位戴口罩且未感染的人(A先生)，PN回答的问题是："如果A先生没有戴口罩，他被感染的概率是多少？"数学表达式为：

PN = P(Y_{X=0}=0 | X=1, Y=1)

而对于一位不戴口罩且被感染的人(B先生)，PS回答的问题是："如果B先生戴了口罩，他不会被感染的概率是多少？"数学表达式为：

PS = P(Y_{X=1}=1 | X=0, Y=0)

1.2 单调性假设下的PN计算

在实际应用中，PN的计算常引入单调性假设：即干预的结果不会比不干预更差(Y_X=1 ≥ Y_X=0)。在单调性假设下，PN可以表示为：

def calculate_pn(p_y1, p_y1_given_do_x0, p_x1y1): """ 计算单调性假设下的必要性概率(PN) 参数: p_y1: P(Y=1) p_y1_given_do_x0: P(Y=1|do(X=0)) p_x1y1: P(X=1,Y=1) 返回: PN值 """ err = (p_y1 - p_y1_given_do_x0) / p_x1y1 return err

这个公式包含两个关键部分：

1. 过剩风险率(ERR)：反映干预与不干预的风险差异
2. 混杂因子(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(Y=1) p_y1_given_do_x0: P(Y=1|do(X=0)) p_x0y0: P(X=0,Y=0) p_x1y1: P(X=1,Y=1) p_y0_given_x1: P(Y=0|X=1) 返回: (下界, 上界) """ # 计算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_size=10000): """生成戴口罩与感染率的模拟数据""" 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)))。我们可以通过两种方式获得干预数据：

1. 随机对照试验数据（理想情况）
2. 通过调整混杂因素估计干预效果

from sklearn.linear_model import LogisticRegression def estimate_do_effect(data): """通过调整混杂因素估计P(Y=1|do(X=0))""" # 拟合感染概率模型 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 结果解读与应用

假设我们得到PN=0.65，这意味着对于戴口罩且未感染的人群中，有65%的人之所以未感染是因为戴了口罩。这个结果可以用于：

1. 公共卫生政策制定：量化戴口罩的防护效果
2. 个人行为决策：帮助个人理解防护措施的重要性
3. 进一步研究设计：识别需要控制的混杂因素

注意：在实际应用中，确保所有相关混杂因素都被纳入模型至关重要。遗漏重要混杂因素会导致估计偏差。

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/(1+np.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(X=0)的估计 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分析揭示了更完整的因果图景。这对法律和商业决策有多方面启示：

1. 法律责任评估：PN值接近1表明药品很可能是死亡原因
2. 药品安全监测：需要更全面的数据收集和分析方法
3. 临床试验设计：考虑如何将观察数据纳入风险评估

4. 案例三：广告点击的归因分析

4.1 多触点归因问题

在数字营销中，用户可能接触多个广告触点后转化。反事实推理可以帮助我们评估每个触点的贡献。

def generate_ad_attribution_data(num_users=10000): """生成广告归因模拟数据""" np.random.seed(456) # 用户属性 user_activity = np.random.exponential(1, num_users) # 用户活跃度 # 广告曝光（受用户活跃度影响） ad1_prob = 1 - 1/(1+np.exp(-(user_activity-1))) ad2_prob = 1 - 1/(1+np.exp(-(user_activity-0.5))) ad1 = np.random.binomial(1, ad1_prob) ad2 = np.random.binomial(1, ad2_prob) # 转化（受广告和用户活跃度影响） conv_prob = 1/(1+np.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值的归因分析可以指导营销决策：

1. 预算优化：将更多资源分配给PN值高的渠道
2. 跨渠道协同：分析不同触点间的交互效应
3. 用户分群：针对不同用户群体优化触点策略

我们可以进一步将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 常见混杂因子识别

在因果推断中，混杂因子是同时影响处理变量和结果的变量。在前面的案例中：

1. 戴口罩案例：健康意识水平
2. 药品案例：患者健康状况
3. 广告案例：用户活跃度

识别混杂因子的关键在于领域知识。以下是一个通用的混杂因子检查流程：

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 confounders

5.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_neighbors=1).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_data

5.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_value

6. 工程实践建议与常见陷阱

6.1 数据需求与准备

可靠的因果推断需要：

1. 处理变量变异：足够的变化以估计效应
2. 混杂因子测量：尽可能测量所有相关混杂
3. 样本量充足：特别是使用匹配方法时

数据检查清单：

def causal_data_check(data, treatment, outcome, confounders): """因果分析数据质量检查""" checks = {} # 处理变量变异检查 treat_balance = data[treatment].value_counts(normalize=True) 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 checks

6.2 模型选择考量

不同场景下的模型选择建议：

6.3 结果解释注意事项

因果推断结果的合理解释需要：

1. 明确反事实框架：清楚说明比较的对象
2. 标注假设条件：特别是识别假设
3. 量化不确定性：提供置信区间或敏感性分析
4. 避免过度解读：区分统计显著与实际显著

结果报告模板：

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 report

7. 扩展应用与前沿方向

7.1 推荐系统中的反事实推理

在推荐系统中，反事实推理可以解决以下问题：

1. 流行度偏差：用户点击可能仅因物品流行
2. 曝光偏差：用户只能看到被推荐的物品
3. 位置偏差：用户倾向于点击靠前物品

解决方案框架：

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, counterfactual=False): # 基础预测 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 pred

7.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 data

7.3 与大语言模型的结合

大语言模型(LLM)可以增强因果推断：

1. 混杂因子识别：利用领域知识建议潜在混杂
2. 因果图构建：帮助构建合理的因果结构
3. 结果解释：用自然语言解释复杂结果

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