R语言mediation包实战：用移民数据手把手教你做中介效应分析（附完整代码）

R语言mediation包实战：从移民数据到中介效应分析的完整指南

在社会科学和医学研究中，我们常常需要探究变量间的复杂关系——不仅仅是"X是否影响Y"，更要理解"X如何通过M影响Y"。这种机制分析正是中介效应模型的核心价值所在。R语言的mediation包由统计学家Kosuke Imai团队开发，为这一分析需求提供了强大而灵活的工具。不同于简单的统计检验，中介分析能帮助我们拆解效应路径，区分直接效应和间接效应，从而更深入地理解变量间的因果关系网络。

本文将基于经典的移民焦虑研究案例，手把手带你完成从数据准备到结果解读的全流程。我们不仅会详细讲解代码实现，更会剖析每一步背后的统计逻辑，帮助你避开常见陷阱，并将这套方法迁移到自己的研究课题中。无论你是心理学研究生、公共卫生研究者还是市场分析师，掌握中介效应分析都将显著提升你的研究深度。

1. 数据准备与变量梳理

1.1 理解研究问题与数据结构

我们使用的数据集来自一项移民行为研究，包含以下关键变量：

• treat：二分变量（0=未接受治疗，1=接受治疗）
• emo：连续变量，通过量表测量的焦虑程度得分
• cong_mesg：二分结果变量（0=未发送移民申请，1=发送申请）

研究假设：心理治疗（treat）通过降低焦虑水平（emo）进而减少移民倾向（cong_mesg）。这里emo就是假设的中介变量。

首先加载必要的包并导入数据：

library(mediation) library(tidyverse) # 假设数据已保存为CSV文件 migration_data <- read_csv("migration_study.csv") # 快速查看数据结构 glimpse(migration_data)

1.2 变量预处理的关键检查点

在建模前必须进行以下诊断检查：

1. 缺失值处理：

   colSums(is.na(migration_data))

   若存在缺失，需决定是删除还是插补。mediation包对缺失值敏感，建议使用mice包进行多重插补。

2. 变量类型验证：

   • 确认treat和cong_mesg是因子类型
   • 检查emo的分布是否近似正态（对线性模型重要）

   migration_data <- migration_data %>% mutate(treat = as.factor(treat), cong_mesg = as.factor(cong_mesg)) ggplot(migration_data, aes(x=emo)) + geom_histogram(bins=30, fill="steelblue")

3. 协变量平衡检查（尤其对实验数据重要）：

   table(migration_data$treat, migration_data$gender)

注意：若中介变量是二分变量，需要修改模型设定（使用glm而非lm）。本文以连续中介变量为例。

2. 模型构建：从理论到统计实现

2.1 建立中介关系方程

中介分析需要构建两个模型：

1. 中介变量模型（M ~ X + C）
2. 结果变量模型（Y ~ X + M + C）

对于我们的案例：

# 模型1：治疗对焦虑的影响（线性回归） med_model <- lm(emo ~ treat + age + educ + gender + income, data = migration_data) # 模型2：治疗和焦虑对移民决定的影响（probit回归） out_model <- glm(cong_mesg ~ emo + treat + age + educ + gender + income, data = migration_data, family = binomial("probit"))

2.2 关键建模决策点

1. 连接函数选择：

   • 连续结果变量：lm()
   • 二分结果变量：glm() with logit/probit
   • 计数数据：glm() with poisson/negative binomial

2. 协变量控制： 应包含所有理论上相关的混淆变量。遗漏重要协变量会导致估计偏差。

3. 交互项考量： 若怀疑X对Y的影响随M变化（调节效应），需在结果模型中添加X*M交互项。

2.3 模型诊断与修正

在运行中介分析前，必须验证模型假设：

# 检查线性模型假设 par(mfrow=c(2,2)) plot(med_model) # 检查probit模型拟合 library(performance) check_model(out_model)

若发现异方差性等问题，可考虑：

• 对连续变量进行变换（如log）
• 使用稳健标准误
• 尝试其他连接函数

3. 中介效应估计与解读

3.1 运行mediation分析

set.seed(123) # 确保结果可重复 med_result <- mediate( model.m = med_model, model.y = out_model, treat = "treat", # 自变量名称 mediator = "emo", # 中介变量名称 robustSE = TRUE, # 使用稳健标准误 sims = 1000 # 建议至少1000次模拟 )

