前门准则与后门准则：混淆变量的调整策略

I. 引言：混淆变量的挑战与因果识别的困境

混淆变量的本质问题

在观察性研究中，混淆变量是同时影响处理变量和结果变量的第三变量，导致观察到的关联不能反映真实的因果效应。混淆变量问题的严重性在于：

• 普遍存在：几乎所有的观察性研究都存在潜在的混淆变量
• 隐蔽性强：许多重要的混淆变量往往难以观测或测量
• 偏误方向不确定：既可能高估也可能低估真实效应

经典案例：在研究吸烟对肺癌的影响时，基因因素可能同时影响吸烟行为和肺癌易感性。如果忽略这一混淆，我们可能错误地将吸烟的效应归因于基因，或者反之。

因果识别的基本框架

为了从观察数据中识别因果效应，我们需要满足某种形式的"可识别性条件"。在潜在结果框架下，这通常意味着条件可忽略性。而在因果图框架中，前门准则和后门准则提供了图形化的识别条件。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import networkx as nx
import statsmodels.api as sm
from sklearn.linear_model import LinearRegression
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

print("前门准则与后门准则：混淆变量的调整策略")
print("=" * 50)

II. 后门准则：处理可观测混淆的标准方法

后门准则的直观理解

后门准则的核心思想很简单：通过阻断所有"后门路径"来隔离真正的因果效应。后门路径是指从处理变量X到结果变量Y的非因果路径，这些路径通常通过混淆变量连接。

数学定义：在因果图G中，一组变量Z满足后门准则相对于(X,Y)，如果：

1. Z阻断所有X和Y之间的后门路径
2. Z不包含X的任何后代节点

# 后门准则的可视化示例
def create_backdoor_example():
    """创建后门准则的典型示例"""
    fig, axes = plt.subplots(1, 2, figsize=(15, 6))
    
    # 示例1：单一混淆变量
    G1 = nx.DiGraph()
    G1.add_edges_from([('X', 'Y'), ('U', 'X'), ('U', 'Y')])
    pos1 = {'U': (1, 1), 'X': (0, 0), 'Y': (2, 0)}
    
    nx.draw(G1, pos1, with_labels=True, node_size=2000, 
            node_color=['lightblue', 'lightcoral', 'lightgreen'], 
            arrowsize=20, ax=axes[0])
    axes[0].set_title('后门路径示例: X ← U → Y\n调整U可阻断后门路径')
    
    # 示例2：多重混淆
    G2 = nx.DiGraph()
    G2.add_edges_from([('X', 'Y'), ('U1', 'X'), ('U1', 'Y'), ('U2', 'X'), ('U2', 'Y')])
    pos2 = {'U1': (0.5, 1), 'U2': (1.5, 1), 'X': (0, 0), 'Y': (2, 0)}
    
    nx.draw(G2, pos2, with_labels=True, node_size=2000,
            node_color=['lightblue', 'lightblue', 'lightcoral', 'lightgreen'],
            arrowsize=20, ax=axes[1])
    axes[1].set_title('多重混淆: 需要调整{U1, U2}\n阻断所有后门路径')
    
    plt.tight_layout()
    plt.show()

create_backdoor_example()

后门准则的数学表达

当满足后门准则时，因果效应可以通过调整公式识别：

P(Y|do(X=x)) = Σᵤ P(Y|X=x, U=u)P(U=u)

其中U是满足后门准则的变量集合。

# 后门准则的数值验证
def simulate_backdoor_scenario(n=10000):
    """模拟后门准则的基本场景"""
    np.random.seed(2024)
    
    # 混淆变量U
    U = np.random.normal(0, 1, n)
    
    # 处理变量X受U影响
    X = 0.5 * U + np.random.normal(0, 0.5, n)
    
    # 结果变量Y受X和U影响
    Y = 2.0 * X + 1.5 * U + np.random.normal(0, 1, n)
    
    data = pd.DataFrame({'U': U, 'X': X, 'Y': Y})
    
    return data

# 生成数据
backdoor_data = simulate_backdoor_scenario()

print("后门准则验证数据描述:")
print(f"样本量: {len(backdoor_data)}")
print(f"X的均值: {backdoor_data['X'].mean():.3f}")
print(f"Y的均值: {backdoor_data['Y'].mean():.3f}")
print(f"U的均值: {backdoor_data['U'].mean():.3f}")

