如何从lme模型列表批量提取time:grup_int交互项p值？

核心思路 #

通过自定义函数提取单个模型中目标交互项的p值，再批量遍历模型列表，最终整合为带因变量名称的数据框，方便导出。

示例代码（含模拟数据） #

1. 准备环境与模拟数据

library(nlme)

# 模拟多因变量的混合效应数据
set.seed(123)
dat <- expand.grid(id = 1:50, time = 0:2, grup_int = c("A", "B"))
dat$y1 <- rnorm(nrow(dat), 1 + 0.5*dat$time + 0.3*dat$time*(dat$grup_int=="B"), 1)
dat$y2 <- rnorm(nrow(dat), 2 + 0.2*dat$time + 0.4*dat$time*(dat$grup_int=="B"), 1)
dat$y3 <- rnorm(nrow(dat), 0.5 + 0.1*dat$time + 0.1*dat$time*(dat$grup_int=="B"), 1)

# 批量拟合模型并存储到列表（列表元素命名为因变量名）
dv_names <- c("y1", "y2", "y3")
model_list <- lapply(dv_names, function(dv) {
  lme(formula = as.formula(paste(dv, "~ time * grup_int")), 
      random = ~1 | id, 
      data = dat)
})
names(model_list) <- dv_names

2. 提取交互项p值（两种实现方式）

方式一：使用base R

# 定义提取函数
get_interaction_p <- function(model, target_term = "time:grup_int") {
  t_table <- summary(model)$tTable
  # 检查目标项是否存在，不存在返回NA并警告
  if (!target_term %in% rownames(t_table)) {
    warning(paste("交互项", target_term, "在模型中未找到"))
    return(NA)
  }
  t_table[target_term, "p-value"]
}

# 批量提取并转换为数据框
p_results <- data.frame(
  因变量名称 = names(model_list),
  交互项p值 = unlist(lapply(model_list, get_interaction_p)),
  stringsAsFactors = FALSE
)

# 查看结果
print(p_results)

# 导出为CSV
write.csv(p_results, "交互项p值结果.csv", row.names = FALSE)

方式二：使用tidyverse（更简洁）

library(tidyverse)

# 批量提取并整合
p_results <- map_dfr(model_list, get_interaction_p, .id = "因变量名称") %>%
  rename(交互项p值 = value)

# 导出
write_csv(p_results, "交互项p值结果_tidy.csv")

关键注意事项 #

• 确认交互项的准确名称：先通过summary(model_list[[1]])$tTable查看单个模型的系数表，确认交互项的行名（比如若grup_int是因子，交互项可能显示为time:grup_intB，需对应修改target_term参数）。
• 缺失处理：函数会对不含目标交互项的模型返回NA并抛出警告，方便排查异常模型。

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