2.5 do-算子与前门准则

当不可观测混淆存在时,因果效应仍然可以识别

do-算子:干预的数学化

do-算子定义

Pearl (2000) 的核心贡献:

含义: 如果我们强制取值为,切断所有指向的箭头,Y的分布是什么?

do-算子 vs 条件概率

这是因果推断中最重要的不等式。

: 观察到后,的条件分布（被动观察）
: 强制设定后,的分布（主动干预）

例子: 吸烟与肺癌

DAG:

观察相关:

可能部分由基因混淆。

因果效应:

通过干预(强制吸烟),切断"基因→吸烟"箭头,得到纯因果效应。

do-算子的图操作

算法

do-算子对应 DAG 上的一个简单操作:

在DAG中,删除所有指向的箭头
在修改后的DAG(记为)中计算

直觉: 干预相当于把从一个受其他变量影响的内生变量,变成一个由外部力量决定的外生变量。

公式

在中。

其中是在原始DAG中的父节点集合。

与后门调整的关系

当后门准则成立时,do-算子可以通过后门调整公式计算:

这意味着:干预分布可以从观察数据中计算出来——不需要真的做实验!

前门准则(Frontdoor Criterion)

为什么需要前门准则?

后门准则要求我们能够观测到足够多的混淆变量来阻断后门路径。但在很多实际情况中,混淆变量是不可观测的。

前门准则: 当存在中介变量时的识别策略——即使存在不可观测混淆!

定义

前门准则成立的条件:

条件1: 截获所有从到的有向路径

条件2: 和之间没有后门路径

条件3: 所有到的后门路径被阻断

前门调整公式

如果满足前门准则:

直觉: 分两步走

先估计对的因果效应: （没有后门路径,直接可估计）
再估计对的因果效应: 通过对边际化,用阻断的后门路径
两步组合得到对的因果效应

经典例子:吸烟与肺癌

DAG

DAG:

         基因(U,不可观测)
        /              \
       ↓                ↓
    吸烟(X) → 焦油(M) → 肺癌(Y)

为什么后门准则不适用?

问题: 基因不可观测,无法用后门准则!

后门路径无法通过控制来阻断。

为什么前门准则适用?

检查: 焦油满足前门准则!

条件1: (焦油)截获了(吸烟)到(肺癌)的所有有向路径 ✓

唯一的有向路径是 ,在路径上

条件2: 和之间没有后门路径 ✓

没有其他变量同时影响吸烟和焦油（基因不直接影响焦油）

条件3: 到的后门路径被阻断 ✓

控制后,这条后门路径被阻断

估计步骤

步骤1: 估计:

这可以直接从数据中估计（因为吸烟和焦油之间没有后门路径）。

步骤2: 估计:

步骤3: 边际化:

意义

即使有不可观测混淆,仍可识别因果效应! 这是 Pearl 因果推断理论最令人惊叹的成果之一。

中介分析(Mediation Analysis)

直接效应 vs 间接效应

设定:

    X → M → Y
     \     /
      ↓   ↓

总效应(Total Effect):

直接效应(Direct Effect):

固定中介,X对Y的效应。

间接效应(Indirect Effect):

通过中介的效应。

自然直接/间接效应

Baron & Kenny (1986) 的经典框架:

自然直接效应(NDE):

自然间接效应(NIE):

其中是当时的值。

分解:

识别条件

需要假设:

1. 无混淆: 外生

2. 无中介-结果混淆: 没有影响和的共同原因

3. 无处理-中介混淆

Python实现:中介分析

python

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

def mediation_analysis(X, M, Y):
    """
    中介分析

    参数:
    ------
    X: (N,) 处理
    M: (N,) 中介变量
    Y: (N,) 结果

    返回:
    ------
    effects: dict, 包含总效应、直接效应、间接效应
    """
    N = len(X)

    # 步骤1:总效应 (Y ~ X)
    model_total = LinearRegression()
    model_total.fit(X.reshape(-1, 1), Y)
    total_effect = model_total.coef_[0]

    # 步骤2:X对M的效应 (M ~ X)
    model_mediator = LinearRegression()
    model_mediator.fit(X.reshape(-1, 1), M)
    a = model_mediator.coef_[0]  # X → M

    # 步骤3:控制M后,X对Y的直接效应 (Y ~ X + M)
    model_direct = LinearRegression()
    model_direct.fit(np.column_stack([X, M]), Y)
    c_prime = model_direct.coef_[0]  # 直接效应
    b = model_direct.coef_[1]  # M → Y

    # 间接效应
    indirect_effect = a * b

    # Sobel检验
    se_indirect = np.sqrt(b**2 * np.var(M) + a**2 * np.var(Y))
    z_score = indirect_effect / se_indirect
    p_value = 2 * (1 - stats.norm.cdf(abs(z_score)))

    results = {
        'Total Effect': total_effect,
        'Direct Effect': c_prime,
        'Indirect Effect': indirect_effect,
        'Proportion Mediated': indirect_effect / total_effect if total_effect != 0 else 0,
        'Sobel z': z_score,
        'p-value': p_value
    }

    return results