# 比较不同调整策略
def compare_adjustment_methods(data):
    """比较不同调整方法的估计效果"""
    
    results = []
    
    # 1. 无调整（有偏估计）
    model_naive = sm.OLS(data['Y'], sm.add_constant(data['X'])).fit()
    naive_estimate = model_naive.params['X']
    results.append({
        'method': '无调整',
        'estimate': naive_estimate,
        'bias': naive_estimate - 2.0,  # 真实效应为2.0
        'bias_percentage': (naive_estimate - 2.0) / 2.0 * 100
    })
    
    # 2. 后门调整（调整U）
    model_adjusted = sm.OLS(data['Y'], sm.add_constant(data[['X', 'U']])).fit()
    adjusted_estimate = model_adjusted.params['X']
    results.append({
        'method': '后门调整(U)',
        'estimate': adjusted_estimate,
        'bias': adjusted_estimate - 2.0,
        'bias_percentage': (adjusted_estimate - 2.0) / 2.0 * 100
    })
    
    # 3. 错误调整（调整Y的后代，如果有的话）
    # 这里简化为调整一个随机变量
    data['W'] = np.random.normal(0, 1, len(data))  # 模拟无关变量
    model_wrong = sm.OLS(data['Y'], sm.add_constant(data[['X', 'W']])).fit()
    wrong_estimate = model_wrong.params['X']
    results.append({
        'method': '错误调整(W)',
        'estimate': wrong_estimate,
        'bias': wrong_estimate - 2.0,
        'bias_percentage': (wrong_estimate - 2.0) / 2.0 * 100
    })
    
    return pd.DataFrame(results)

# 比较调整方法
adjustment_comparison = compare_adjustment_methods(backdoor_data)
print("\n不同调整方法的比较:")
print(adjustment_comparison.round(4))

# 可视化结果
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 估计值比较
methods = adjustment_comparison['method']
estimates = adjustment_comparison['estimate']
biases = np.abs(adjustment_comparison['bias'])

colors = ['red', 'green', 'orange']
axes[0].bar(methods, estimates, color=colors, alpha=0.7)
axes[0].axhline(y=2.0, color='blue', linestyle='--', linewidth=2, label='真实效应 (2.0)')
axes[0].set_ylabel('因果效应估计')
axes[0].set_title('不同调整方法的估计比较')
axes[0].legend()
axes[0].tick_params(axis='x', rotation=45)

# 偏差比较
axes[1].bar(methods, biases, color=colors, alpha=0.7)
axes[1].set_ylabel('绝对偏差')
axes[1].set_title('估计偏差比较')
axes[1].tick_params(axis='x', rotation=45)

plt.tight_layout()
plt.show()

后门准则的应用条件

后门准则的成功应用依赖于以下条件：

III. 前门准则：处理未观测混淆的创新方法

前门准则的直觉理解

当存在未观测的混淆变量时，后门准则无法直接应用。前门准则通过利用中介变量提供了另一种识别策略：

核心思想：即使X和Y之间存在未观测的混淆，如果存在一个完全中介变量M，且满足特定条件，我们仍然可以识别因果效应。

前门准则的正式定义

在因果图G中，一组变量M满足前门准则相对于(X,Y)，如果：

1. M完全中介X对Y的效应（所有从X到Y的有向路径都经过M）
2. X和M之间没有未观测的混淆
3. M和Y之间所有未观测的混淆都被X阻断

