中介效应分析：探究因果机制的有效方法

I. 引言：从总效应到机制探索

为什么需要中介效应分析

在传统的因果推断中，我们主要关注处理变量（X）对结果变量（Y）的总效应。然而，这种"黑箱"式的分析存在明显局限：

• 政策制定：只知道教育提高收入不够，还需要知道通过什么途径
• 理论发展：机制理解推动学科理论进步
• 干预优化：识别关键中介路径，提高干预效率

经典案例：在研究教育对收入的影响时，我们不仅想知道教育是否提高收入，更希望了解这种影响是通过提高技能、扩大社交网络，还是通过信号传递实现的。这种机制理解对教育政策制定至关重要。

中介效应分析的基本框架

中介效应分析的核心思想是将总效应分解为直接效应和间接效应：

• 总效应：X对Y的总影响
• 直接效应：X不通过中介变量直接影响Y的路径
• 间接效应：X通过中介变量M影响Y的路径

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.api as sm
from statsmodels.formula.api import ols
from sklearn.linear_model import LinearRegression
from sklearn.utils import resample
import warnings
warnings.filterwarnings('ignore')

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

print("中介效应分析：探究因果机制的有效方法")
print("=" * 50)

II. 中介效应的理论基础

中介效应模型的发展历程

中介效应分析的概念最早可以追溯到20世纪初的路径分析，但现代中介分析框架主要由Baron和Kenny（1986）系统提出，后来经过Preacher和Hayes（2004，2008）等人的发展，特别是引入Bootstrap方法后变得更加完善。

三类主要效应定义

在中介效应分析中，我们主要关注三类效应：

中介效应的识别假设

要获得无偏的中介效应估计，需要满足以下关键假设：

1. 无混淆假设：没有未观测的混淆变量影响X→Y、X→M、M→Y关系
2. 一致性假设：观测到的处理水平就是实际接受的处理
3. 模型设定正确：中介模型的形式设定正确
4. 无交互作用：处理变量和中介变量没有交互作用（可放松）

# 模拟基本的中介效应数据
np.random.seed(2024)
n = 2000

# 生成自变量X（如教育年限）
X = np.random.normal(12, 2, n)  # 平均教育年限12年

# 中介变量M（如工作技能）
a = 0.5  # X对M的效应
M = 10 + a * X + np.random.normal(0, 1, n)

# 结果变量Y（如收入）
c_prime = 2.0  # X对Y的直接效应
b = 1.5       # M对Y的效应
Y = 20 + c_prime * X + b * M + np.random.normal(0, 3, n)

# 计算总效应
total_effect = c_prime + a * b

# 创建数据框
mediation_data = pd.DataFrame({
    'education': X,
    'skills': M,
    'income': Y
})

print("中介效应模拟数据描述:")
print(f"样本量: {n}")
print(f"教育年限均值: {X.mean():.2f} ± {X.std():.2f}")
print(f"工作技能均值: {M.mean():.2f} ± {M.std():.2f}")
print(f"收入均值: {Y.mean():.2f} ± {Y.std():.2f}")
print(f"\n真实参数:")
print(f"X→M效应 (a): {a:.2f}")
print(f"M→Y效应 (b): {b:.2f}")
print(f"直接效应 (c'): {c_prime:.2f}")
print(f"间接效应 (ab): {a*b:.2f}")
print(f"总效应 (c): {total_effect:.2f}")

