R中的快速方差分析计算

Question

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

dim(b)  
[1]    974 433685

列表示我想运行方差分析的变量(即，我想运行433,685个方差分析)。样本量为974。最后一列是"group"变量。

我已经提出了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 

方法2(基于'mclapply'从'parallel'包):

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 

方法3(基于'cbind'-ing LHS):

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集群。

是否有更快的计算方法?

Answer 1

为了使用相同的设计矩阵同时拟合许多一般线性模型(如ANOVA)， Bioconductor/R limma包提供了非常快速的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秒内完成，数据的维度与您的问题中的数据相同。