# 示例
np.random.seed(42)
N = 1000
X = np.random.binomial(1, 0.5, N)
M = 0.5*X + np.random.normal(0, 0.5, N)  # X → M
Y = 0.3*X + 0.4*M + np.random.normal(0, 0.5, N)  # X → Y, M → Y

results = mediation_analysis(X, M, Y)
for key, value in results.items():
    print(f"{key}: {value:.4f}")

️ DoWhy库:因果推断统一框架

DoWhy的四步框架

Microsoft开发的DoWhy统一了Pearl和Rubin:

步骤1: 建模(Model) - 构造因果图

步骤2: 识别(Identify) - 使用图准则识别因果效应

步骤3: 估计(Estimate) - 选择估计方法

步骤4: 反驳(Refute) - 敏感性分析

完整示例

python

import numpy as np
import pandas as pd
import dowhy
from dowhy import CausalModel
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']

# =================================================================
# 第1步:生成数据
# =================================================================

np.random.seed(123)
N = 2000

# 混淆变量
Z1 = np.random.normal(0, 1, N)
Z2 = np.random.binomial(1, 0.5, N)

# 处理(受Z影响)
logit_D = -0.5 + 0.8*Z1 + 0.5*Z2
p_D = 1 / (1 + np.exp(-logit_D))
D = np.random.binomial(1, p_D, N)

# 结果(受D和Z影响)
Y = 2*D + 1.5*Z1 + 0.8*Z2 + np.random.normal(0, 1, N)

# 创建DataFrame
data = pd.DataFrame({
    'D': D,
    'Y': Y,
    'Z1': Z1,
    'Z2': Z2
})

print(f"数据概览:")
print(data.describe())
print(f"\n简单均值差(有偏): {data[data['D']==1]['Y'].mean() - data[data['D']==0]['Y'].mean():.3f}")
print(f"真实ATE: 2.0")

# =================================================================
# 第2步:构造因果模型
# =================================================================

# 方法1:指定DAG
model = CausalModel(
    data=data,
    treatment='D',
    outcome='Y',
    graph="""
    digraph {
        Z1 -> D;
        Z2 -> D;
        Z1 -> Y;
        Z2 -> Y;
        D -> Y;
    }
    """
)

# 可视化DAG
model.view_model()
plt.title('因果图', fontsize=14, fontweight='bold')
plt.show()

# 方法2:自动识别混淆
model_auto = CausalModel(
    data=data,
    treatment='D',
    outcome='Y',
    common_causes=['Z1', 'Z2']
)

# =================================================================
# 第3步:识别因果效应
# =================================================================

identified_estimand = model.identify_effect(proceed_when_unidentifiable=True)
print("\n识别的因果效应:")
print(identified_estimand)

# =================================================================
# 第4步:估计因果效应(多种方法)
# =================================================================

# 4a. 后门调整(线性回归)
estimate_backdoor = model.estimate_effect(
    identified_estimand,
    method_name="backdoor.linear_regression"
)
print(f"\n后门调整(回归): ATE = {estimate_backdoor.value:.3f}")

# 4b. 倾向得分匹配
estimate_psm = model.estimate_effect(
    identified_estimand,
    method_name="backdoor.propensity_score_matching"
)
print(f"倾向得分匹配: ATE = {estimate_psm.value:.3f}")

# 4c. 倾向得分加权
estimate_ipw = model.estimate_effect(
    identified_estimand,
    method_name="backdoor.propensity_score_weighting"
)
print(f"倾向得分加权: ATE = {estimate_ipw.value:.3f}")

# =================================================================
# 第5步:反驳/稳健性检验
# =================================================================

print("\n稳健性检验:")

# 5a. 随机共同原因
refute_random = model.refute_estimate(
    identified_estimand,
    estimate_backdoor,
    method_name="random_common_cause"
)
print(f"随机共同原因检验: p = {refute_random.refutation_result['p_value']:.3f}")

# 5b. 安慰剂处理
refute_placebo = model.refute_estimate(
    identified_estimand,
    estimate_backdoor,
    method_name="placebo_treatment_refuter"
)
print(f"安慰剂处理检验: 新估计 = {refute_placebo.estimated_effect:.3f}")

# 5c. 数据子集验证
refute_subset = model.refute_estimate(
    identified_estimand,
    estimate_backdoor,
    method_name="data_subset_refuter",
    subset_fraction=0.8
)
print(f"子集验证: 新估计 = {refute_subset.estimated_effect:.3f}")

# =================================================================
# 第6步:CATE估计(可选)
# =================================================================