# 前门准则的可视化示例
def create_frontdoor_example():
    """创建前门准则的典型示例"""
    fig, axes = plt.subplots(1, 2, figsize=(15, 6))
    
    # 示例1：经典吸烟-焦油-癌症例子
    G1 = nx.DiGraph()
    G1.add_edges_from([('Smoking', 'Tar'), ('Tar', 'Cancer')])
    # 未观测混淆U
    G1.add_edges_from([('U', 'Smoking'), ('U', 'Cancer')])
    
    pos1 = {'U': (1, 1), 'Smoking': (0, 0), 'Tar': (1, 0), 'Cancer': (2, 0)}
    
    nx.draw(G1, pos1, with_labels=True, node_size=2000, 
            node_color=['lightgray', 'lightcoral', 'lightyellow', 'lightgreen'],
            arrowsize=20, ax=axes[0])
    axes[0].set_title('前门准则示例: 吸烟→焦油→癌症\n存在未观测混淆U')
    
    # 示例2：更复杂的前门路径
    G2 = nx.DiGraph()
    G2.add_edges_from([('X', 'M1'), ('M1', 'M2'), ('M2', 'Y')])
    G2.add_edges_from([('U1', 'X'), ('U1', 'Y'), ('U2', 'M2'), ('U2', 'Y')])
    
    pos2 = {'U1': (0.5, 1), 'U2': (1.5, 1), 'X': (0, 0), 
            'M1': (1, 0), 'M2': (2, 0), 'Y': (3, 0)}
    
    nx.draw(G2, pos2, with_labels=True, node_size=2000,
            node_color=['lightgray', 'lightgray', 'lightcoral', 
                       'lightyellow', 'lightyellow', 'lightgreen'],
            arrowsize=20, ax=axes[1])
    axes[1].set_title('复杂前门路径: X→M1→M2→Y\n存在多重未观测混淆')
    
    plt.tight_layout()
    plt.show()

create_frontdoor_example()

前门准则的估计公式

当满足前门准则时，因果效应可以通过以下公式识别：

P(Y|do(X=x)) = Σₘ P(m|X=x) Σₓ′ P(Y|X=x′, m)P(X=x′)

这个公式通过两个阶段估计因果效应：

1. X对M的效应
2. M对Y的效应（控制X）

# 前门准则的数值实现
def frontdoor_estimation(data, x_var, m_var, y_var):
    """
    实现前门准则的因果效应估计
    """
    # 第一阶段: 估计X对M的效应
    stage1 = sm.OLS(data[m_var], sm.add_constant(data[x_var])).fit()
    alpha = stage1.params[x_var]  # X -> M的效应
    
    # 第二阶段: 估计M对Y的效应，控制X
    stage2 = sm.OLS(data[y_var], sm.add_constant(data[[x_var, m_var]])).fit()
    beta = stage2.params[m_var]   # M -> Y的效应（控制X）
    
    # 前门估计: alpha * beta
    frontdoor_estimate = alpha * beta
    
    # 计算标准误（使用Delta方法）
    alpha_se = stage1.bse[x_var]
    beta_se = stage2.bse[m_var]
    frontdoor_se = np.sqrt((beta**2 * alpha_se**2) + (alpha**2 * beta_se**2))
    
    return {
        'estimate': frontdoor_estimate,
        'std_error': frontdoor_se,
        'ci_lower': frontdoor_estimate - 1.96 * frontdoor_se,
        'ci_upper': frontdoor_estimate + 1.96 * frontdoor_se,
        'alpha': alpha,
        'beta': beta
    }

# 模拟前门准则场景
def simulate_frontdoor_scenario(n=10000):
    """模拟前门准则的应用场景"""
    np.random.seed(2024)
    
    # 未观测混淆变量U
    U = np.random.normal(0, 1, n)
    
    # 处理变量X受U影响
    X = 0.6 * U + np.random.normal(0, 0.5, n)
    
    # 中介变量M完全由X决定（满足前门条件1）
    M = 0.8 * X + np.random.normal(0, 0.3, n)
    
    # 结果变量Y受M和U影响（满足前门条件3）
    Y = 1.5 * M + 1.2 * U + np.random.normal(0, 1, n)
    
    data = pd.DataFrame({'X': X, 'M': M, 'Y': Y, 'U': U})
    
    # 真实因果效应: X->M->Y = 0.8 * 1.5 = 1.2
    true_effect = 0.8 * 1.5
    
    return data, true_effect

# 生成数据并应用前门准则
frontdoor_data, true_effect = simulate_frontdoor_scenario()

print("前门准则验证数据描述:")
print(f"样本量: {len(frontdoor_data)}")
print(f"X的均值: {frontdoor_data['X'].mean():.3f}")
print(f"M的均值: {frontdoor_data['M'].mean():.3f}")
print(f"Y的均值: {frontdoor_data['Y'].mean():.3f}")
print(f"真实因果效应: {true_effect:.3f}")

# 应用前门准则估计
frontdoor_result = frontdoor_estimation(frontdoor_data, 'X', 'M', 'Y')

