如何在R中用emmeans获取线性混合模型中Visit的平均效应

问题背景 #

用lme4构建含重复测量的线性混合模型时，需要确定Visit（时间）对结局的效应（每次随访的结局变化量）。但模型包含多个与Visit的交互项，固定效应表中的Visit系数仅对应女性、20-24岁年龄组、暴露量为0的亚组，无法反映整体平均效应。尝试用emmeans包时仅得到Visit=3的边际均值，需要正确提取Visit的平均效应。

示例数据与模型 #

library(lme4)
library(emmeans)

# 示例数据
example <-  data.frame(
  "subject" = c(rep(1,5), rep(2,5), rep(3, 5)),
  "outcome" = c(15,13,10,10,9,12,12,11,9,9,18,14,12,11,10), 
  "visit" = c(rep(sequence(5), 3)), 
  "sex"= c(rep("male", 10), rep("female", 5)), 
  current_age = c(rep("20-24", 3), rep("25-29", 2), rep("30-34", 4), rep("35-39", 1), rep("20-24", 1), rep("25-29", 4)), 
  "exposure" = c(10, 9, 9, 8, 6, 12, 12, 12, 11, 10, 13, 12, 11, 10, 9)
)

# 混合模型（忽略奇异拟合错误）
example_mod <- lmer(outcome ~ (visit | subject) + visit + sex + visit:current_age + visit:exposure, data = example)

模型固定效应摘要：

Fixed effects:
                       Estimate Std. Error t value
(Intercept)            20.08153    2.60105   7.721
visit                  -0.59361    0.87432  -0.679
sexmale                -4.39004    2.99455  -1.466
visit:current_age25-29  0.11279    0.35448   0.318
visit:current_age30-34  1.42966    0.58610   2.439
visit:current_age35-39  1.34885    0.51105   2.639
visit:exposure         -0.18092    0.08711  -2.077

解决方案：用emtrends提取平均斜率 #

要获取Visit的平均效应（即每次随访结局的平均变化量），需要提取Visit作为连续变量的斜率，而非特定Visit值的边际均值，这里用emtrends()函数实现：

1. 获取所有亚组的Visit斜率 #

emtrends()会计算每个协变量组合下，Visit每增加1单位时结局的变化量（斜率）：

# 按性别、年龄组、暴露水平分组计算Visit斜率
visit_subgroup_trends <- emtrends(example_mod, ~ sex + current_age + exposure, var = "visit")
visit_subgroup_trends

输出会包含每个亚组的Visit斜率、标准误等信息。

2. 计算所有亚组的平均斜率 #

通过summary()指定by = NULL，对所有亚组的斜率求平均，得到整体平均时间效应：

# 未加权平均（各亚组权重相同）
average_trend_unweighted <- summary(visit_subgroup_trends, by = NULL)
print(average_trend_unweighted)

# 按样本中各亚组的比例加权平均
average_trend_weighted <- summary(visit_subgroup_trends, by = NULL, weights = "proportional")
print(average_trend_weighted)

• 未加权平均：所有协变量组合的斜率取算术平均
• 加权平均：按样本中各亚组的出现比例加权，更贴近实际数据分布

3. 仅按部分协变量分组求平均 #

如果只需要按某几个协变量分组的平均效应（比如仅按性别），可以指定by参数：

# 按性别分组的平均Visit效应
sex_average_trend <- summary(visit_subgroup_trends, by = "sex")
print(sex_average_trend)

其他方法：手动计算平均效应 #

如果更倾向手动实现，可以用predict()生成所有协变量组合的预测值，再计算斜率后平均：

# 生成所有协变量组合的预测数据
new_data <- expand.grid(
  visit = 1:5,
  sex = unique(example$sex),
  current_age = unique(example$current_age),
  exposure = unique(example$exposure),
  subject = 1  # 个体随机效应设为0，仅计算固定效应部分
)

# 生成预测值
new_data$pred <- predict(example_mod, newdata = new_data, re.form = ~0)

# 计算每个协变量组合的Visit斜率
library(dplyr)
slopes <- new_data %>%
  group_by(sex, current_age, exposure) %>%
  summarize(slope = lm(pred ~ visit)$coefficients[2], .groups = "drop")

# 计算平均斜率
average_slope <- mean(slopes$slope)
weighted_average_slope <- weighted.mean(slopes$slope, w = table(interaction(example$sex, example$current_age, example$exposure)))

内容的提问来源于stack exchange，提问作者kjb17