参数说明：

• sims：Bootstrap重抽样次数，影响结果精度
• boot.ci.type：置信区间类型（"perc"或"bca"）
• control.value和treat.value：可自定义比较组

3.2 结果解读框架

查看汇总报告：

summary(med_result)

输出包含三个关键效应量：

可视化呈现：

plot(med_result)

3.3 效应量计算与比例分析

计算中介效应占比：

prop_mediated <- med_result$d1/med_result$tau print(paste("中介效应占比:", round(prop_mediated*100, 1), "%"))

提示：当总效应不显著而中介效应显著时，可能存在"不一致中介"现象，需要谨慎解释。

4. 高级分析与敏感性检验

4.1 敏感性分析：隐藏混淆的影响

中介分析的关键假设是"无未测量混淆"。我们可以检验这一假设的稳健性：

sens_result <- medsens(med_result, rho.by = 0.1, effect.type = "indirect", sims = 500) summary(sens_result) plot(sens_result, sens.par = "rho")

解读要点：

• Rho=0表示无隐藏混淆时的估计
• 红线显示使效应不再显著所需的混淆强度
• 绝对值越大，结论越稳健

4.2 多类别处理与多中介分析

当处理变量有多个类别时（如不同剂量组）：

# 对比治疗组1 vs 对照组 med_result1 <- mediate(..., treat.value = "0", control.value = "1") # 对比治疗组2 vs 对照组 med_result2 <- mediate(..., treat.value = "0", control.value = "2")

对于并行多重中介：

# 分别建立各中介模型 med1_model <- lm(emo1 ~ treat + covariates, data) med2_model <- lm(emo2 ~ treat + covariates, data) # 结果模型包含所有中介 out_multi_model <- glm(outcome ~ treat + emo1 + emo2 + covariates, family = binomial("probit")) # 使用multimed函数 multi_result <- multimed(..., mediator = c("emo1","emo2"))

4.3 纵向数据与交叉滞后模型

对于时间序列数据，可考虑：

library(lavaan) clpm <- ' # 自回归路径 emo_t2 ~ a1*emo_t1 + b1*treat_t1 cong_t2 ~ c1*cong_t1 + d1*emo_t1 + d2*treat_t1 # 交叉滞后路径 emo_t2 ~ e1*cong_t1 cong_t2 ~ f1*emo_t1 # 间接效应 indirect := b1 * d1 ' fit <- sem(clpm, data = panel_data)

5. 迁移应用到其他研究场景

5.1 医学研究案例：药物疗效机制

假设研究抗抑郁药（X）通过改善睡眠质量（M）影响患者生活质量（Y）：

# 数据准备 clinic_data <- read_csv("depression_trial.csv") %>% mutate(drug = factor(drug, levels = c("placebo","SSRI"))) # 模型构建 med_model <- lm(sleep_quality ~ drug + age + severity, data = clinic_data) out_model <- glm(qol ~ sleep_quality + drug + age + severity, family = gaussian, data = clinic_data) # 中介分析 med_result <- mediate(med_model, out_model, treat = "drug", mediator = "sleep_quality", sims = 1000)

5.2 市场营销应用：广告投放效果

分析广告支出（X）通过品牌认知度（M）影响销售额（Y）：

# 考虑非线性关系 med_model <- lm(brand_awareness ~ log(ad_spend) + market_size, data = marketing) out_model <- glm(sales ~ brand_awareness + log(ad_spend) + market_size, family = poisson, data = marketing) # 使用BCa置信区间 med_result <- mediate(med_model, out_model, treat = "ad_spend", mediator = "brand_awareness", boot.ci.type = "bca")

5.3 教育研究示例

探究教学方法（X）通过学生参与度（M）影响学业成绩（Y）：

# 多水平模型考虑班级嵌套 library(lme4) med_model <- lmer(engagement ~ method + pretest + (1|class), data = edu_data) out_model <- glmer(score ~ engagement + method + pretest + (1|class), family = Gamma(link="log"), data = edu_data) # 使用贝叶斯方法处理复杂模型 library(mediation) med_result <- mediate(med_model, out_model, treat = "method", mediator = "engagement", sims = 500, boot = TRUE)

