在R中按Pair分组对多列数据执行线性回归与ANOVA分析

我来帮你搞定这个头疼的多分组多变量回归问题！针对你这种7056个pair水平、100个V列的场景，最省心且高效的方法是用tidyverse生态工具链（适合代码可读性优先），或者data.table（适合大数据集速度优先），毕竟手动循环很容易出错，而且效率低下。

一、推荐方案1：tidyverse + broom（代码直观，结果整洁） #

这个方法用dplyr做分组、purrr做模型映射、broom把复杂的模型输出转换成标准DataFrame，完美适配你的需求：

步骤1：加载依赖包 #

library(tidyverse)
library(broom) # 专门用来整理模型输出的神器

步骤2：把宽格式数据转成长格式 #

把你的V1-V100从多列转成“变量名+变量值”的两列，这样就能用分组操作覆盖所有V列：

df_long <- df %>%
  pivot_longer(
    cols = starts_with("V"), # 筛选所有以V开头的列
    names_to = "var_name",   # 新列：存储V1/V2...
    values_to = "var_value"  # 新列：存储对应V列的数值
  )

步骤3：分组拟合模型并提取结果 #

按pair和var_name（即每个V列）分组，分别拟合lm(var_value ~ g)，同时提取回归系数和ANOVA统计量：

# 提取线性回归的系数、p值等核心结果
reg_results <- df_long %>%
  group_by(pair, var_name) %>%
  nest() %>% # 按分组把数据嵌套成列表，方便后续拟合模型
  mutate(
    # 对每个分组数据拟合线性回归
    model = map(data, ~ lm(var_value ~ g, data = .x)),
    # 把模型系数转换成整洁的DataFrame
    tidy_coefs = map(model, tidy),
    # 把ANOVA结果转换成整洁的DataFrame
    tidy_anova = map(model, ~ tidy(anova(.x)))
  ) %>%
  # 展开系数结果，得到每行对应一个分组、一个变量、一个系数的结构
  unnest(tidy_coefs) %>%
  select(pair, var_name, term, estimate, std.error, statistic, p.value)

# 如果需要单独查看ANOVA结果，可以单独展开tidy_anova
anova_results <- df_long %>%
  group_by(pair, var_name) %>%
  nest() %>%
  mutate(
    model = map(data, ~ lm(var_value ~ g, data = .x)),
    tidy_anova = map(model, ~ tidy(anova(.x)))
  ) %>%
  unnest(tidy_anova)

这样得到的reg_results和anova_results都是标准的DataFrame，你可以直接导出成CSV或者做后续分析，完全不需要手动处理复杂的模型输出。

二、推荐方案2：data.table（大数据集速度优先） #

如果你的数据集特别大（7056组+100列确实不小），data.table的分组操作速度会比tidyverse快很多，内存占用也更高效：

步骤1：加载依赖包 #

library(data.table)
library(broom)

步骤2：转长格式+分组拟合 #

setDT(df) # 把普通DataFrame转换成data.table格式

# 转长格式
dt_long <- melt(
  df,
  id.vars = c("pair", "g"), # 保留分组列和因变量列
  variable.name = "var_name",
  value.name = "var_value"
)

# 按pair和var_name分组，拟合模型并提取结果
dt_results <- dt_long[, 
  {
    mod <- lm(var_value ~ g)
    # 分别提取回归系数和ANOVA结果
    list(
      coefs = list(tidy(mod)),
      anova = list(tidy(anova(mod)))
    )
  },
  by = .(pair, var_name) # 分组依据
]

# 展开结果成扁平的DataFrame
final_coefs <- rbindlist(dt_results$coefs, idcol = c("pair", "var_name"))
final_anova <- rbindlist(dt_results$anova, idcol = c("pair", "var_name"))

三、你之前代码的问题梳理 #

你之前的循环没得到理想结果，主要是这几个问题：

1. 用索引3:ncol(df)来定位V列很不稳定，如果列顺序变了就会出错，不如直接用列名筛选（比如starts_with("V")）更可靠；
2. lapply(levels(df$pair), function(x) {df=x...})这里的df=x是完全错误的，应该先筛选出当前pair的子集再拟合模型；
3. 手动处理coef(summary(lm(...)))的输出很容易出错，broom::tidy()能自动把模型输出转换成标准结构，省掉很多麻烦。

小提示：加速处理 #

如果想进一步提速，可以用furrr包开启并行计算，把map()换成future_map()，利用电脑的多核心同时处理多个分组，大幅节省时间。

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