如何针对多结局变量运行混合线性回归并生成规范结果？

问题背景 #

现有包含3个组别、2个时间点、49个结局变量的纵向数据，需要针对每个结局变量分析组内不同时间点的变化（采用混合线性回归，纳入个体随机截距）。手动完成147次分析不现实，现有代码能生成模型列表，但无法提取p值，直接用tidy()处理模型列表时报错：Error: No tidy method recognized for this list.

解决方案 #

核心思路 #

• 避免拆分数据集，直接通过映射处理组别+结局变量的所有组合
• 使用broom.mixed包的tidy()函数提取模型结果（适配lmerTest输出的模型，可直接获取p值）
• 将所有结果整合成结构化表格，便于查看和后续分析

完整代码 #

# 加载所需包
library(lme4)
library(lmerTest)
library(tidyverse)
library(broom.mixed)

# 模拟示例数据（替换为你的真实数据）
set.seed(123) # 保证结果可重复
df <- data.frame(
  id = rep(1:66, each = 2),
  visit = rep(0:1, 66),
  rand = rep(0:2, each = 44), # 修正分组逻辑，确保每个组含44个个体（共88行数据）
  x1 = sample(4000:9000, 132),
  x2 = sample(1200:3400, 132),
  x3 = sample(220:400, 132)
)

# 自动提取结局变量和组别列表
outcome_vars <- colnames(df)[4:ncol(df)] # 适配你的49个结局变量
groups <- unique(df$rand)

# 批量运行模型并整理结果
results <- crossing(group = groups, outcome = outcome_vars) %>%
  mutate(
    # 对每个组别+结局变量组合运行混合线性回归
    model = map2(group, outcome, function(g, y) {
      lmer(
        formula = paste(y, "~ visit + (1|id)"),
        data = df %>% filter(rand == g),
        na.action = na.omit
      )
    }),
    # 提取模型固定效应结果（含p值、置信区间）
    model_tidy = map(model, tidy, effects = "fixed", conf.int = TRUE)
  ) %>%
  # 展开结果，仅保留时间点（visit）的效应结果
  unnest(model_tidy) %>%
  filter(term == "visit") %>%
  select(group, outcome, estimate, std.error, statistic, p.value, conf.low, conf.high)

# 查看最终汇总表格
print(results)

关键说明 #

1. 分组逻辑优化：用crossing()生成所有组别与结局变量的组合，无需手动拆分数据集，减少冗余代码
2. p值提取：broom.mixed::tidy()专门适配lmerTest输出的模型对象，可直接提取固定效应的p值、置信区间等核心指标
3. 结果精简：最终表格仅保留我们关注的visit项结果，包含组别、结局变量、效应量、统计量、p值和置信区间，结构清晰
4. 原报错原因：直接对模型列表调用tidy()（默认调用tidymodels方法）不支持，需用broom.mixed的tidy()逐个处理模型对象

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