print("\n前门准则估计结果:")
print(f"因果效应估计: {frontdoor_result['estimate']:.4f}")
print(f"标准误: {frontdoor_result['std_error']:.4f}")
print(f"95%置信区间: [{frontdoor_result['ci_lower']:.4f}, {frontdoor_result['ci_upper']:.4f}]")
print(f"X→M效应 (α): {frontdoor_result['alpha']:.4f}")
print(f"M→Y效应 (β): {frontdoor_result['beta']:.4f}")

# 比较不同方法
def compare_frontdoor_backdoor(data, true_effect):
    """比较前门准则和后门准则的表现"""
    
    results = []
    
    # 1. 无调整（有偏）
    naive_model = sm.OLS(data['Y'], sm.add_constant(data['X'])).fit()
    results.append({
        'method': '无调整',
        'estimate': naive_model.params['X'],
        'bias': naive_model.params['X'] - true_effect
    })
    
    # 2. 后门调整（需要U，但U通常未观测）
    # 这里我们假设U可观测作为基准
    backdoor_model = sm.OLS(data['Y'], sm.add_constant(data[['X', 'U']])).fit()
    results.append({
        'method': '后门调整(U)',
        'estimate': backdoor_model.params['X'],
        'bias': backdoor_model.params['X'] - true_effect
    })
    
    # 3. 前门准则
    frontdoor_result = frontdoor_estimation(data, 'X', 'M', 'Y')
    results.append({
        'method': '前门准则',
        'estimate': frontdoor_result['estimate'],
        'bias': frontdoor_result['estimate'] - true_effect
    })
    
    return pd.DataFrame(results)

# 执行比较
method_comparison = compare_frontdoor_backdoor(frontdoor_data, true_effect)
print("\n方法比较结果:")
print(method_comparison.round(4))

# 可视化比较
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

methods = method_comparison['method']
estimates = method_comparison['estimate']
biases = np.abs(method_comparison['bias'])

colors = ['red', 'green', 'blue']
axes[0].bar(methods, estimates, color=colors, alpha=0.7)
axes[0].axhline(y=true_effect, color='black', linestyle='--', 
                linewidth=2, label=f'真实效应 ({true_effect})')
axes[0].set_ylabel('因果效应估计')
axes[0].set_title('不同方法的估计比较')
axes[0].legend()
axes[0].tick_params(axis='x', rotation=45)

axes[1].bar(methods, biases, color=colors, alpha=0.7)
axes[1].set_ylabel('绝对偏差')
axes[1].set_title('估计偏差比较')
axes[1].tick_params(axis='x', rotation=45)

plt.tight_layout()
plt.show()

前门准则的应用条件

前门准则的成功应用需要严格的条件：

IV. 实例分析：教育回报的因果识别

研究背景与问题设定

教育对收入的影响是社会科学中的经典问题。传统方法面临严重混淆问题：

• 能力偏差：高能力个体选择更多教育且获得更高收入
• 家庭背景：优越家庭提供教育机会和就业优势
• 动机因素：个人动机影响教育选择和职业成就

构建因果图

# 构建教育回报的因果图
def create_education_causal_graph():
    """创建教育回报的因果图"""
    G = nx.DiGraph()
    
    # 添加变量和边
    # 未观测变量
    G.add_node('Ability', color='lightgray')      # 能力
    G.add_node('Family', color='lightgray')       # 家庭背景
    G.add_node('Motivation', color='lightgray')   # 动机
    
    # 处理变量
    G.add_node('Education', color='lightcoral')   # 教育
    
    # 中介变量（前门路径）
    G.add_node('Skills', color='lightyellow')     # 工作技能
    G.add_node('Network', color='lightyellow')    # 社会网络
    G.add_node('Signaling', color='lightyellow')  # 信号效应
    
    # 结果变量
    G.add_node('Income', color='lightgreen')      # 收入
    
    # 添加因果边
    # 混淆路径
    G.add_edges_from([
        ('Ability', 'Education'), ('Ability', 'Income'),
        ('Family', 'Education'), ('Family', 'Income'), 
        ('Motivation', 'Education'), ('Motivation', 'Income')
    ])
    
    # 前门路径（教育通过多个机制影响收入）
    G.add_edges_from([
        ('Education', 'Skills'), ('Skills', 'Income'),
        ('Education', 'Network'), ('Network', 'Income'),
        ('Education', 'Signaling'), ('Signaling', 'Income')
    ])
    
    return G

# 可视化教育回报因果图
education_graph = create_education_causal_graph()

