中介效应分析:探究因果机制的有效方法
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方法后变得更加完善。
三类主要效应定义
在中介效应分析中,我们主要关注三类效应:
效应类型 | 符号表示 | 定义解释 | 数学表达式 |
---|---|---|---|
总效应 | c | X对Y的总影响 | Y = i₁ + cX + e₁ |
直接效应 | c’ | X对Y的直接影响 | Y = i₂ + c’X + bM + e₂ |
间接效应 | ab | X通过M影响Y的效应 | M = i₃ + aX + e₃ |
中介效应的识别假设
要获得无偏的中介效应估计,需要满足以下关键假设:
- 无混淆假设:没有未观测的混淆变量影响X→Y、X→M、M→Y关系
- 一致性假设:观测到的处理水平就是实际接受的处理
- 模型设定正确:中介模型的形式设定正确
- 无交互作用:处理变量和中介变量没有交互作用(可放松)
# 模拟基本的中介效应数据
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()
中介模型的理论基础
中介效应分析建立在因果推断的潜在结果框架上,近年来发展的因果中介分析框架提供了更严谨的理论基础。
Lexical error on line 10. Unrecognized text. ... D[直接效应 c'] -.->|X → Y| C E[间接效应 ab -----------------------^III. 中介效应的检验方法
逐步检验法(Baron & Kenny方法)
这是最传统的中介效应检验方法,包含四个步骤:
- 检验X对Y的总效应(路径c)
- 检验X对M的效应(路径a)
- 检验M对Y的效应(路径b)
- 检验加入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()
检验方法比较
不同中介效应检验方法各有优缺点:
检验方法 | 原理 | 优点 | 缺点 | 适用场景 |
---|---|---|---|---|
逐步检验法 | 分步检验各路径显著性 | 直观易懂,解释性强 | 统计功效较低,Type I错误率控制不佳 | 初步探索,样本量充足 |
Sobel检验 | 直接检验ab乘积 | 比逐步法功效高 | 依赖正态分布假设 | 中等样本量,正态性满足 |
Bootstrap法 | 重抽样获得经验分布 | 不依赖分布假设,准确性高 | 计算成本高 | 推荐方法,各种样本量 |
IV. 实例分析:教育对收入的机制研究
研究背景与问题
我们继续深入探讨教育对收入影响的中介机制。假设我们怀疑教育通过三个主要途径影响收入:
- 技能路径:教育提高工作技能
- 网络路径:教育扩展社会网络
- 信号路径:教育作为能力信号
数据生成与模型设定
# 生成更复杂的多中介数据
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()
结果解释与理论意义
通过多重中介分析,我们获得了以下重要发现:
- 技能路径是最重要的机制:教育主要通过提高工作技能来增加收入
- 网络效应的稳健性:社会网络扩展也是显著的中介路径
- 信号效应的有限性:教育信号的作用相对较小且不确定性较高
- 直接效应的持续性:即使在控制所有中介后,教育仍有直接效应
这些发现对教育政策有重要启示:应该注重教育质量的提升(技能路径)而不仅仅是教育年限的延长。
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()
最佳实践总结
基于理论研究和实践经验,我们总结出中介效应分析的最佳实践:
分析阶段 | 最佳实践 | 避免的误区 |
---|---|---|
研究设计 | 理论驱动的中介选择,预注册分析计划 | 数据驱动的中介探索,p-hacking |
模型设定 | 基于理论选择模型形式,考虑时序关系 | 忽略变量间因果顺序,错误设定 |
估计方法 | 使用Bootstrap法,报告置信区间 | 仅依赖逐步法,不报告不确定性 |
结果解释 | 谨慎解释因果机制,考虑竞争性解释 | 过度解读相关关系为因果关系 |
稳健性检验 | 进行敏感性分析,检验不同设定 | 忽略模型假设和局限性 |
- 点赞
- 收藏
- 关注作者
评论(0)