from dowhy.causal_estimators.econml import Econml

# 使用EconML的DML估计CATE
estimate_dml = model.estimate_effect(
    identified_estimand,
    method_name="backdoor.econml.dml.DML",
    method_params={
        "init_params": {
            'model_y': GradientBoostingRegressor(),
            'model_t': GradientBoostingRegressor()
        },
        "fit_params": {}
    }
)
print(f"\nDML估计: ATE = {estimate_dml.value:.3f}")

# 异质性分析
if hasattr(estimate_dml, 'estimator'):
    cate = estimate_dml.estimator.effect(data[['Z1', 'Z2']].values)
    print(f"CATE范围: [{cate.min():.3f}, {cate.max():.3f}]")

# =================================================================
# 第7步:可视化
# =================================================================

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 子图1:不同方法的估计
methods = ['后门调整', 'PSM', 'IPW', 'DML']
estimates = [
    estimate_backdoor.value,
    estimate_psm.value,
    estimate_ipw.value,
    estimate_dml.value
]

ax1 = axes[0]
ax1.barh(methods, estimates, color='steelblue', alpha=0.7)
ax1.axvline(x=2.0, color='red', linestyle='--', linewidth=2, label='真实ATE')
ax1.set_xlabel('估计的ATE', fontsize=12)
ax1.set_title('不同方法比较', fontsize=14, fontweight='bold')
ax1.legend()
ax1.grid(axis='x', alpha=0.3)

# 子图2:CATE分布
if 'cate' in locals():
    ax2 = axes[1]
    ax2.hist(cate, bins=50, alpha=0.7, color='green', edgecolor='black')
    ax2.axvline(x=cate.mean(), color='red', linestyle='--', linewidth=2,
                label=f'平均CATE = {cate.mean():.3f}')
    ax2.set_xlabel('CATE', fontsize=12)
    ax2.set_ylabel('频数', fontsize=12)
    ax2.set_title('异质性处理效应分布', fontsize=14, fontweight='bold')
    ax2.legend()
    ax2.grid(axis='y', alpha=0.3)

plt.tight_layout()
plt.show()

print("\nDoWhy分析完成!")

DoWhy的优势

1. 统一框架: 融合Pearl(图)和Rubin(潜在结果)

2. 明确假设: 通过DAG清晰展示识别假设

3. 多种方法: 支持后门、前门、IV等

4. 自动化: 自动识别调整集

5. 稳健性: 内置多种反驳检验

️ 常见陷阱

陷阱1:错误的DAG

问题: 如果DAG画错,所有推断都错!

缓解:

领域知识
敏感性分析
比较多个DAG

陷阱2:忽略时间顺序

问题: DAG必须反映时间因果顺序

错误:

收入 → 教育  (不合理!)

正确:

教育 → 收入

陷阱3:前门准则的适用范围

注意: 前门准则需要严格的条件,在实际应用中相对少见。大多数情况下,后门准则是更常用的识别策略。

本章总结

核心要点

1. 因果阶梯:

Level 1: 关联
Level 2: 干预
Level 3: 反事实

2. DAG的三大工具:

d-分离: 识别独立性
后门准则: 通过调整识别因果效应
前门准则: 利用中介识别

3. do-算子: 形式化干预

4. 中介分析: 分解总效应 = 直接效应 + 间接效应

5. DoWhy: 统一Pearl和Rubin的实用工具

关键公式

do-算子:

后门调整:

前门调整:

中介分解:

Python工具

dowhy: 因果推断统一框架
causalgraphicalmodels: DAG操作
networkx: 图论算法
pgmpy: 概率图模型

上一节: 2.4 碰撞偏差与控制变量选择 | 下一节: Module 3 Python与回归基础

因果图:从图形到因果的认知桥梁!

2.5 do-算子与前门准则 ​

do-算子:干预的数学化 ​

do-算子定义 ​

do-算子 vs 条件概率 ​

例子: 吸烟与肺癌 ​

do-算子的图操作 ​

算法 ​

公式 ​

与后门调整的关系 ​

前门准则(Frontdoor Criterion) ​

为什么需要前门准则? ​

定义 ​

前门调整公式 ​

经典例子:吸烟与肺癌 ​

DAG ​

为什么后门准则不适用? ​

为什么前门准则适用? ​

估计步骤 ​

意义 ​

中介分析(Mediation Analysis) ​

直接效应 vs 间接效应 ​

自然直接/间接效应 ​

识别条件 ​

Python实现:中介分析 ​

️ DoWhy库:因果推断统一框架 ​

DoWhy的四步框架 ​

完整示例 ​

DoWhy的优势 ​

️ 常见陷阱 ​

陷阱1:错误的DAG ​

陷阱2:忽略时间顺序 ​

陷阱3:前门准则的适用范围 ​

本章总结 ​

核心要点 ​

关键公式 ​

Python工具 ​