plt.figure(figsize=(12, 8))
pos = {
    'Ability': (0, 2), 'Family': (1, 2), 'Motivation': (2, 2),
    'Education': (1, 1),
    'Skills': (0, 0), 'Network': (1, 0), 'Signaling': (2, 0),
    'Income': (1, -1)
}

node_colors = [education_graph.nodes[node]['color'] for node in education_graph.nodes()]
nx.draw(education_graph, pos, with_labels=True, node_color=node_colors, 
        node_size=2500, arrowsize=20, font_size=10, alpha=0.8)

plt.title("教育回报的因果图\n灰色: 未观测混淆, 红色: 处理变量, 黄色: 中介变量, 绿色: 结果变量")
plt.show()

# 分析因果识别策略
print("教育回报的因果识别分析:")
print("=" * 50)

print("后门准则挑战:")
print("- 主要混淆变量(能力、家庭背景、动机)通常难以完全观测")
print("- 即使测量，也可能存在测量误差")
print("- 需要收集大量协变量数据")

print("\n前门准则机会:")
print("- 教育通过技能、网络、信号等机制影响收入")
print("- 这些中介变量相对容易观测和测量") 
print("- 可以分别估计各机制的效应")

print("\n推荐策略:")
print("1. 后门调整: 控制所有可观测混淆变量")
print("2. 前门估计: 利用中介变量处理未观测混淆")
print("3. 工具变量: 寻找外生教育变异来源")
print("4. 多重验证: 比较不同方法的结果")

数据生成与效应估计

# 生成教育回报的模拟数据
def generate_education_data(n=10000):
    """生成教育回报的模拟数据"""
    np.random.seed(2024)
    
    # 未观测混淆变量
    ability = np.random.normal(0, 1, n)
    family_bg = np.random.normal(0, 1, n) 
    motivation = np.random.normal(0, 1, n)
    
    # 教育决策（受混淆变量影响）
    education_prob = 1 / (1 + np.exp(-(
        -2.0 + 0.6 * ability + 0.5 * family_bg + 0.4 * motivation
    )))
    education = np.random.binomial(1, education_prob)  # 0=高中, 1=大学
    
    # 中介变量（前门路径）
    skills = 0.7 * education + 0.3 * ability + np.random.normal(0, 0.5, n)
    network = 0.6 * education + 0.4 * family_bg + np.random.normal(0, 0.5, n)
    signaling = 0.8 * education + np.random.normal(0, 0.3, n)
    
    # 收入生成
    base_income = 30000
    income = (
        base_income +
        5000 * education +  # 直接效应
        3000 * ability +
        2000 * family_bg + 
        1500 * motivation +
        4000 * skills +     # 技能路径
        3000 * network +    # 网络路径
        2000 * signaling +  # 信号路径
        np.random.normal(0, 4000, n)
    )
    
    data = pd.DataFrame({
        'education': education,
        'skills': skills, 'network': network, 'signaling': signaling,
        'income': income,
        'ability': ability, 'family_bg': family_bg, 'motivation': motivation
    })
    
    # 真实因果效应
    true_effect = 5000  # 直接效应
    
    return data, true_effect

# 生成数据
edu_data, true_effect_edu = generate_education_data()

print("教育回报模拟数据描述:")
print(f"样本量: {len(edu_data)}")
print(f"大学教育比例: {edu_data['education'].mean():.3f}")
print(f"大学组平均收入: ${edu_data[edu_data['education']==1]['income'].mean():.2f}")
print(f"高中组平均收入: ${edu_data[edu_data['education']==0]['income'].mean():.2f}")
print(f"简单收入差异: ${edu_data[edu_data['education']==1]['income'].mean() - edu_data[edu_data['education']==0]['income'].mean():.2f}")
print(f"真实因果效应: ${true_effect_edu:.2f}")

