R语言mediation包进阶：GLMM处理分类变量中介效应的完整解决方案

1. 社会科学研究中的分类变量中介分析挑战

在心理学、教育学等社会科学研究中，我们经常遇到自变量、中介变量和因变量均为分类变量的情况。例如，研究学校类型（天主教/非天主教）是否通过学生对学校的依恋程度（喜欢/不喜欢）影响校园暴力行为（发生/未发生）。这类研究设计需要特殊的中介效应分析方法。

传统的中介效应分析通常假设连续型变量和线性关系，但社会科学数据往往具有以下特征：

• 嵌套结构：学生嵌套在班级/学校中
• 分类变量：二分或多分类的测量尺度
• 非正态分布：特别是当因变量为罕见事件时

GLMM（广义线性混合模型）与mediation包的组合为解决这些问题提供了理想方案。这种组合能够：

• 处理分类变量
• 考虑组内相关性
• 提供准确的效应量估计

2. 数据准备与变量编码

2.1 示例数据集结构

我们使用一个虚构的学生数据集，包含以下关键变量：

变量名	类型	描述	编码
SCH_ID	分类	学校ID	数值
catholic	二分	是否天主教学校	0=非天主教,1=天主教
attachment	二分	学校依恋程度	0=不喜欢,1=喜欢
fight	二分	是否发生打架	0=否,1=是
gender	二分	学生性别	0=男,1=女
income	有序	家庭收入等级	1-13级

# 数据导入与查看
library(tidyverse)
student <- read_csv("student_data.csv") %>% 
  mutate(across(c(catholic, attachment, fight, gender), as.factor))
glimpse(student)

2.2 分类变量编码原则

对于中介分析中的分类变量，需要特别注意编码方式：

1. 二分变量：建议使用0/1编码而非1/2编码

   • 确保结果解释时logit值方向正确
   • 例如：relevel(factor(catholic), ref = "0")

2. 多分类变量：

   • 有序分类：可视为连续或使用多项式模型
   • 无序分类：必须创建虚拟变量

注意：当使用GLMM时，分类预测变量的参考水平设置会直接影响结果解释。建议在分析前使用contrasts()函数检查对比矩阵。

3. GLMM模型构建与验证

3.1 两阶段模型设定

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

1. 中介变量模型：预测中介变量(attachment)与自变量(catholic)的关系
2. 结果变量模型：预测结果变量(fight)与自变量、中介变量的关系

library(lme4)

# 模型1：中介变量模型 (M ~ X + covariates)
med_formula <- attachment ~ catholic + gender + income + (1|SCH_ID)
med.fit <- glmer(med_formula, 
                family = binomial(link = "logit"),
                data = student,
                control = glmerControl(optimizer = "bobyqa"))

# 模型2：结果变量模型 (Y ~ X + M + X*M + covariates)
out_formula <- fight ~ catholic * attachment + gender + income + (1 + attachment|SCH_ID)
out.fit <- glmer(out_formula,
                family = binomial(link = "logit"),
                data = student,
                control = glmerControl(optimizer = "bobyqa"))

3.2 模型诊断关键指标

在运行模型后，必须检查以下诊断指标：

1. 收敛性：确保模型已收敛

   • 检查optimizer警告
   • 尝试不同优化算法如bobyqa或Nelder_Mead

2. 奇异拟合：检查随机效应方差是否接近0

   • 使用isSingular()函数
   • 可能需要简化随机效应结构

3. 多重共线性：

   # 检查VIF值
   car::vif(med.fit)
   car::vif(out.fit)
   

4. 模型比较：

   • 使用ANOVA比较嵌套模型
   • 或通过AIC/BIC评估模型拟合

4. 中介效应分析与结果解读

4.1 mediation包实现

library(mediation)
set.seed(1234)  # 保证结果可重复

med.out <- mediate(
  med.fit, 
  out.fit,
  treat = "catholic",  # 处理变量
  mediator = "attachment",  # 中介变量
  sims = 1000,  # 建议至少1000次模拟
  boot = TRUE,  # 使用bootstrap
  boot.ci.type = "perc"  # 百分位置信区间
)

summary(med.out)

4.2 结果解读框架

典型输出包含三部分效应：

1. ACME (Average Causal Mediation Effect)：

   • 中介效应（间接效应）
   • 表示X通过M影响Y的效应量

2. ADE (Average Direct Effect)：

   • 直接效应
   • 表示X直接影响Y的部分

3. Total Effect：

   • 总效应 = ACME + ADE

对于分类变量，效应量以概率尺度报告。解读时需注意：

• 效应方向：系数符号表示效应方向
• 统计显著性：95%CI不包含0表示显著
• 效应大小：比较ACME与ADE的相对比例

