在R中基于两个行匹配的data.frame执行多重t检验的实现方案

现有两个按id匹配的data.frame（df1、df2）：

• df1包含g1-g10共10个分组列，取值为0、1或NA
• df2包含b1-b17共17个数值列，存在NA值

需要对每个g列，检验其0组和1组在每个b列的数值是否存在显著差异，最终生成包含以下字段的结果data.frame df3：

• g：分组列名（如g1、g2）
• b：数值列名（如b1、b2）
• mean_0：0组对应b列的均值（忽略NA）
• mean_1：1组对应b列的均值（忽略NA）
• p.val：t检验的p值
• p.adjust：按每个g列校正后的p值

并将df3按g列和p.val排序。

以下提供两种实现方式，可根据习惯选择：

方式一：使用tidyverse（简洁高效） #

# 加载包
library(tidyverse)

# 合并df1和df2，保留id用于匹配
combined_df <- inner_join(df1, df2, by = "id")

# 转换为长格式，方便分组处理
long_df <- combined_df %>%
  pivot_longer(cols = starts_with("g"), names_to = "g", values_to = "group") %>%
  pivot_longer(cols = starts_with("b"), names_to = "b", values_to = "value") %>%
  filter(group %in% c(0, 1))  # 仅保留有效分组（0/1），排除NA分组

# 批量执行t检验并整理结果
df3 <- long_df %>%
  group_by(g, b) %>%
  summarise(
    mean_0 = mean(value[group == 0], na.rm = TRUE),
    mean_1 = mean(value[group == 1], na.rm = TRUE),
    # 处理组内无有效数据的情况，避免t.test报错
    p.val = ifelse(
      all(is.na(value[group == 0])) | all(is.na(value[group == 1])),
      NA,
      t.test(value ~ group, na.action = na.omit)$p.value
    ),
    .groups = "drop"
  ) %>%
  # 按每个g列进行FDR校正，可修改method参数更换校正方法（如"bonferroni"）
  group_by(g) %>%
  mutate(p.adjust = p.adjust(p.val, method = "fdr")) %>%
  ungroup() %>%
  # 按分组列和p值排序
  arrange(g, p.val)

方式二：使用基础R循环（无需额外包） #

# 获取g列和b列的名称
g_cols <- names(df1)[startsWith(names(df1), "g")]
b_cols <- names(df2)[startsWith(names(df2), "b")]

# 初始化结果列表
result_list <- list()

# 遍历每个分组列g
for (g in g_cols) {
  group_vec <- df1[[g]]
  # 筛选出分组为0和1的有效行索引
  idx_0 <- which(group_vec == 0 & !is.na(group_vec))
  idx_1 <- which(group_vec == 1 & !is.na(group_vec))
  
  # 遍历每个数值列b
  for (b in b_cols) {
    val_0 <- df2[[b]][idx_0]
    val_1 <- df2[[b]][idx_1]
    
    # 计算两组均值
    mean_0 <- mean(val_0, na.rm = TRUE)
    mean_1 <- mean(val_1, na.rm = TRUE)
    
    # 计算p值，处理无有效数据的情况
    if (all(is.na(val_0)) | all(is.na(val_1))) {
      p_val <- NA
    } else {
      p_val <- t.test(val_0, val_1, na.rm = TRUE)$p.value
    }
    
    # 将结果存入列表
    result_list[[length(result_list) + 1]] <- data.frame(
      g = g,
      b = b,
      mean_0 = mean_0,
      mean_1 = mean_1,
      p.val = p_val,
      stringsAsFactors = FALSE
    )
  }
}

# 合并结果为data.frame
df3 <- do.call(rbind, result_list)

# 按每个g列进行FDR校正
df3 <- df3 %>%
  group_by(g) %>%
  mutate(p.adjust = p.adjust(p.val, method = "fdr")) %>%
  ungroup() %>%
  # 排序
  arrange(g, p.val)

关键说明 #

• 两种方式都自动处理了分组内全为NA的情况，避免t.test函数报错
• p值校正默认使用FDR方法，可通过修改p.adjust函数的method参数更换为其他方法（如"bonferroni"）
• 最终结果按分组列g和p值p.val升序排列，便于快速查看最显著的差异结果

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