如何在R中对含多因变量的大数据集执行三因素ANOVA？

问题背景 #

我拥有一个包含3个自变量（Sex、Tratment1、Treatment2）及150余个因变量（如var1、var2等）的数据集，希望为每个因变量执行三因素ANOVA，且仅展示显著结果。了解可使用for循环或apply系列函数，但不知如何在该数据集上应用；虽听说循环口碑不佳，但不确定其他处理方式。

示例数据集 #

df <- structure(list(ID = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 
29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39), Sex = c(1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), Tratment1 = c(1, 
1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 
0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1), Treatment2 = c(1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0), var1 = c(38, 
34, 97, 6, 99, 44, 4, 74, 12, 80, 61, 44, 100, 45, 31, 46, 4, 
92, 42, 52, 30, 55, 40, 31, 11, 7, 15, 67, 19, 14, 78, 19, 7, 
48, 39, 41, 29, 37, 10), var2 = c(38, 34, 47, 39, 83, 21, 66, 
87, 94, 55, 28, 52, 90, 96, 18, 36, 89, 77, 70, 49, 58, 16, 57, 
55, 50, 82, 88, 90, 54, 56, 27, 99, 50, 13, 94, 63, 84, 17, 71
), var3 = c(13, 11, 16, 38, 37, 36, 1, 26, 28, 61, 29, 12, 91, 
50, 34, 43, 100, 65, 14, 11, 23, 36, 31, 54, 72, 62, 31, 66, 
62, 73, 56, 31, 53, 23, 67, 86, 71, 96, 92), var4 = c(72, 78, 
32, 11, 1, 72, 95, 44, 23, 72, 67, 73, 100, 83, 2, 54, 67, 76, 
20, 79, 29, 75, 17, 96, 91, 35, 52, 71, 48, 15, 79, 6, 8, 57, 
78, 84, 87, 37, 31), var5 = c(22, 85, 60, 20, 53, 26, 65, 5, 
75, 20, 5, 47, 59, 13, 22, 90, 19, 51, 98, 93, 76, 17, 9, 5, 
57, 85, 21, 95, 40, 5, 30, 48, 83, 17, 90, 47, 72, 6, 36), var6 = c(34, 
9, 49, 7, 58, 55, 82, 65, 1, 36, 74, 48, 40, 26, 21, 37, 49, 
100, 36, 68, 78, 47, 77, 62, 38, 52, 36, 77, 46, 1, 99, 67, 40, 
43, 86, 17, 96, 33, 70), var7 = c(48, 94, 72, 13, 58, 4, 76, 
13, 45, 11, 41, 14, 36, 53, 2, 37, 1, 57, 23, 59, 7, 33, 48, 
52, 78, 64, 22, 63, 13, 56, 14, 28, 99, 74, 72, 39, 29, 49, 18
), var8 = c(98, 94, 53, 60, 72, 61, 62, 31, 28, 51, 94, 27, 8, 
46, 72, 91, 77, 2, 78, 45, 97, 39, 61, 12, 96, 59, 5, 60, 42, 
32, 42, 22, 50, 23, 32, 40, 18, 73, 52), var10 = c(90, 9, 51, 
88, 62, 36, 94, 96, 24, 21, 64, 82, 55, 33, 85, 98, 54, 68, 51, 
86, 22, 26, 60, 13, 29, 21, 4, 77, 64, 20, 28, 65, 97, 65, 61, 
20, 5, 74, 91)), row.names = c(NA, -39L), class = c("tbl_df", 
"tbl", "data.frame"))

解决方案 #

方法1：使用tidyverse生态（推荐，代码更整洁） #

利用purrr的map函数配合broom包整理ANOVA结果，全程保持数据框格式，便于筛选和查看。

# 加载所需包
library(tidyverse)
library(broom)

# 1. 将数据转为长格式，每个因变量单独一行
long_df <- df %>%
  select(-ID) %>%  # 移除不必要的ID列
  pivot_longer(cols = starts_with("var"), names_to = "dependent_var", values_to = "value")

# 2. 按因变量分组，执行三因素ANOVA，提取显著结果
significant_results <- long_df %>%
  group_by(dependent_var) %>%
  nest() %>%
  mutate(
    anova_model = map(data, ~aov(value ~ Sex * Tratment1 * Treatment2, data = .x)),
    anova_tidy = map(anova_model, tidy),
    # 筛选p值<0.05的项，同时保留因变量名
    significant_terms = map(anova_tidy, ~filter(.x, p.value < 0.05))
  ) %>%
  # 移除没有显著项的因变量
  filter(map_lgl(significant_terms, ~nrow(.x) > 0)) %>%
  unnest(significant_terms) %>%
  select(dependent_var, term, p.value)

# 查看结果
print(significant_results)

方法2：基础R的lapply实现 #

无需额外包，用基础R的lapply遍历所有因变量列，收集结果后筛选。

# 提取所有因变量列名
dep_vars <- grep("^var", names(df), value = TRUE)

# 定义ANOVA函数
run_anova <- function(dep_var) {
  formula <- as.formula(paste(dep_var, "~ Sex * Tratment1 * Treatment2"))
  model <- aov(formula, data = df)
  # 提取方差分析表
  anova_table <- summary(model)[[1]]
  # 整理成数据框，保留显著项
  res <- data.frame(
    dependent_var = dep_var,
    term = rownames(anova_table),
    p.value = anova_table[, "Pr(>F)"]
  ) %>%
  filter(p.value < 0.05)
  return(res)
}

# 批量执行并合并结果
significant_results <- do.call(rbind, lapply(dep_vars, run_anova))

# 查看结果
print(significant_results)

方法3：for循环（完全可用，不必排斥） #

循环并非洪水猛兽，对于这种明确的批量任务，逻辑清晰易调试，处理150个变量完全没问题。

# 提取因变量列名
dep_vars <- grep("^var", names(df), value = TRUE)
# 初始化空数据框存储结果
significant_results <- data.frame()

for(var in dep_vars) {
  # 构建公式
  formula <- as.formula(paste(var, "~ Sex * Tratment1 * Treatment2"))
  # 执行ANOVA
  model <- aov(formula, data = df)
  anova_table <- summary(model)[[1]]
  # 筛选显著项
  sig_rows <- anova_table[, "Pr(>F)"] < 0.05
  if(any(sig_rows)) {
    temp <- data.frame(
      dependent_var = var,
      term = rownames(anova_table)[sig_rows],
      p.value = anova_table[sig_rows, "Pr(>F)"]
    )
    significant_results <- rbind(significant_results, temp)
  }
}

# 查看结果
print(significant_results)

注意事项 #

• 显著性水平默认设为0.05，可根据需求修改p.value < 0.05中的阈值。
• 若需要校正多重比较（比如150个因变量容易出现假阳性），可在提取p值后添加p.adjust()函数，例如p.value = p.adjust(p.value, method = "bonferroni")。
• 确保自变量Sex、Tratment1、Treatment2是因子类型（若当前是数值型，可先用df$Sex <- as.factor(df$Sex)转换），否则ANOVA结果可能不准确。

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