4.3 可视化呈现

# 基础可视化
plot(med.out)

# 使用ggplot2增强可视化
library(ggplot2)
effects <- data.frame(
  Effect = c("ACME", "ADE", "Total"),
  Estimate = c(med.out$d0, med.out$z0, med.out$tau.coef),
  CI_lower = c(med.out$d0.ci[1], med.out$z0.ci[1], med.out$tau.ci[1]),
  CI_upper = c(med.out$d0.ci[2], med.out$z0.ci[2], med.out$tau.ci[2])
)

ggplot(effects, aes(x = Effect, y = Estimate)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = CI_lower, ymax = CI_upper), width = 0.1) +
  labs(title = "中介效应分解", y = "效应量估计", x = "") +
  theme_minimal()

5. 高级主题与问题排查

5.1 常见问题解决方案

1. 模型不收敛：

   • 增加迭代次数：control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e5))
   • 标准化连续预测变量
   • 简化随机效应结构

2. 零膨胀问题：

   • 考虑零膨胀模型
   • 使用brms包进行贝叶斯估计

3. 小样本校正：

   • 使用Kenward-Roger或Satterthwaite自由度近似
   • 考虑贝叶斯方法提供更稳定的估计

5.2 多类别中介分析

当自变量或中介变量为多分类时：

# 对于三分类中介变量
med.fit <- glmer(attachment_level ~ catholic + (1|SCH_ID),
                family = binomial(link = "logit"),
                data = student %>% 
                  mutate(attachment_level = case_when(
                    attachment == 0 ~ "low",
                    attachment == 1 ~ "medium",
                    TRUE ~ "high") %>% 
                      ordered()))

# 需要使用multinomial模型
library(nnet)
med.fit <- multinom(attachment_level ~ catholic, data = student)

5.3 敏感性分析

评估未测量混杂因素的影响：

# 使用mediation包的sens参数
med.sens <- medsens(med.out, rho.by = 0.1, effect.type = "indirect")
plot(med.sens)

6. 完整案例代码与数据获取

6.1 完整分析流程

# 加载必要包
library(tidyverse)
library(lme4)
library(mediation)

# 数据准备
student <- read_csv("student_data.csv") %>% 
  mutate(across(c(catholic, attachment, fight, gender), 
                ~ factor(.x, levels = c(0, 1), labels = c("No", "Yes"))))

# 模型构建
med.fit <- glmer(attachment ~ catholic + gender + income + (1|SCH_ID),
                family = binomial,
                data = student,
                control = glmerControl(optimizer = "bobyqa"))

out.fit <- glmer(fight ~ catholic * attachment + gender + income + (1|SCH_ID),
                family = binomial,
                data = student,
                control = glmerControl(optimizer = "bobyqa"))

# 中介分析
med.out <- mediate(med.fit, out.fit, 
                  treat = "catholic", mediator = "attachment",
                  sims = 1000, boot = TRUE)

# 结果输出
summary(med.out)
plot(med.out)

# 效应量转换
effects <- summary(med.out)
logit_to_prob <- function(x) exp(x)/(1+exp(x))
data.frame(
  Effect = c("ACME", "ADE", "Total"),
  Logit_Scale = c(effects$d0, effects$z0, effects$tau.coef),
  Probability_Scale = logit_to_prob(c(effects$d0, effects$z0, effects$tau.coef))
)

6.2 模拟数据生成

如果无法获取真实数据，可以使用以下代码生成模拟数据集：

set.seed(123)
n_schools <- 20
n_students <- 200

sim_data <- data.frame(
  SCH_ID = rep(1:n_schools, each = n_students/n_schools),
  catholic = rbinom(n_students, 1, 0.5),
  gender = rbinom(n_students, 1, 0.5),
  income = sample(1:13, n_students, replace = TRUE)
) %>% 
  group_by(SCH_ID) %>% 
  mutate(
    # 学校随机效应
    school_effect = rnorm(1, 0, 0.5),
    # 中介变量模型
    logit_attachment = -0.5 + 0.8*catholic + 0.3*gender - 0.1*income + school_effect,
    prob_attachment = plogis(logit_attachment),
    attachment = rbinom(n(), 1, prob_attachment),
    # 结果变量模型
    logit_fight = -1 + 0.5*catholic - 0.8*attachment + 0.4*catholic*attachment + 
      0.2*gender - 0.05*income + school_effect,
    prob_fight = plogis(logit_fight),
    fight = rbinom(n(), 1, prob_fight)
  ) %>% 
  ungroup()

write_csv(sim_data, "student_sim_data.csv")