# 可视化基本关系
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# X与Y的关系
axes[0,0].scatter(mediation_data['education'], mediation_data['income'], alpha=0.5)
z = np.polyfit(mediation_data['education'], mediation_data['income'], 1)
p = np.poly1d(z)
axes[0,0].plot(mediation_data['education'], p(mediation_data['education']), "r-", linewidth=2)
axes[0,0].set_xlabel('教育年限 (X)')
axes[0,0].set_ylabel('收入 (Y)')
axes[0,0].set_title('X → Y: 总效应')
axes[0,0].text(0.05, 0.95, f'斜率: {z[0]:.2f}', transform=axes[0,0].transAxes, 
               fontsize=12, verticalalignment='top',
               bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

# X与M的关系
axes[0,1].scatter(mediation_data['education'], mediation_data['skills'], alpha=0.5)
z = np.polyfit(mediation_data['education'], mediation_data['skills'], 1)
p = np.poly1d(z)
axes[0,1].plot(mediation_data['education'], p(mediation_data['education']), "r-", linewidth=2)
axes[0,1].set_xlabel('教育年限 (X)')
axes[0,1].set_ylabel('工作技能 (M)')
axes[0,1].set_title('X → M: 路径a')
axes[0,1].text(0.05, 0.95, f'斜率: {z[0]:.2f}', transform=axes[0,1].transAxes, 
               fontsize=12, verticalalignment='top',
               bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

# M与Y的关系（控制X）
# 这里使用残差图展示偏相关
model_m_y = ols('income ~ education + skills', data=mediation_data).fit()
mediation_data['income_resid'] = model_m_y.resid
model_m = ols('skills ~ education', data=mediation_data).fit()
mediation_data['skills_resid'] = model_m.resid

axes[1,0].scatter(mediation_data['skills_resid'], mediation_data['income_resid'], alpha=0.5)
z = np.polyfit(mediation_data['skills_resid'], mediation_data['income_resid'], 1)
p = np.poly1d(z)
axes[1,0].plot(mediation_data['skills_resid'], p(mediation_data['skills_resid']), "r-", linewidth=2)
axes[1,0].set_xlabel('技能残差 (控制X)')
axes[1,0].set_ylabel('收入残差 (控制X)')
axes[1,0].set_title('M → Y: 路径b (控制X)')
axes[1,0].text(0.05, 0.95, f'斜率: {z[0]:.2f}', transform=axes[1,0].transAxes, 
               fontsize=12, verticalalignment='top',
               bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

# 效应分解
effects = ['总效应', '直接效应', '间接效应']
values = [total_effect, c_prime, a*b]
colors = ['lightblue', 'lightgreen', 'lightcoral']

axes[1,1].bar(effects, values, color=colors, alpha=0.7)
axes[1,1].set_ylabel('效应大小')
axes[1,1].set_title('效应分解')
for i, v in enumerate(values):
    axes[1,1].text(i, v + 0.1, f'{v:.2f}', ha='center', va='bottom')

plt.tight_layout()
plt.show()

中介模型的理论基础

中介效应分析建立在因果推断的潜在结果框架上，近年来发展的因果中介分析框架提供了更严谨的理论基础。

III. 中介效应的检验方法

逐步检验法（Baron & Kenny方法）

这是最传统的中介效应检验方法，包含四个步骤：

1. 检验X对Y的总效应（路径c）
2. 检验X对M的效应（路径a）
3. 检验M对Y的效应（路径b）
4. 检验加入M后X对Y的直接效应（路径c’）

# 逐步检验法实现
def baron_kenny_test(data, x_var, m_var, y_var):
    """
    执行Baron & Kenny逐步检验法
    """
    results = {}
    
    # 步骤1: X -> Y (总效应c)
    model1 = ols(f'{y_var} ~ {x_var}', data=data).fit()
    c = model1.params[x_var]
    c_pvalue = model1.pvalues[x_var]
    results['step1'] = {
        'model': 'Y ~ X',
        'coefficient': c,
        'pvalue': c_pvalue,
        'significant': c_pvalue < 0.05
    }
    
    # 步骤2: X -> M (路径a)
    model2 = ols(f'{m_var} ~ {x_var}', data=data).fit()
    a = model2.params[x_var]
    a_pvalue = model2.pvalues[x_var]
    results['step2'] = {
        'model': 'M ~ X',
        'coefficient': a,
        'pvalue': a_pvalue,
        'significant': a_pvalue < 0.05
    }
    
    # 步骤3: M -> Y (路径b) 控制X
    model3 = ols(f'{y_var} ~ {x_var} + {m_var}', data=data).fit()
    b = model3.params[m_var]
    b_pvalue = model3.pvalues[m_var]
    c_prime = model3.params[x_var]  # 直接效应
    
    results['step3'] = {
        'model': 'Y ~ X + M',
        'coefficient_b': b,
        'pvalue_b': b_pvalue,
        'significant_b': b_pvalue < 0.05,
        'coefficient_c_prime': c_prime,
        'pvalue_c_prime': model3.pvalues[x_var]
    }
    
    # 计算间接效应和总效应
    indirect_effect = a * b
    total_effect = c
    proportion_mediated = indirect_effect / total_effect if total_effect != 0 else 0
    
    results['effects'] = {
        'total_effect': total_effect,
        'direct_effect': c_prime,
        'indirect_effect': indirect_effect,
        'proportion_mediated': proportion_mediated
    }
    
    # 判断中介效应是否存在
    results['mediation_present'] = (
        results['step1']['significant'] and 
        results['step2']['significant'] and 
        results['step3']['significant_b']
    )
    
    return results, model1, model2, model3

# 执行Baron & Kenny检验
bk_results, model1, model2, model3 = baron_kenny_test(
    mediation_data, 'education', 'skills', 'income'
)

print("Baron & Kenny逐步检验结果:")
print("=" * 50)

print("\n步骤1: X → Y (总效应)")
step1 = bk_results['step1']
print(f"模型: {step1['model']}")
print(f"系数 c: {step1['coefficient']:.4f}")
print(f"p值: {step1['pvalue']:.4f}")
print(f"是否显著: {'是' if step1['significant'] else '否'}")

print("\n步骤2: X → M (路径a)")
step2 = bk_results['step2']
print(f"模型: {step2['model']}")
print(f"系数 a: {step2['coefficient']:.4f}")
print(f"p值: {step2['pvalue']:.4f}")
print(f"是否显著: {'是' if step2['significant'] else '否'}")

print("\n步骤3: M → Y (路径b) 控制X")
step3 = bk_results['step3']
print(f"模型: {step3['model']}")
print(f"系数 b: {step3['coefficient_b']:.4f}")
print(f"p值: {step3['pvalue_b']:.4f}")
print(f"是否显著: {'是' if step3['significant_b'] else '否'}")
print(f"直接效应 c': {step3['coefficient_c_prime']:.4f}")

print("\n效应分解:")
effects = bk_results['effects']
print(f"总效应: {effects['total_effect']:.4f}")
print(f"直接效应: {effects['direct_effect']:.4f}")
print(f"间接效应: {effects['indirect_effect']:.4f}")
print(f"中介比例: {effects['proportion_mediated']:.4f}")

print(f"\n中介效应是否存在: {'是' if bk_results['mediation_present'] else '否'}")

# 可视化逐步检验结果
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# 步骤1: X -> Y
axes[0].scatter(mediation_data['education'], mediation_data['income'], alpha=0.5)
z = np.polyfit(mediation_data['education'], mediation_data['income'], 1)
p = np.poly1d(z)
axes[0].plot(mediation_data['education'], p(mediation_data['education']), "r-", linewidth=2)
axes[0].set_xlabel('教育年限 (X)')
axes[0].set_ylabel('收入 (Y)')
axes[0].set_title(f'步骤1: X → Y\nc = {step1["coefficient"]:.2f} (p = {step1["pvalue"]:.3f})')
axes[0].grid(True, alpha=0.3)

# 步骤2: X -> M
axes[1].scatter(mediation_data['education'], mediation_data['skills'], alpha=0.5)
z = np.polyfit(mediation_data['education'], mediation_data['skills'], 1)
p = np.poly1d(z)
axes[1].plot(mediation_data['education'], p(mediation_data['education']), "r-", linewidth=2)
axes[1].set_xlabel('教育年限 (X)')
axes[1].set_ylabel('工作技能 (M)')
axes[1].set_title(f'步骤2: X → M\na = {step2["coefficient"]:.2f} (p = {step2["pvalue"]:.3f})')
axes[1].grid(True, alpha=0.3)

# 步骤3: M -> Y (控制X)
axes[2].scatter(mediation_data['skills_resid'], mediation_data['income_resid'], alpha=0.5)
z = np.polyfit(mediation_data['skills_resid'], mediation_data['income_resid'], 1)
p = np.poly1d(z)
axes[2].plot(mediation_data['skills_resid'], p(mediation_data['skills_resid']), "r-", linewidth=2)
axes[2].set_xlabel('技能残差 (控制X)')
axes[2].set_ylabel('收入残差 (控制X)')
axes[2].set_title(f'步骤3: M → Y (控制X)\nb = {step3["coefficient_b"]:.2f} (p = {step3["pvalue_b"]:.3f})')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

系数乘积检验法（Sobel检验）

Sobel检验直接检验间接效应ab是否显著不为零，比逐步检验法有更高的统计功效。

# Sobel检验实现
def sobel_test(a, a_se, b, b_se):
    """
    执行Sobel检验
    """
    # 计算间接效应
    indirect_effect = a * b
    
    # 计算标准误 (Sobel公式)
    sobel_se = np.sqrt(a**2 * b_se**2 + b**2 * a_se**2)
    
    # 计算z统计量
    z_score = indirect_effect / sobel_se
    
    # 计算p值 (双侧检验)
    from scipy import stats
    p_value = 2 * (1 - stats.norm.cdf(np.abs(z_score)))
    
    return {
        'indirect_effect': indirect_effect,
        'sobel_se': sobel_se,
        'z_score': z_score,
        'p_value': p_value,
        'significant': p_value < 0.05
    }

# 从Baron-Kenny结果中获取参数
a = bk_results['step2']['coefficient']
a_se = model2.bse['education']
b = bk_results['step3']['coefficient_b']
b_se = model3.bse['skills']

# 执行Sobel检验
sobel_results = sobel_test(a, a_se, b, b_se)

print("Sobel检验结果:")
print("=" * 50)
print(f"间接效应 (ab): {sobel_results['indirect_effect']:.4f}")
print(f"Sobel标准误: {sobel_results['sobel_se']:.4f}")
print(f"Z统计量: {sobel_results['z_score']:.4f}")
print(f"p值: {sobel_results['p_value']:.4f}")
print(f"间接效应是否显著: {'是' if sobel_results['significant'] else '否'}")

# 比较两种方法
print(f"\n方法比较:")
print(f"逐步检验法结论: {'中介效应存在' if bk_results['mediation_present'] else '中介效应不存在'}")
print(f"Sobel检验结论: {'中介效应存在' if sobel_results['significant'] else '中介效应不存在'}")

Bootstrap法

Bootstrap法是目前推荐的中介效应检验方法，它不依赖正态分布假设，有更好的统计性质。

# Bootstrap中介效应检验
def bootstrap_mediation(data, x_var, m_var, y_var, n_bootstrap=5000, ci_level=0.95):
    """
    使用Bootstrap法检验中介效应
    """
    bootstrap_indirect = []
    bootstrap_direct = []
    bootstrap_total = []
    
    for i in range(n_bootstrap):
        # Bootstrap重抽样
        bootstrap_sample = resample(data, replace=True, random_state=i)
        
        try:
            # 估计路径a: X -> M
            model_a = ols(f'{m_var} ~ {x_var}', data=bootstrap_sample).fit()
            a_est = model_a.params[x_var]
            
            # 估计路径b和c': M -> Y 控制X
            model_b = ols(f'{y_var} ~ {x_var} + {m_var}', data=bootstrap_sample).fit()
            b_est = model_b.params[m_var]
            c_prime_est = model_b.params[x_var]
            
            # 估计总效应
            model_total = ols(f'{y_var} ~ {x_var}', data=bootstrap_sample).fit()
            c_est = model_total.params[x_var]
            
            # 计算间接效应
            indirect_est = a_est * b_est
            
            bootstrap_indirect.append(indirect_est)
            bootstrap_direct.append(c_prime_est)
            bootstrap_total.append(c_est)
            
        except:
            # 如果模型拟合失败，跳过这次迭代
            continue
    
    # 计算置信区间
    alpha = 1 - ci_level
    lower_idx = int(n_bootstrap * alpha / 2)
    upper_idx = int(n_bootstrap * (1 - alpha / 2))
    
    bootstrap_indirect_sorted = np.sort(bootstrap_indirect)
    bootstrap_direct_sorted = np.sort(bootstrap_direct)
    bootstrap_total_sorted = np.sort(bootstrap_total)
    
    results = {
        'indirect_effect': {
            'mean': np.mean(bootstrap_indirect),
            'ci_lower': bootstrap_indirect_sorted[lower_idx],
            'ci_upper': bootstrap_indirect_sorted[upper_idx],
            'significant': (bootstrap_indirect_sorted[lower_idx] > 0) or (bootstrap_indirect_sorted[upper_idx] < 0)
        },
        'direct_effect': {
            'mean': np.mean(bootstrap_direct),
            'ci_lower': bootstrap_direct_sorted[lower_idx],
            'ci_upper': bootstrap_direct_sorted[upper_idx]
        },
        'total_effect': {
            'mean': np.mean(bootstrap_total),
            'ci_lower': bootstrap_total_sorted[lower_idx],
            'ci_upper': bootstrap_total_sorted[upper_idx]
        }
    }
    
    return results, bootstrap_indirect

# 执行Bootstrap检验
print("执行Bootstrap中介效应检验...")
bootstrap_results, bootstrap_distribution = bootstrap_mediation(
    mediation_data, 'education', 'skills', 'income', n_bootstrap=1000  # 减少次数用于演示
)

print("Bootstrap检验结果:")
print("=" * 50)

print("\n间接效应 (ab):")
indirect = bootstrap_results['indirect_effect']
print(f"均值: {indirect['mean']:.4f}")
print(f"{95}%置信区间: [{indirect['ci_lower']:.4f}, {indirect['ci_upper']:.4f}]")
print(f"是否显著: {'是' if indirect['significant'] else '否'}")

print("\n直接效应 (c'):")
direct = bootstrap_results['direct_effect']
print(f"均值: {direct['mean']:.4f}")
print(f"{95}%置信区间: [{direct['ci_lower']:.4f}, {direct['ci_upper']:.4f}]")

print("\n总效应 (c):")
total = bootstrap_results['total_effect']
print(f"均值: {total['mean']:.4f}")
print(f"{95}%置信区间: [{total['ci_lower']:.4f}, {total['ci_upper']:.4f}]")

# 可视化Bootstrap分布
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# 间接效应分布
axes[0].hist(bootstrap_distribution, bins=30, alpha=0.7, color='lightblue', density=True)
axes[0].axvline(indirect['mean'], color='red', linestyle='-', label='均值')
axes[0].axvline(indirect['ci_lower'], color='red', linestyle='--', label='95% CI')
axes[0].axvline(indirect['ci_upper'], color='red', linestyle='--')
axes[0].axvline(0, color='black', linestyle='-', alpha=0.5)
axes[0].set_xlabel('间接效应 (ab)')
axes[0].set_ylabel('密度')
axes[0].set_title('Bootstrap间接效应分布')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# 三种效应比较
effects = ['总效应', '直接效应', '间接效应']
means = [total['mean'], direct['mean'], indirect['mean']]
ci_lowers = [total['ci_lower'], direct['ci_lower'], indirect['ci_lower']]
ci_uppers = [total['ci_upper'], direct['ci_upper'], indirect['ci_upper']]

axes[1].bar(effects, means, yerr=[np.array(means) - np.array(ci_lowers), 
                                 np.array(ci_uppers) - np.array(means)],
           capsize=5, alpha=0.7, color=['lightblue', 'lightgreen', 'lightcoral'])
axes[1].set_ylabel('效应大小')
axes[1].set_title('三种效应的Bootstrap估计')
axes[1].grid(True, alpha=0.3)

# 方法比较
methods = ['逐步检验', 'Sobel检验', 'Bootstrap']
conclusions = [
    '存在' if bk_results['mediation_present'] else '不存在',
    '存在' if sobel_results['significant'] else '不存在', 
    '存在' if indirect['significant'] else '不存在'
]
colors = ['green' if c == '存在' else 'red' for c in conclusions]

axes[2].barh(methods, [1, 1, 1], color=colors, alpha=0.7)
axes[2].set_xlim(0, 1)
axes[2].set_xlabel('结论')
axes[2].set_title('三种检验方法结论比较')
for i, (method, conclusion) in enumerate(zip(methods, conclusions)):
    axes[2].text(0.5, i, f'{conclusion}', va='center', ha='center', fontsize=12, color='white')

plt.tight_layout()
plt.show()

检验方法比较

不同中介效应检验方法各有优缺点：

IV. 实例分析：教育对收入的机制研究

研究背景与问题

我们继续深入探讨教育对收入影响的中介机制。假设我们怀疑教育通过三个主要途径影响收入：

1. 技能路径：教育提高工作技能
2. 网络路径：教育扩展社会网络
3. 信号路径：教育作为能力信号

数据生成与模型设定

# 生成更复杂的多中介数据
np.random.seed(2024)
n = 3000

# 自变量：教育年限
education = np.random.normal(12, 2, n)

# 中介变量1：工作技能（连续变量）
skills = 10 + 0.6 * education + np.random.normal(0, 1.5, n)

# 中介变量2：社会网络（连续变量）
network = 5 + 0.4 * education + 0.3 * skills + np.random.normal(0, 1.2, n)

# 中介变量3：教育信号（二元变量，是否名校）
prestige = np.random.binomial(1, 1/(1 + np.exp(-(0.5 * education - 6))))

# 结果变量：收入
income = (20 + 1.5 * education +  # 直接效应
          2.0 * skills +         # 技能路径
          1.2 * network +        # 网络路径  
          3.0 * prestige +       # 信号路径
          np.random.normal(0, 4, n))

# 创建数据框
multi_mediation_data = pd.DataFrame({
    'education': education,
    'skills': skills,
    'network': network,
    'prestige': prestige,
    'income': income
})

print("多中介数据描述:")
print(f"样本量: {n}")
print(f"教育年限: {education.mean():.2f} ± {education.std():.2f}")
print(f"工作技能: {skills.mean():.2f} ± {skills.std():.2f}")
print(f"社会网络: {network.mean():.2f} ± {network.std():.2f}")
print(f"名校比例: {prestige.mean():.3f}")
print(f"收入: {income.mean():.2f} ± {income.std():.2f}")

# 计算真实效应
true_effects = {
    'total': 1.5 + 0.6*2.0 + 0.4*1.2 + 0.5*3.0/(np.log(2)),  # 近似计算
    'direct': 1.5,
    'skills_indirect': 0.6 * 2.0,
    'network_indirect': 0.4 * 1.2,
    'prestige_indirect': 0.5 * 3.0 / (np.log(2))  # 逻辑回归系数近似
}

print(f"\n真实效应值:")
for effect, value in true_effects.items():
    print(f"{effect}: {value:.3f}")

并行多重中介分析

当存在多个中介变量时，我们需要使用并行多重中介模型。

# 并行多重中介分析
def parallel_mediation_analysis(data, x_var, m_vars, y_var, n_bootstrap=1000):
    """
    执行并行多重中介分析
    """
    results = {}
    bootstrap_distributions = {}
    
    # 总效应模型
    total_model = ols(f'{y_var} ~ {x_var}', data=data).fit()
    total_effect = total_model.params[x_var]
    
    # 直接效应模型（包含所有中介变量）
    mediators_formula = ' + '.join(m_vars)
    direct_model = ols(f'{y_var} ~ {x_var} + {mediators_formula}', data=data).fit()
    direct_effect = direct_model.params[x_var]
    
    # 存储各路径系数
    path_coefficients = {}
    
    # 路径a: X -> M
    for m_var in m_vars:
        model_a = ols(f'{m_var} ~ {x_var}', data=data).fit()
        path_coefficients[f'a_{m_var}'] = model_a.params[x_var]
        path_coefficients[f'a_se_{m_var}'] = model_a.bse[x_var]
    
    # 路径b: M -> Y (控制X和其他M)
    for m_var in m_vars:
        path_coefficients[f'b_{m_var}'] = direct_model.params[m_var]
        path_coefficients[f'b_se_{m_var}'] = direct_model.bse[m_var]
    
    # 计算间接效应
    indirect_effects = {}
    for m_var in m_vars:
        indirect = path_coefficients[f'a_{m_var}'] * path_coefficients[f'b_{m_var}']
        indirect_effects[m_var] = indirect
    
    total_indirect = sum(indirect_effects.values())
    
    # Bootstrap检验
    bootstrap_indirects = {m_var: [] for m_var in m_vars}
    bootstrap_directs = []
    bootstrap_totals = []
    
    for i in range(n_bootstrap):
        bootstrap_sample = resample(data, replace=True, random_state=i)
        
        try:
            # 总效应
            total_bs = ols(f'{y_var} ~ {x_var}', data=bootstrap_sample).fit()
            bootstrap_totals.append(total_bs.params[x_var])
            
            # 直接效应
            direct_bs = ols(f'{y_var} ~ {x_var} + {mediators_formula}', data=bootstrap_sample).fit()
            bootstrap_directs.append(direct_bs.params[x_var])
            
            # 各中介路径
            for m_var in m_vars:
                # 路径a
                a_bs = ols(f'{m_var} ~ {x_var}', data=bootstrap_sample).fit()
                a_coef = a_bs.params[x_var]
                
                # 路径b
                b_coef = direct_bs.params[m_var]
                
                bootstrap_indirects[m_var].append(a_coef * b_coef)
                
        except:
            continue
    
    # 计算Bootstrap置信区间
    ci_results = {}
    alpha = 0.05
    lower_idx = int(n_bootstrap * alpha / 2)
    upper_idx = int(n_bootstrap * (1 - alpha / 2))
    
    for m_var in m_vars:
        sorted_effects = np.sort(bootstrap_indirects[m_var])
        ci_results[m_var] = {
            'point_estimate': indirect_effects[m_var],
            'ci_lower': sorted_effects[lower_idx],
            'ci_upper': sorted_effects[upper_idx],
            'significant': (sorted_effects[lower_idx] > 0) or (sorted_effects[upper_idx] < 0)
        }
    
    # 总间接效应
    total_indirect_bs = []
    for i in range(min(len(bootstrap_indirects[m_vars[0]]), n_bootstrap)):
        total = sum(bootstrap_indirects[m_var][i] for m_var in m_vars)
        total_indirect_bs.append(total)
    
    sorted_total = np.sort(total_indirect_bs)
    ci_results['total_indirect'] = {
        'point_estimate': total_indirect,
        'ci_lower': sorted_total[lower_idx],
        'ci_upper': sorted_total[upper_idx],
        'significant': (sorted_total[lower_idx] > 0) or (sorted_total[upper_idx] < 0)
    }
    
    results = {
        'total_effect': total_effect,
        'direct_effect': direct_effect,
        'indirect_effects': indirect_effects,
        'path_coefficients': path_coefficients,
        'bootstrap_results': ci_results,
        'bootstrap_distributions': bootstrap_indirects
    }
    
    return results

# 执行并行多重中介分析
print("执行并行多重中介分析...")
mediators = ['skills', 'network', 'prestige']
multi_mediation_results = parallel_mediation_analysis(
    multi_mediation_data, 'education', mediators, 'income', n_bootstrap=500
)

print("多重中介分析结果:")
print("=" * 50)

print(f"\n总效应: {multi_mediation_results['total_effect']:.4f}")
print(f"直接效应: {multi_mediation_results['direct_effect']:.4f}")

print(f"\n各路径间接效应:")
for mediator in mediators:
    indirect = multi_mediation_results['indirect_effects'][mediator]
    print(f"{mediator}: {indirect:.4f}")

total_indirect = sum(multi_mediation_results['indirect_effects'].values())
print(f"总间接效应: {total_indirect:.4f}")

print(f"\nBootstrap检验结果 (95%置信区间):")
for mediator in mediators + ['total_indirect']:
    bs_result = multi_mediation_results['bootstrap_results'][mediator]
    sig_flag = "***" if bs_result['significant'] else ""
    print(f"{mediator}: {bs_result['point_estimate']:.4f} [{bs_result['ci_lower']:.4f}, {bs_result['ci_upper']:.4f}] {sig_flag}")

# 可视化多重中介结果
fig, axes = plt.subplots(2, 2, figsize=(15, 10))

# 1. 效应分解饼图
effects = ['直接效应', '技能中介', '网络中个', '信号中介']
effect_values = [
    multi_mediation_results['direct_effect'],
    multi_mediation_results['indirect_effects']['skills'],
    multi_mediation_results['indirect_effects']['network'], 
    multi_mediation_results['indirect_effects']['prestige']
]

axes[0,0].pie(effect_values, labels=effects, autopct='%1.1f%%', startangle=90)
axes[0,0].set_title('教育对收入影响的效应分解')

# 2. 各路径间接效应比较
mediator_names = ['工作技能', '社会网络', '教育信号']
indirect_values = [multi_mediation_results['indirect_effects'][m] for m in mediators]
ci_lowers = [multi_mediation_results['bootstrap_results'][m]['ci_lower'] for m in mediators]
ci_uppers = [multi_mediation_results['bootstrap_results'][m]['ci_upper'] for m in mediators]

y_pos = np.arange(len(mediator_names))
axes[0,1].barh(y_pos, indirect_values, xerr=[np.array(indirect_values) - np.array(ci_lowers), 
                                           np.array(ci_uppers) - np.array(indirect_values)],
              capsize=5, alpha=0.7, color=['lightblue', 'lightgreen', 'lightcoral'])
axes[0,1].set_yticks(y_pos)
axes[0,1].set_yticklabels(mediator_names)
axes[0,1].set_xlabel('间接效应大小')
axes[0,1].set_title('各中介路径的间接效应')
axes[0,1].axvline(x=0, color='black', linestyle='-', alpha=0.3)

# 3. Bootstrap分布
for i, mediator in enumerate(mediators):
    axes[1,0].hist(multi_mediation_results['bootstrap_distributions'][mediator], 
                   bins=20, alpha=0.6, label=mediator_names[i], density=True)

axes[1,0].axvline(x=0, color='black', linestyle='-', alpha=0.5)
axes[1,0].set_xlabel('间接效应')
axes[1,0].set_ylabel('密度')
axes[1,0].set_title('各路径间接效应的Bootstrap分布')
axes[1,0].legend()

# 4. 路径系数可视化
path_coeffs = multi_mediation_results['path_coefficients']
path_labels = ['X→M1\n(技能)', 'X→M2\n(网络)', 'X→M3\n(信号)', 
               'M1→Y\n(技能)', 'M2→Y\n(网络)', 'M3→Y\n(信号)', 'X→Y\n(直接)']
path_values = [
    path_coeffs['a_skills'], path_coeffs['a_network'], path_coeffs.get('a_prestige', 0),
    path_coeffs['b_skills'], path_coeffs['b_network'], path_coeffs['b_prestige'],
    multi_mediation_results['direct_effect']
]

colors = ['blue']*3 + ['green']*3 + ['red']
axes[1,1].bar(path_labels, path_values, color=colors, alpha=0.7)
axes[1,1].set_ylabel('系数大小')
axes[1,1].set_title('各路径系数估计')
axes[1,1].tick_params(axis='x', rotation=45)

plt.tight_layout()
plt.show()

结果解释与理论意义

通过多重中介分析，我们获得了以下重要发现：

1. 技能路径是最重要的机制：教育主要通过提高工作技能来增加收入
2. 网络效应的稳健性：社会网络扩展也是显著的中介路径
3. 信号效应的有限性：教育信号的作用相对较小且不确定性较高
4. 直接效应的持续性：即使在控制所有中介后，教育仍有直接效应

这些发现对教育政策有重要启示：应该注重教育质量的提升（技能路径）而不仅仅是教育年限的延长。

V. 复杂中介模型扩展

链式中介模型（序列中介）

当中介变量之间存在先后顺序时，我们需要使用链式中介模型。

# 链式中介模型分析
def serial_mediation_analysis(data, x_var, m_sequence, y_var, n_bootstrap=1000):
    """
    执行链式中介分析
    m_sequence: 中介变量序列，如 ['M1', 'M2', 'M3']
    """
    
    results = {}
    
    # 路径系数估计
    path_coefficients = {}
    
    # 路径1: X -> M1
    model1 = ols(f'{m_sequence[0]} ~ {x_var}', data=data).fit()
    path_coefficients[f'a1'] = model1.params[x_var]
    
    # 路径2: X -> M2, M1 -> M2
    model2 = ols(f'{m_sequence[1]} ~ {x_var} + {m_sequence[0]}', data=data).fit()
    path_coefficients[f'a2'] = model2.params[x_var]  # X -> M2 (直接)
    path_coefficients[f'd21'] = model2.params[m_sequence[0]]  # M1 -> M2
    
    # 路径3: X -> M3, M1 -> M3, M2 -> M3 (如果有更多中介)
    if len(m_sequence) > 2:
        model3 = ols(f'{m_sequence[2]} ~ {x_var} + {m_sequence[0]} + {m_sequence[1]}', data=data).fit()
        path_coefficients[f'a3'] = model3.params[x_var]
        path_coefficients[f'd31'] = model3.params[m_sequence[0]]
        path_coefficients[f'd32'] = model3.params[m_sequence[1]]
    
    # 结果模型: Y ~ X + 所有M
    mediators_formula = ' + '.join(m_sequence)
    outcome_model = ols(f'{y_var} ~ {x_var} + {mediators_formula}', data=data).fit()
    
    # 存储b路径系数
    for i, m_var in enumerate(m_sequence):
        path_coefficients[f'b{i+1}'] = outcome_model.params[m_var]
    
    path_coefficients['c_prime'] = outcome_model.params[x_var]
    
    # 计算各种间接效应
    indirect_effects = {}
    
    # 特定路径间接效应
    if len(m_sequence) == 2:
        # X -> M1 -> Y
        indirect_effects['X_M1_Y'] = path_coefficients['a1'] * path_coefficients['b1']
        # X -> M2 -> Y  
        indirect_effects['X_M2_Y'] = path_coefficients['a2'] * path_coefficients['b2']
        # X -> M1 -> M2 -> Y
        indirect_effects['X_M1_M2_Y'] = path_coefficients['a1'] * path_coefficients['d21'] * path_coefficients['b2']
        
    elif len(m_sequence) == 3:
        # 更复杂的路径...
        pass
    
    total_indirect = sum(indirect_effects.values())
    
    # Bootstrap检验 (简化版)
    bootstrap_results = {}
    for effect_name in indirect_effects.keys():
        bootstrap_results[effect_name] = []
    
    for i in range(n_bootstrap):
        bootstrap_sample = resample(data, replace=True, random_state=i)
        
        try:
            # 重新估计所有模型
            # 这里简化实现，实际需要完整重新估计
            pass
        except:
            continue
    
    results = {
        'path_coefficients': path_coefficients,
        'indirect_effects': indirect_effects,
        'total_indirect': total_indirect,
        'direct_effect': path_coefficients['c_prime']
    }
    
    return results

# 模拟链式中介数据
np.random.seed(123)
n_serial = 2000

education_serial = np.random.normal(12, 2, n_serial)

# 链式: 教育 -> 认知能力 -> 工作技能 -> 收入
cognition = 100 + 2.0 * education_serial + np.random.normal(0, 5, n_serial)
skills_serial = 10 + 0.3 * education_serial + 0.5 * cognition + np.random.normal(0, 2, n_serial)
income_serial = (20 + 1.0 * education_serial + 0.1 * cognition + 2.5 * skills_serial + 
                np.random.normal(0, 3, n_serial))

serial_data = pd.DataFrame({
    'education': education_serial,
    'cognition': cognition,
    'skills': skills_serial,
    'income': income_serial
})

print("链式中介分析示例:")
serial_sequence = ['cognition', 'skills']
serial_results = serial_mediation_analysis(serial_data, 'education', serial_sequence, 'income')

print("路径系数:")
for path, coef in serial_results['path_coefficients'].items():
    print(f"{path}: {coef:.4f}")

print("\n间接效应分解:")
for effect, value in serial_results['indirect_effects'].items():
    print(f"{effect}: {value:.4f}")

print(f"总间接效应: {serial_results['total_indirect']:.4f}")
print(f"直接效应: {serial_results['direct_effect']:.4f}")

调节中介模型

当中介效应在不同群体或条件下不同时，我们需要考虑调节中介模型。

# 调节中介模型分析
def moderated_mediation_analysis(data, x_var, m_var, y_var, moderator_var, n_bootstrap=1000):
    """
    执行调节中介分析
    """
    results = {}
    
    # 1. 检验调节效应在路径a (X -> M)
    model_a = ols(f'{m_var} ~ {x_var} * {moderator_var}', data=data).fit()
    
    # 2. 检验调节效应在路径b (M -> Y)
    model_b = ols(f'{y_var} ~ {x_var} + {m_var} * {moderator_var}', data=data).fit()
    
    # 提取系数
    a1 = model_a.params[x_var]  # X -> M 的主效应
    a3 = model_a.params[f'{x_var}:{moderator_var}']  # 调节效应在路径a
    
    b1 = model_b.params[m_var]  # M -> Y 的主效应  
    b3 = model_b.params[f'{m_var}:{moderator_var}']  # 调节效应在路径b
    
    # 条件间接效应
    # 通常在调节变量的不同水平上计算
    moderator_values = np.percentile(data[moderator_var], [10, 50, 90])
    conditional_indirect_effects = {}
    
    for i, mod_val in enumerate(moderator_values):
        # 条件路径a
        conditional_a = a1 + a3 * mod_val
        # 条件路径b
        conditional_b = b1 + b3 * mod_val
        # 条件间接效应
        conditional_indirect = conditional_a * conditional_b
        
        conditional_indirect_effects[f'percentile_{[10,50,90][i]}'] = {
            'moderator_value': mod_val,
            'indirect_effect': conditional_indirect,
            'path_a': conditional_a,
            'path_b': conditional_b
        }
    
    results = {
        'model_a': model_a,
        'model_b': model_b,
        'conditional_indirect_effects': conditional_indirect_effects,
        'moderator_values': moderator_values
    }
    
    return results

# 模拟调节中介数据
np.random.seed(456)
n_mod = 2000

education_mod = np.random.normal(12, 2, n_mod)
motivation = np.random.normal(0, 1, n_mod)  # 调节变量: 学习动机

# 调节路径a: 教育对技能的影响受动机调节
skills_mod = (10 + 0.5 * education_mod + 0.3 * motivation + 
              0.2 * education_mod * motivation + np.random.normal(0, 1, n_mod))

# 调节路径b: 技能对收入的影响受动机调节
income_mod = (20 + 1.0 * education_mod + 2.0 * skills_mod + 0.5 * motivation +
              0.3 * skills_mod * motivation + np.random.normal(0, 3, n_mod))

moderated_data = pd.DataFrame({
    'education': education_mod,
    'motivation': motivation,
    'skills': skills_mod,
    'income': income_mod
})

print("调节中介分析示例:")
moderated_results = moderated_mediation_analysis(
    moderated_data, 'education', 'skills', 'income', 'motivation'
)

print("条件间接效应:")
for percentile, effects in moderated_results['conditional_indirect_effects'].items():
    print(f"{percentile}:")
    print(f"  调节变量值: {effects['moderator_value']:.3f}")
    print(f"  路径a: {effects['path_a']:.3f}")
    print(f"  路径b: {effects['path_b']:.3f}")
    print(f"  间接效应: {effects['indirect_effect']:.3f}")

# 可视化调节中介效应
mod_values = np.linspace(moderated_data['motivation'].min(), 
                        moderated_data['motivation'].max(), 100)

conditional_indirect = []
for mod_val in mod_values:
    a_cond = (moderated_results['model_a'].params['education'] + 
              moderated_results['model_a'].params['education:motivation'] * mod_val)
    b_cond = (moderated_results['model_b'].params['skills'] + 
              moderated_results['model_b'].params['skills:motivation'] * mod_val)
    conditional_indirect.append(a_cond * b_cond)

plt.figure(figsize=(10, 6))
plt.plot(mod_values, conditional_indirect, 'b-', linewidth=2)
plt.xlabel('学习动机 (调节变量)')
plt.ylabel('条件间接效应')
plt.title('调节中介效应: 间接效应随调节变量变化')
plt.grid(True, alpha=0.3)

# 标记特定百分位点
percentiles = [10, 50, 90]
for p in percentiles:
    mod_val = np.percentile(moderated_data['motivation'], p)
    idx = np.argmin(np.abs(mod_values - mod_val))
    plt.plot(mod_val, conditional_indirect[idx], 'ro', markersize=8)
    plt.text(mod_val, conditional_indirect[idx] + 0.1, f'{p}%', 
             ha='center', va='bottom')

plt.show()

复杂模型选择指南

根据研究问题的复杂性，可以选择不同的中介模型：

VI. 中介效应分析的注意事项

常见误区与解决方法

中介效应分析在实践中存在多个常见误区，需要特别注意：

# 中介分析误区诊断
def diagnose_mediation_issues(data, x_var, m_vars, y_var, covariates=None):
    """
    诊断中介分析中的常见问题
    """
    issues = []
    recommendations = []
    
    # 1. 样本量检查
    n = len(data)
    if n < 100:
        issues.append("样本量可能不足")
        recommendations.append("考虑增加样本量或使用Bootstrap方法")
    
    # 2. 多重共线性检查
    if covariates is not None:
        all_vars = [x_var] + m_vars + covariates + [y_var]
        correlation_matrix = data[all_vars].corr().abs()
        high_corr = (correlation_matrix > 0.8) & (correlation_matrix < 1.0)
        if high_corr.any().any():
            issues.append("存在高度相关的变量")
            recommendations.append("检查多重共线性，考虑变量选择或正则化")
    
    # 3. 测量误差敏感性
    # 模拟测量误差的影响（简化演示）
    print("进行测量误差敏感性分析...")
    
    # 4. 模型设定检查
    for m_var in m_vars:
        # 检查X->M关系
        model_a = ols(f'{m_var} ~ {x_var}', data=data).fit()
        if model_a.rsquared < 0.01:
            issues.append(f"X对{M}的解释力很低")
            recommendations.append(f"检查{M}作为中介变量的合理性")
    
    return issues, recommendations

# 执行诊断
issues, recommendations = diagnose_mediation_issues(
    multi_mediation_data, 'education', ['skills', 'network', 'prestige'], 'income'
)

print("中介分析问题诊断:")
if issues:
    print("发现的问题:")
    for i, issue in enumerate(issues, 1):
        print(f"  {i}. {issue}")
    
    print("\n改进建议:")
    for i, rec in enumerate(recommendations, 1):
        print(f"  {i}. {rec}")
else:
    print("未发现明显问题")

# 敏感性分析函数
def sensitivity_analysis(data, x_var, m_var, y_var, confounder_strength_range):
    """
    中介效应估计对未观测混淆的敏感性分析
    """
    sensitivity_results = []
    
    for rho in confounder_strength_range:
        # 模拟未观测混淆的影响
        # 这里使用简化的敏感性分析公式
        original_indirect = multi_mediation_results['indirect_effects'][m_var]
        
        # 敏感性调整（简化模型）
        # 实际敏感性分析更复杂，需要更多假设
        adjusted_indirect = original_indirect * (1 - rho**2)
        
        sensitivity_results.append({
            'confounder_strength': rho,
            'original_indirect': original_indirect,
            'adjusted_indirect': adjusted_indirect,
            'bias': original_indirect - adjusted_indirect,
            'bias_percentage': (original_indirect - adjusted_indirect) / original_indirect * 100
        })
    
    return pd.DataFrame(sensitivity_results)

# 执行敏感性分析
confounder_strengths = [0.1, 0.2, 0.3, 0.4, 0.5]
sensitivity_df = sensitivity_analysis(
    multi_mediation_data, 'education', 'skills', 'income', confounder_strengths
)

print("\n敏感性分析结果:")
print(sensitivity_df.round(4))

# 可视化敏感性分析
plt.figure(figsize=(10, 6))
plt.plot(sensitivity_df['confounder_strength'], 
         sensitivity_df['adjusted_indirect'], 'bo-', linewidth=2, markersize=8)
plt.axhline(y=sensitivity_df['original_indirect'].iloc[0], color='red', 
           linestyle='--', label='原始估计')
plt.xlabel('未观测混淆的强度 (ρ)')
plt.ylabel('调整后的间接效应')
plt.title('中介效应估计对未观测混淆的敏感性')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

最佳实践总结

基于理论研究和实践经验，我们总结出中介效应分析的最佳实践：