# 实施不同的因果识别策略
def implement_education_strategies(data, true_effect):
    """实施教育回报的不同因果识别策略"""
    
    results = []
    
    # 1. 无调整
    naive = sm.OLS(data['income'], sm.add_constant(data['education'])).fit()
    results.append({
        'method': '无调整',
        'estimate': naive.params['education'],
        'bias': naive.params['education'] - true_effect,
        'ci_width': (naive.conf_int().loc['education', 1] - naive.conf_int().loc['education', 0])
    })
    
    # 2. 后门调整（部分可观测变量）
    # 假设我们只能观测到家庭背景
    backdoor_partial = sm.OLS(data['income'], 
                             sm.add_constant(data[['education', 'family_bg']])).fit()
    results.append({
        'method': '后门调整(家庭背景)',
        'estimate': backdoor_partial.params['education'],
        'bias': backdoor_partial.params['education'] - true_effect,
        'ci_width': (backdoor_partial.conf_int().loc['education', 1] - 
                    backdoor_partial.conf_int().loc['education', 0])
    })
    
    # 3. 后门调整（更多可观测变量）
    # 假设我们能观测家庭背景和动机
    backdoor_more = sm.OLS(data['income'], 
                          sm.add_constant(data[['education', 'family_bg', 'motivation']])).fit()
    results.append({
        'method': '后门调整(家庭+动机)',
        'estimate': backdoor_more.params['education'],
        'bias': backdoor_more.params['education'] - true_effect,
        'ci_width': (backdoor_more.conf_int().loc['education', 1] - 
                    backdoor_more.conf_int().loc['education', 0])
    })
    
    # 4. 前门准则（通过技能路径）
    frontdoor_skills = frontdoor_estimation(data, 'education', 'skills', 'income')
    results.append({
        'method': '前门准则(技能)',
        'estimate': frontdoor_skills['estimate'],
        'bias': frontdoor_skills['estimate'] - true_effect,
        'ci_width': frontdoor_skills['ci_upper'] - frontdoor_skills['ci_lower']
    })
    
    # 5. 理想情况：所有混淆可观测（基准）
    ideal = sm.OLS(data['income'], 
                  sm.add_constant(data[['education', 'ability', 'family_bg', 'motivation']])).fit()
    results.append({
        'method': '理想调整(所有混淆)',
        'estimate': ideal.params['education'],
        'bias': ideal.params['education'] - true_effect,
        'ci_width': (ideal.conf_int().loc['education', 1] - 
                    ideal.conf_int().loc['education', 0])
    })
    
    return pd.DataFrame(results)

# 实施策略比较
strategy_results = implement_education_strategies(edu_data, true_effect_edu)
print("\n教育回报的不同识别策略结果:")
print(strategy_results.round(2))

# 可视化策略比较
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

methods = strategy_results['method']
estimates = strategy_results['estimate']
biases = np.abs(strategy_results['bias'])
ci_widths = strategy_results['ci_width']

# 估计值比较
axes[0].bar(methods, estimates, alpha=0.7, color=['red', 'orange', 'yellow', 'lightgreen', 'green'])
axes[0].axhline(y=true_effect_edu, color='blue', linestyle='--', 
                linewidth=2, label=f'真实效应 (${true_effect_edu})')
axes[0].set_ylabel('教育回报估计 ($)')
axes[0].set_title('不同方法的因果效应估计')
axes[0].legend()
axes[0].tick_params(axis='x', rotation=45)

# 偏差比较
axes[1].bar(methods, biases, alpha=0.7, color=['red', 'orange', 'yellow', 'lightgreen', 'green'])
axes[1].set_ylabel('绝对偏差 ($)')
axes[1].set_title('估计偏差比较')
axes[1].tick_params(axis='x', rotation=45)

# 置信区间宽度比较
axes[2].bar(methods, ci_widths, alpha=0.7, color=['red', 'orange', 'yellow', 'lightgreen', 'green'])
axes[2].set_ylabel('95%置信区间宽度 ($)')
axes[2].set_title('估计不确定性比较')
axes[2].tick_params(axis='x', rotation=45)

plt.tight_layout()
plt.show()

实例分析洞察

通过系统的因果识别分析，我们获得以下重要发现：

1. 混淆变量的严重影响：未调整的估计严重高估教育回报
2. 部分调整的局限性：即使控制部分混淆，估计仍有明显偏误
3. 前门准则的价值：在存在未观测混淆时提供可行的识别策略
4. 方法互补性：不同方法相互验证提高结论可靠性

V. 准则比较与选择指南

后门准则 vs 前门准则

两种准则在理论基础、应用条件和实践可行性上存在重要差异：