6. 常见问题与解决方案

6.1 模型收敛问题

问题表现：

• 警告信息："maximum number of iterations reached"
• 估计值异常大/小

解决方案：

1. 增加迭代次数：

   control = glm.control(maxit = 100)

2. 重新缩放连续预测变量：

   data <- data %>% mutate(age_s = scale(age))

3. 尝试不同的优化算法

6.2 结果不稳定

问题表现：

• 每次运行结果差异大
• 置信区间异常宽

解决方案：

1. 增加sims参数（至少1000次）

   med_result <- mediate(..., sims = 5000)

2. 检查多重共线性：

   car::vif(med_model)

3. 使用更稳健的置信区间方法：

   boot.ci.type = "bca"

6.3 效应量解释困难

实用解读框架：

1. 标准化效应：

   # 计算标准化中介效应 std_effect <- med_result$d1/sd(data$cong_mesg, na.rm=TRUE)

2. 临床/实践显著性：

   • 结合领域知识判断效应大小
   • 计算需治疗人数(NNT)等指标

3. 对比其他研究：

   • 在文献背景下解释结果
   • 进行元分析比较

6.4 代码优化技巧

1. 并行计算加速：

   library(parallel) cl <- makeCluster(4) med_result <- mediate(..., parallel = TRUE, cl = cl) stopCluster(cl)

2. 结果自动化报告：

   library(rmarkdown) render("mediation_report.Rmd", params = list(result = med_result))

3. 创建可复用函数：

   run_mediation <- function(data, treat, mediator, outcome) { # 封装完整分析流程 ... return(list(result = med_result, plots = med_plots)) }

7. 完整代码模板与扩展资源

7.1 可复用的代码模板

# 中介分析通用模板 mediation_analysis <- function(data, treat_var, mediator_var, outcome_var, covariates, outcome_type = "continuous", sims = 1000) { # 预处理数据 data <- data %>% mutate(across(all_of(c(treat_var, outcome_var)), as.factor)) # 构建中介模型 med_formula <- paste(mediator_var, "~", treat_var, paste(covariates, collapse = " + "), sep = " ") med_model <- if(outcome_type == "continuous") { lm(as.formula(med_formula), data = data) } else { glm(as.formula(med_formula), family = binomial, data = data) } # 构建结果模型 out_formula <- paste(outcome_var, "~", mediator_var, "+", treat_var, paste(covariates, collapse = " + "), sep = " ") out_model <- if(outcome_type == "continuous") { lm(as.formula(out_formula), data = data) } else { glm(as.formula(out_formula), family = binomial, data = data) } # 运行中介分析 set.seed(123) mediate(med_model, out_model, treat = treat_var, mediator = mediator_var, sims = sims, robustSE = TRUE) } # 示例使用 result <- mediation_analysis( data = migration_data, treat_var = "treat", mediator_var = "emo", outcome_var = "cong_mesg", covariates = c("age","educ","gender","income"), outcome_type = "binary" )

7.2 扩展学习资源

1. 官方文档：

   • ?mediation::mediate
   • mediation包vignette

2. 可视化进阶：

   library(ggplot2) library(broom) # 绘制路径系数图 path_plot <- function(med_result) { tidy_m <- tidy(med_result$model.m) tidy_y <- tidy(med_result$model.y) # 构建绘图数据 # ...具体绘图代码 }

3. 替代方案比较：

   • lavaan：结构方程模型框架
   • brms：贝叶斯中介分析
   • mediation：本文介绍的经典方法

7.3 研究设计建议

1. 样本量规划：

   • 使用power4mediation包进行功效分析
   • 一般规则：每个参数至少20-30个样本

2. 测量时点设计：

   • 中介变量测量时机至关重要
   • 考虑交叉滞后面板设计

3. 预注册分析计划：

   • 预先明确假设和模型设定
   • 减少研究者自由度问题

在实际分析中，我发现最常被忽视的步骤是模型诊断和敏感性分析。许多研究者直接跳转到结果解读，而忽略了检查模型假设是否满足。特别是在处理非连续变量时，连接函数的选择会显著影响结论。另一个实用技巧是将sims参数设置为至少1000次——虽然会增加计算时间，但能获得更稳定的结果。