在R中快速计算ANOVA

Question

我有一个具有以下维度的数据框：

dim(b)  
[1]    974 433685

这些列代表我想在上面运行方差分析的变量（即，我想运行433,685个方差分析）。样本大小为974。最后一列是“组”变量。

我想出了3种不同的方法，但由于测试数量太多，它们都太慢了。

首先，让我们生成一个小的练习数据集来玩：

dat = as.data.frame(matrix(runif(10000*500), ncol = 10000, nrow = 500))
dat$group = rep(letters[1:10], 5000)

方法1（基于'sapply'）：

system.time(sapply(dat[,-length(dat)], function(x) aov(x~group, data=dat) ))

   user  system elapsed 
 143.76    0.33  151.79 

方法二（基于“parallel”包中的“mclapply”）：

library(parallel)
options(mc.cores=3)
system.time(mclapply(dat[,-length(dat)], function(x) aov(x~group, data=dat) ))

   user  system elapsed 
 141.76    0.21  142.58 

第三种方法（基于 'cbind' 绑定左侧）：

formula = as.formula( paste0("cbind(", paste(names(dat)[-length(dat)],collapse=","), ")~group") ) 
system.time(aov(formula, data=dat))

  user  system elapsed 
  10.00    0.22   10.25 

在实践数据集中，方法3是明显的赢家。然而，当我在我的实际数据上执行此操作时，仅使用方法3计算10个（共433,685个）列需要这么长时间：

   user  system elapsed
119.028   5.430 124.414

我不确定为什么在我的实际数据上需要更长时间。我可以访问一个拥有16个核心和72GB RAM的Linux集群。

有没有任何方法可以更快地计算这个问题？

- Chad Johnson

---

"433,685个ANOVA？你这样做的目的是什么？对于你的问题来说肯定有更好的统计方法。" - Roland

这是一个经过Bonferonni校正的p值，为1.152911e-07。 - Paul Lemmens

Roland，我正在尝试通过我的“组”变量来量化“批处理效应”，即来自不需要的技术因素的变异。我有433K个单个探针来测量特定基因组位点的DNA甲基化。我的想法是将每个探针的ANOVA的单个F统计量相加。然后，我会将这个数字与从不同预处理流程生成的类似数据集的数字进行比较，以找到去除更多“批处理效应”的那个。 - Chad Johnson

Answer 1

对于使用相同的 设计矩阵 同时拟合许多通用线性模型（如方差分析），Bioconductor/R limma package 提供了一个非常快速的 lmFit() 函数。以下是使用 limma 拟合 ANOVA 模型的方法：

library(limma)

# generate some data 
# (same dimensions as in your question)
nrows <- 1e4
ncols <- 5e2
nlevels <- 10
dat <- matrix(
  runif(nrows * ncols), 
  nrow = nrows, 
  ncol = ncols
)
group <- factor(rep(
  letters[1:nlevels], 
  ncols / nlevels
))

# construct the design matrix
# (same as implicitly used in your question)
dmat <- model.matrix(~ group)
# fit the ANOVA model
fit <- lmFit(dat, dmat)

在我的笔记本电脑上，它在与您问题中的数据相同维度的数据上完成时间为0.4-0.45秒。

- Johann de Jong