如何加速大数据量下多变量组合的线性回归计算？

我有一个包含36365760行、10列的data.frame，示例数据如下：

dat3 <- data.frame("Region"=rep(c("R1","R2","R3","R1","R2"),20),
                   "Phase"=rep(c("S1","S2"),50),
                   "Treatment"=rep(c("P","D"),50),
                   "Region_ID"=rep(1:2,50),
                   "Signal"=rnorm(100),
                   "Bin"=rep(1,100))

随后我针对Region、Phase、Treatment、Region_ID的每种组合拟合线性模型，代码如下：

res <- lapply(unique(dat3$Region), function(i){
  lapply(unique(dat3$Phase), function(j){
    lapply(unique(dat3$Treatment),function(k){
      lapply(unique(dat3$Region_ID[dat3$Region == i & dat3$Phase == j & dat3$Treatment == k]),function(l){
        y=dat3[dat3$Region==i & dat3$Phase==j & dat3$Treatment ==k & dat3$Region_ID ==l & dat3$Bin %in% c(1:10,90:100),]$Signal
        x=dat3[dat3$Region==i & dat3$Phase==j & dat3$Treatment ==k & dat3$Region_ID ==l & dat3$Bin %in% c(1:10,90:100),]$Bin
        lm(y~x)
      })
    })
  })
})

这段代码在计算集群上运行一晚上都未完成，但子集化完整的data.frame时可以正常运行，求优化方案。

1. 提前过滤数据，减少重复运算 #

每次嵌套循环都重复筛选Bin %in% c(1:10,90:100)的行，会产生大量冗余计算。先一次性过滤出符合条件的数据，后续所有操作基于这个子集：

# 提前过滤Bin条件，缩小处理数据集
filtered_dat <- dat3[dat3$Bin %in% c(1:10, 90:100), ]

2. 用分组操作替代嵌套lapply #

嵌套lapply会反复执行相同的分组筛选，效率极低。推荐用dplyr或data.table的分组功能，自动处理所有分组组合，内部优化了数据访问逻辑：

基于dplyr的实现 #

library(dplyr)
library(tidyr)
library(purrr)

model_results <- filtered_dat %>%
  group_by(Region, Phase, Treatment, Region_ID) %>%
  nest() %>% # 将每组数据打包成列表列
  mutate(model = map(data, ~lm(Signal ~ Bin, data = .x))) %>% # 对每组拟合模型
  ungroup()

基于data.table的实现（更适合大数据） #

data.table的分组运算速度和内存效率远高于原生data.frame，适合千万级数据：

library(data.table)
setDT(filtered_dat)

model_results <- filtered_dat[, 
  .(model = list(lm(Signal ~ Bin, .SD))), 
  by = .(Region, Phase, Treatment, Region_ID)
]

3. 替换lm为轻量拟合函数 #

lm会生成大量冗余信息（如残差、拟合值），如果只需要系数、R²等核心结果，用更高效的函数替代：

用fastLm加速拟合 #

RcppEigen包的fastLm比原生lm快数倍，适合批量小模型拟合：

library(RcppEigen)

# 结合data.table使用
model_results <- filtered_dat[, 
  .(model = list(fastLm(Signal ~ Bin, .SD))), 
  by = .(Region, Phase, Treatment, Region_ID)
]

手动计算回归系数（最快方案） #

如果只需要截距和斜率，直接用矩阵运算跳过模型对象，速度和内存占用最优：

model_results <- filtered_dat[, {
  X <- cbind(1, Bin) # 构造含截距项的设计矩阵
  beta <- solve(crossprod(X), crossprod(X, Signal)) # 最小二乘解
  list(intercept = beta[1], slope = beta[2])
}, by = .(Region, Phase, Treatment, Region_ID)]

4. 并行计算利用集群多核资源 #

集群的核心优势是多线程并行，将分组任务分配到多个核心同时执行：

基于furrr的并行实现（dplyr生态） #

library(furrr)
plan(multisession) # 根据集群环境调整，如multicore、cluster

model_results <- filtered_dat %>%
  group_by(Region, Phase, Treatment, Region_ID) %>%
  nest() %>%
  mutate(model = future_map(data, ~lm(Signal ~ Bin, data = .x))) %>%
  ungroup()

基于parallel+data.table的并行实现 #

library(parallel)
cl <- makeCluster(detectCores()) # 根据集群可用核心数设置
clusterExport(cl, c("filtered_dat")) # 将数据导出到集群节点

model_results <- filtered_dat[, 
  .(model = list(parLapply(cl, .SD, function(x) lm(Signal ~ Bin, x)))), 
  by = .(Region, Phase, Treatment, Region_ID)
]

stopCluster(cl)

5. 内存优化 #

3600万行的数据内存压力大，通过以下方式减少内存占用：

• 将data.frame转为data.table，内存占用可降低约50%；
• 将字符型分组列转为因子，减少内存消耗：

setDT(dat3)
dat3[, c("Region", "Phase", "Treatment") := lapply(.SD, as.factor), .SDcols = c("Region", "Phase", "Treatment")]

• 处理完后及时删除大对象并释放内存：

rm(dat3)
gc()

6. 提前生成有效分组组合 #

嵌套循环中反复计算unique(dat3$Region_ID[...])会重复筛选数据，提前生成所有非空的分组组合，避免冗余：

# 提前获取所有存在的分组组合
valid_groups <- unique(filtered_dat[, .(Region, Phase, Treatment, Region_ID)])

# 遍历每个分组拟合模型
res <- lapply(1:nrow(valid_groups), function(row_idx){
  group <- valid_groups[row_idx, ]
  sub_dat <- filtered_dat[Region == group$Region & Phase == group$Phase & 
                           Treatment == group$Treatment & Region_ID == group$Region_ID, ]
  lm(Signal ~ Bin, data = sub_dat)
})

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