# 准则比较分析
def compare_criteria():
    """系统比较前后门准则"""
    
    comparison_data = {
        'Criterion': ['后门准则', '前门准则'],
        '核心思想': [
            '阻断所有后门路径', 
            '利用完全中介路径'
        ],
        '主要条件': [
            '所有混淆变量可观测', 
            '存在完全中介且满足前门条件'
        ],
        '优势': [
            '直观易懂，应用广泛', 
            '能处理未观测混淆'
        ],
        '局限': [
            '需要所有混淆变量可观测', 
            '前门条件在现实中较难满足'
        ],
        '适用场景': [
            '随机化实验的观察性模仿', 
            '存在明确中介机制的复杂场景'
        ],
        '估计方法': [
            '回归调整、匹配、加权', 
            '两阶段乘积估计'
        ]
    }
    
    return pd.DataFrame(comparison_data)

# 显示准则比较
criteria_comparison = compare_criteria()
print("后门准则与前门准则的系统比较:")
print("=" * 60)
print(criteria_comparison.to_string(index=False))

# 可视化准则选择流程
def plot_criteria_selection_flow():
    """绘制准则选择流程图"""
    fig, ax = plt.subplots(1, 1, figsize=(10, 8))
    
    # 创建流程图元素
    elements = [
        ('开始: 因果识别问题', 'start', 0.5, 0.9),
        ('所有混淆变量可观测?', 'decision', 0.5, 0.7),
        ('应用后门准则', 'process', 0.3, 0.5),
        ('存在满足前门条件的中介?', 'decision', 0.5, 0.5),
        ('应用前门准则', 'process', 0.7, 0.3),
        ('考虑工具变量或其他方法', 'process', 0.5, 0.1),
        ('结束', 'end', 0.5, -0.1)
    ]
    
    # 绘制元素
    for text, etype, x, y in elements:
        if etype == 'start':
            ax.add_patch(plt.Circle((x, y), 0.05, fill=True, color='green', alpha=0.7))
        elif etype == 'end':
            ax.add_patch(plt.Circle((x, y), 0.05, fill=True, color='red', alpha=0.7))
        elif etype == 'decision':
            ax.add_patch(plt.Rectangle((x-0.15, y-0.05), 0.3, 0.1, fill=True, 
                                     color='lightblue', alpha=0.7))
        else:  # process
            ax.add_patch(plt.Rectangle((x-0.15, y-0.05), 0.3, 0.1, fill=True, 
                                     color='lightyellow', alpha=0.7))
        
        ax.text(x, y, text, ha='center', va='center', fontsize=9, 
               bbox=dict(boxstyle="round,pad=0.3", facecolor='white', alpha=0.8))
    
    # 绘制连接线
    arrows = [
        ((0.5, 0.85), (0.5, 0.75)),  # 开始 -> 决策1
        ((0.5, 0.65), (0.3, 0.55)),  # 是 -> 后门
        ((0.3, 0.45), (0.5, 0.35)),  # 后门 -> 结束前
        ((0.5, 0.65), (0.7, 0.55)),  # 否 -> 决策2
        ((0.7, 0.45), (0.7, 0.35)),  # 是 -> 前门
        ((0.7, 0.25), (0.5, 0.15)),  # 前门 -> 其他
        ((0.5, 0.65), (0.5, 0.55)),  # 否 -> 决策2
        ((0.5, 0.45), (0.5, 0.35)),  # 否 -> 其他
        ((0.5, 0.05), (0.5, -0.05))  # 其他 -> 结束
    ]
    
    for (x1, y1), (x2, y2) in arrows:
        ax.arrow(x1, y1, x2-x1, y2-y1, head_width=0.02, head_length=0.03, 
                fc='black', ec='black', alpha=0.6)
    
    ax.set_xlim(0, 1)
    ax.set_ylim(-0.2, 1)
    ax.set_aspect('equal')
    ax.axis('off')
    ax.set_title('前后门准则选择流程图', fontsize=14, pad=20)
    
    plt.tight_layout()
    plt.show()

plot_criteria_selection_flow()

选择指南与实践建议

基于理论分析和实践经验，我们总结出以下选择指南：

VI. 实践中的挑战与解决方案

未观测混淆的敏感性分析

即使使用前门准则，我们仍需要评估结论对未观测混淆的稳健性。

