因果推断实战：Python 计算 PN/PS 概率，从观测数据到反事实归因的 3 步流程
在数据分析领域，我们常常需要回答"如果当时做了不同的选择，结果会怎样"这类问题。这种思考方式被称为反事实推理，它是因果推断的核心工具之一。本文将带你用 Python 实现 PN（必要性概率）和 PS（充分性概率）的计算，通过 3 个步骤从观测数据得出因果结论。
1. 理解 PN 和 PS 概率
在因果推断中，PN（Probability of Necessity）和 PS（Probability of Sufficiency）是两个关键指标：
PN（必要性概率）：在事件发生的情况下，如果没有采取某个行动，事件不会发生的概率。例如："如果A先生没有戴口罩，他会被感染的概率是多少？"
PS（充分性概率）：在事件未发生的情况下，如果采取了某个行动，事件会发生的概率。例如："如果B先生戴了口罩，他不会感染的概率是多少？"
这两个概率的计算需要区分观测数据和干预数据：
# 观测数据示例（调查数据） observed_data = { '戴口罩': [1, 1, 0, 0, 1, 0, 1, 1], '感染': [0, 0, 1, 1, 0, 1, 0, 0] } # 干预数据示例（随机试验数据） experimental_data = { '强制不戴口罩': [0, 0, 1, 1, 0, 1, 0, 1], '感染': [1, 0, 1, 1, 0, 1, 0, 1] }2. 数据准备与预处理
计算 PN/PS 需要三类数据：
- 观测数据：自然状态下收集的数据（如调查问卷）
- 干预数据：通过随机对照试验获得的数据
- 混杂因子：影响因和果的其他变量
2.1 构建模拟数据集
我们使用 pandas 创建模拟数据集：
import pandas as pd import numpy as np # 设置随机种子保证可复现 np.random.seed(42) # 生成观测数据 n_samples = 1000 obs_df = pd.DataFrame({ '戴口罩': np.random.binomial(1, 0.6, n_samples), '年龄': np.random.normal(45, 15, n_samples).astype(int), '卫生习惯': np.random.beta(2, 5, n_samples) }) obs_df['感染'] = np.where( (obs_df['戴口罩'] == 0) & (obs_df['卫生习惯'] < 0.3), 1, np.random.binomial(1, 0.1, n_samples) ) # 生成干预数据（随机试验） exp_df = pd.DataFrame({ '强制不戴口罩': np.random.binomial(1, 0.5, n_samples), '年龄': np.random.normal(45, 15, n_samples).astype(int), '卫生习惯': np.random.beta(2, 5, n_samples) }) exp_df['感染'] = np.where( (exp_df['强制不戴口罩'] == 1) & (exp_df['卫生习惯'] < 0.4), 1, np.random.binomial(1, 0.15, n_samples) )2.2 数据探索分析
计算基础概率：
# 观测数据中的条件概率 p_y1_x1 = obs_df[obs_df['戴口罩']==1]['感染'].mean() # P(Y=1|X=1) p_y1_x0 = obs_df[obs_df['戴口罩']==0]['感染'].mean() # P(Y=1|X=0) # 干预数据中的概率 p_y1_do_x0 = exp_df[exp_df['强制不戴口罩']==1]['感染'].mean() # P(Y=1|do(X=0)) print(f"P(Y=1|X=1): {p_y1_x1:.3f}") print(f"P(Y=1|X=0): {p_y1_x0:.3f}") print(f"P(Y=1|do(X=0)): {p_y1_do_x0:.3f}")3. 计算 PN 和 PS 概率
3.1 单调性假设下的计算
当满足单调性假设（即干预不会产生反效果）时，PN 和 PS 的计算相对简单：
def calculate_pn_monotonic(p_y1, p_y1_do_x0, p_x1_y1): """计算单调性假设下的PN""" return (p_y1 - p_y1_do_x0) / p_x1_y1 def calculate_ps_monotonic(p_y0_do_x1, p_y0_x0, p_x0_y0): """计算单调性假设下的PS""" return (p_y0_do_x1 - p_y0_x0) / p_x0_y0 # 计算必要概率PN p_x1_y1 = len(obs_df[(obs_df['戴口罩']==1)&(obs_df['感染']==1)])/len(obs_df[obs_df['感染']==1]) pn = calculate_pn_monotonic(p_y1, p_y1_do_x0, p_x1_y1) # 计算充分概率PS (需要补充P(Y=0|do(X=1))的数据) # 这里假设我们有这个数据 p_y0_do_x1 = 0.85 # 假设值 p_x0_y0 = len(obs_df[(obs_df['戴口罩']==0)&(obs_df['感染']==0)])/len(obs_df[obs_df['感染']==0]) p_y0_x0 = 1 - p_y1_x0 ps = calculate_ps_monotonic(p_y0_do_x1, p_y0_x0, p_x0_y0) print(f"PN (单调性假设): {pn:.3f}") print(f"PS (单调性假设): {ps:.3f}")3.2 非单调情况下的上下界
当不满足单调性假设时，我们计算 PN 的上下界：
def calculate_pn_bounds(p_y1, p_y1_x1, p_y1_x0, p_y1_do_x0, p_x1_y1): """计算PN的上下界""" # 计算ERR (Excess Risk Ratio) err = (p_y1_x1 - p_y1_x0) / p_y1_x1 # 计算混杂因子CF (Confounding Factor) cf = (p_y1_x0 - p_y1_do_x0) / p_x1_y1 # 计算q值 q = (1 - p_y1_x1) / p_y1_x1 # 下界 lb = err + cf # 上界 ub = err + q + cf return max(0, lb), min(1, ub) # 计算P(Y=1) p_y1 = obs_df['感染'].mean() # 计算PN的上下界 lb, ub = calculate_pn_bounds(p_y1, p_y1_x1, p_y1_x0, p_y1_do_x0, p_x1_y1) print(f"PN 下界: {lb:.3f}") print(f"PN 上界: {ub:.3f}")3.3 完整计算流程封装
我们将上述步骤封装成一个完整的类：
class CounterfactualAnalyzer: def __init__(self, observed_data, experimental_data): self.obs_df = observed_data self.exp_df = experimental_data def calculate_probabilities(self): """计算所有需要的概率""" # 观测数据概率 self.p_y1_x1 = self.obs_df[self.obs_df['戴口罩']==1]['感染'].mean() self.p_y1_x0 = self.obs_df[self.obs_df['戴口罩']==0]['感染'].mean() self.p_x1_y1 = len(self.obs_df[(self.obs_df['戴口罩']==1)&(self.obs_df['感染']==1)])/len(self.obs_df[self.obs_df['感染']==1]) self.p_x0_y0 = len(self.obs_df[(self.obs_df['戴口罩']==0)&(self.obs_df['感染']==0)])/len(self.obs_df[self.obs_df['感染']==0]) # 干预数据概率 self.p_y1_do_x0 = self.exp_df[self.exp_df['强制不戴口罩']==1]['感染'].mean() # 其他概率 self.p_y1 = self.obs_df['感染'].mean() self.p_y0_x0 = 1 - self.p_y1_x0 def calculate_pn_monotonic(self): """计算单调性假设下的PN""" return (self.p_y1 - self.p_y1_do_x0) / self.p_x1_y1 def calculate_ps_monotonic(self, p_y0_do_x1): """计算单调性假设下的PS""" return (p_y0_do_x1 - self.p_y0_x0) / self.p_x0_y0 def calculate_pn_bounds(self): """计算PN的上下界""" err = (self.p_y1_x1 - self.p_y1_x0) / self.p_y1_x1 cf = (self.p_y1_x0 - self.p_y1_do_x0) / self.p_x1_y1 q = (1 - self.p_y1_x1) / self.p_y1_x1 lb = err + cf ub = err + q + cf return max(0, lb), min(1, ub) def analyze(self, p_y0_do_x1=None): """执行完整分析""" self.calculate_probabilities() results = { 'PN_monotonic': self.calculate_pn_monotonic(), 'PN_lower_bound': None, 'PN_upper_bound': None } if p_y0_do_x1 is not None: results['PS_monotonic'] = self.calculate_ps_monotonic(p_y0_do_x1) lb, ub = self.calculate_pn_bounds() results['PN_lower_bound'] = lb results['PN_upper_bound'] = ub return results # 使用示例 analyzer = CounterfactualAnalyzer(obs_df, exp_df) results = analyzer.analyze(p_y0_do_x1=0.85) # 打印结果 print("分析结果:") for k, v in results.items(): print(f"{k}: {v:.3f}" if isinstance(v, (int, float)) else f"{k}: {v}")4. 结果可视化与解释
4.1 关键指标可视化
import matplotlib.pyplot as plt # 准备数据 labels = ['P(Y=1|X=1)', 'P(Y=1|X=0)', 'P(Y=1|do(X=0))'] values = [p_y1_x1, p_y1_x0, p_y1_do_x0] # 创建条形图 plt.figure(figsize=(10, 6)) bars = plt.bar(labels, values, color=['blue', 'orange', 'green']) plt.ylabel('概率值') plt.title('关键概率指标比较') plt.ylim(0, 0.5) # 在每个条形上添加数值标签 for bar in bars: height = bar.get_height() plt.text(bar.get_x() + bar.get_width()/2., height, f'{height:.3f}', ha='center', va='bottom') plt.show()4.2 PN 和 PS 的概率解释
根据计算结果：
PN 值：如果 PN=0.7，意味着在戴口罩且未感染的人群中，有70%的人如果不戴口罩会被感染。这表明戴口罩对预防感染有较强的必要性。
PS 值：如果 PS=0.6，意味着在不戴口罩且感染的人群中，有60%的人如果戴了口罩可以避免感染。这表明戴口罩对预防感染有中等程度的充分性。
4.3 敏感性分析
我们可以测试不同假设对结果的影响：
# 测试不同p_y0_do_x1值对PS的影响 ps_values = [] p_y0_do_x1_range = np.linspace(0.7, 0.95, 10) for p in p_y0_do_x1_range: ps = analyzer.calculate_ps_monotonic(p) ps_values.append(ps) # 绘制敏感性分析图 plt.figure(figsize=(10, 6)) plt.plot(p_y0_do_x1_range, ps_values, marker='o') plt.xlabel('P(Y=0|do(X=1))') plt.ylabel('PS 值') plt.title('PS 对 P(Y=0|do(X=1)) 的敏感性分析') plt.grid(True) plt.show()5. 实际应用中的注意事项
在实际应用中，有几个关键点需要考虑：
- 数据质量：确保观测数据和干预数据来自相同的总体
- 混杂控制：尽可能测量和调整所有重要的混杂因素
- 假设验证：检验单调性假设是否合理
- 不确定性量化：使用自助法等方法计算置信区间
# 自助法计算PN的置信区间 def bootstrap_pn(data, n_bootstrap=1000): pn_samples = [] n = len(data) for _ in range(n_bootstrap): sample = data.sample(n, replace=True) analyzer = CounterfactualAnalyzer(sample, exp_df) # 假设exp_df不变 res = analyzer.analyze() pn_samples.append(res['PN_monotonic']) return np.percentile(pn_samples, [2.5, 97.5]) # 计算PN的95%置信区间 ci = bootstrap_pn(obs_df) print(f"PN的95%置信区间: [{ci[0]:.3f}, {ci[1]:.3f}]")6. 扩展应用：医药案例
让我们模拟一个医药案例，分析某种药物是否导致患者死亡：
# 医药案例数据 drug_data = { '服药': [1, 1, 0, 0, 1, 0, 1, 1, 0, 1], '死亡': [0, 0, 1, 1, 0, 1, 0, 0, 1, 1], '年龄': [45, 60, 75, 68, 52, 70, 48, 55, 80, 65] } drug_df = pd.DataFrame(drug_data) # 干预数据（临床试验） drug_exp = { '强制服药': [1, 1, 0, 0, 1, 0, 1, 0], '死亡': [0, 1, 1, 0, 0, 1, 0, 1], '年龄': [50, 65, 70, 55, 60, 72, 58, 68] } drug_exp_df = pd.DataFrame(drug_exp) # 分析 drug_analyzer = CounterfactualAnalyzer(drug_df, drug_exp_df) drug_results = drug_analyzer.analyze(p_y0_do_x1=0.9) # 假设P(Y=0|do(X=1))=0.9 print("\n医药案例分析结果:") for k, v in drug_results.items(): print(f"{k}: {v:.3f}" if isinstance(v, (int, float)) else f"{k}: {v}")7. 总结与最佳实践
通过本文的 Python 实现，我们掌握了计算 PN 和 PS 概率的完整流程。在实际项目中，建议遵循以下最佳实践：
数据收集阶段：
- 确保有高质量的观测数据和干预数据
- 尽可能收集潜在的混杂因素
分析阶段：
- 先进行探索性数据分析，理解数据分布
- 检验关键假设（如单调性）
- 计算点估计和区间估计
结果解释阶段：
- 结合领域知识解释结果
- 进行敏感性分析，评估结果稳健性
- 明确结论的局限性
# 最佳实践示例：检查单调性假设 def check_monotonicity(obs_df, exp_df): """检查数据是否满足单调性假设""" p_y1_x1 = obs_df[obs_df['戴口罩']==1]['感染'].mean() p_y1_x0 = obs_df[obs_df['戴口罩']==0]['感染'].mean() p_y1_do_x0 = exp_df[exp_df['强制不戴口罩']==1]['感染'].mean() # 简单检查：干预效果不应比观测效果差 return p_y1_do_x0 <= p_y1_x0 is_monotonic = check_monotonicity(obs_df, exp_df) print(f"数据是否满足单调性假设: {is_monotonic}")