# 敏感性分析实现
def sensitivity_analysis_continuous(data, x_var, y_var, rho_range=np.linspace(0, 0.8, 9)):
    """
    连续处理变量的敏感性分析
    """
    results = []
    
    # 基准估计（无调整）
    baseline_model = sm.OLS(data[y_var], sm.add_constant(data[x_var])).fit()
    baseline_effect = baseline_model.params[x_var]
    
    for rho in rho_range:
        if rho == 0:
            # 无混淆情况
            adjusted_effect = baseline_effect
        else:
            # 模拟未观测混淆的影响
            # 简化方法：根据相关性调整估计
            # 实际敏感性分析更复杂
            adjustment_factor = 1 / (1 - rho**2)
            adjusted_effect = baseline_effect / adjustment_factor
        
        results.append({
            'confounding_strength': rho,
            'estimated_effect': adjusted_effect,
            'bias': adjusted_effect - true_effect_edu if 'true_effect_edu' in locals() else np.nan
        })
    
    return pd.DataFrame(results)

# 执行敏感性分析
sensitivity_results = sensitivity_analysis_continuous(edu_data, 'education', 'income')

print("未观测混淆的敏感性分析:")
print("=" * 50)
print(sensitivity_results.round(3))

# 可视化敏感性分析
plt.figure(figsize=(10, 6))
plt.plot(sensitivity_results['confounding_strength'], 
         sensitivity_results['estimated_effect'], 'bo-', linewidth=2, markersize=6)
if 'true_effect_edu' in locals():
    plt.axhline(y=true_effect_edu, color='red', linestyle='--', 
                linewidth=2, label=f'真实效应 (${true_effect_edu})')
plt.axhline(y=sensitivity_results['estimated_effect'].iloc[0], color='green', 
            linestyle='--', linewidth=2, label='基准估计')
plt.xlabel('未观测混淆的强度 (ρ)')
plt.ylabel('调整后的因果效应估计')
plt.title('因果效应估计对未观测混淆的敏感性')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

模型设定稳健性检验

不同的模型设定可能影响因果效应估计，需要进行稳健性检验。

# 模型稳健性检验
def robustness_checks(data, x_var, y_var, adjustment_sets):
    """
    进行模型设定稳健性检验
    """
    robustness_results = []
    
    for i, adj_set in enumerate(adjustment_sets):
        if adj_set:  # 有调整变量
            formula = f"{y_var} ~ {x_var} + " + " + ".join(adj_set)
        else:  # 无调整变量
            formula = f"{y_var} ~ {x_var}"
        
        model = sm.OLS.from_formula(formula, data=data).fit()
        
        robustness_results.append({
            'model': f"模型{i+1}",
            'adjustment_set': str(adj_set) if adj_set else "无",
            'estimate': model.params[x_var],
            'std_error': model.bse[x_var],
            'r_squared': model.rsquared
        })
    
    return pd.DataFrame(robustness_results)

# 定义不同的调整集
adjustment_sets = [
    [],  # 无调整
    ['family_bg'],  # 仅家庭背景
    ['family_bg', 'motivation'],  # 家庭背景+动机
    ['skills'],  # 仅技能（前门路径）
    ['family_bg', 'skills']  # 混合策略
]

# 执行稳健性检验
robustness_results = robustness_checks(edu_data, 'education', 'income', adjustment_sets)
print("模型设定稳健性检验:")
print(robustness_results.round(4))

# 可视化稳健性结果
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

models = robustness_results['model']
estimates = robustness_results['estimate']
std_errors = robustness_results['std_error']

# 估计值比较
axes[0].errorbar(models, estimates, yerr=std_errors, fmt='o', capsize=5, 
                capthick=2, markersize=8, color='blue', alpha=0.7)
if 'true_effect_edu' in locals():
    axes[0].axhline(y=true_effect_edu, color='red', linestyle='--', 
                   linewidth=2, label='真实效应')
axes[0].set_ylabel('因果效应估计')
axes[0].set_title('不同模型设定的估计比较')
axes[0].tick_params(axis='x', rotation=45)
if 'true_effect_edu' in locals():
    axes[0].legend()

# R²比较
r_squared = robustness_results['r_squared']
axes[1].bar(models, r_squared, alpha=0.7, color='green')
axes[1].set_ylabel('R²')
axes[1].set_title('模型拟合优度比较')
axes[1].tick_params(axis='x', rotation=45)

plt.tight_layout()
plt.show()

实践挑战的解决方案

在实际应用中，我们面临多种挑战，需要相应的解决方案：