7

我在研究中遇到的一个反复出现的分析范式是需要根据所有不同的组 id 值进行子集化，依次对每个组执行统计分析，并将结果放入输出矩阵中以供进一步处理/汇总。

我通常如何在 R 中执行此操作类似于以下内容：

data.mat <- read.csv("...")  
groupids <- unique(data.mat$ID)  #Assume there are then 100 unique groups
  
results <- matrix(rep("NA",300),ncol=3,nrow=100)  

for(i in 1:100) {  
  tempmat <- subset(data.mat,ID==groupids[i])  

  # Run various stats on tempmat (correlations, regressions, etc), checking to  
  # make sure this specific group doesn't have NAs in the variables I'm using  
  # and assign results to x, y, and z, for example.  

  results[i,1] <- x  
  results[i,2] <- y  
  results[i,3] <- z  
}

这最终对我有用，但根据数据的大小和我正在使用的组的数量，这可能需要长达三天的时间。

除了扩展到并行处理之外，还有什么“技巧”可以让这样的东西运行得更快吗？例如，将循环转换为其他东西（比如应用包含我想在循环内运行的统计信息的函数），或者消除将数据子集实际分配给变量的需要？

编辑：

也许这只是常识（或抽样错误），但我尝试在我的一些代码中使用括号进行子集化，而不是使用子集命令，它似乎提供了轻微的性能提升，这让我感到惊讶。我有一些使用与上面相同的对象名称在下面输出的代码：

system.time(for(i in 1:1000){data.mat[data.mat$ID==groupids[i],]})  

   user  system elapsed  
 361.41   92.62  458.32

system.time(for(i in 1:1000){subset(data.mat,ID==groupids[i])})  

   user  system elapsed   
 378.44  102.03  485.94

更新：

在其中一个答案中，jorgusch 建议我使用 data.table 包来加快我的子集化。因此，我将其应用于本周早些时候遇到的一个问题。在一个有超过 1,500,000 行和 4 列（ID、Var1、Var2、Var3）的数据集中，我想计算每组中的两个相关性（由“ID”变量索引）。有50,000多个组。下面是我的初始代码（与上面非常相似）：

data.mat <- read.csv("//home....")  
groupids <- unique(data.mat$ID)
  
results <- matrix(rep("NA",(length(groupids) * 3)),ncol=3,nrow=length(groupids))  

for(i in 1:length(groupids)) {  
  tempmat <- data.mat[data.mat$ID==groupids[i],] 

  results[i,1] <- groupids[i]  
  results[i,2] <- cor(tempmat$Var1,tempmat$Var2,use="pairwise.complete.obs")  
  results[i,3] <- cor(tempmat$Var1,tempmat$Var3,use="pairwise.complete.obs")    

}

我现在正在重新运行它，以准确衡量这需要多长时间，但据我所知，我早上进入办公室时开始运行它，并在下午中旬的某个时间结束。图 5-7 小时。

重组我的代码以使用 data.table....

data.mat <- read.csv("//home....")  
data.mat <- data.table(data.mat)  
  
testfunc <- function(x,y,z) {  
  temp1 <- cor(x,y,use="pairwise.complete.obs")  
  temp2 <- cor(x,z,use="pairwise.complete.obs")  
  res <- list(temp1,temp2)  
  res  
}  

system.time(test <- data.mat[,testfunc(Var1,Var2,Var3),by="ID"])  

 user  system  elapsed  
16.41    0.05    17.44

将使用 data.table 的结果与我使用 for 循环对所有 ID 进行子集化并手动记录结果的结果进行比较，它们似乎给了我相同的答案（尽管我必须更彻底地检查一下）。这看起来是一个相当大的速度提升。

更新 2：

使用子集运行代码终于又完成了：

   user     system   elapsed  
17575.79  4247.41   23477.00

更新 3：

我想看看使用同样推荐的 plyr 包是否有任何不同的结果。这是我第一次使用它，所以我的工作效率可能有些低，但与带有子集的 for 循环相比，它仍然有很大帮助。

使用与以前相同的变量和设置...

data.mat <- read.csv("//home....")  
system.time(hmm <- ddply(data.mat,"ID",function(df)c(cor(df$Var1,df$Var2,  use="pairwise.complete.obs"),cor(df$Var1,df$Var3,use="pairwise.complete.obs"))))  

  user  system elapsed  
250.25    7.35  272.09
4

4 回答 4

6

这几乎正​​是该plyr软件包旨在简化的内容。然而，它不太可能让事情变得更快——大部分时间可能都花在了统计上。

于 2010-03-27T01:31:15.787 回答
3

此外plyr，您可以尝试使用foreach包来排除显式循环计数器，但我不知道它是否会给您带来任何性能优势。

Foreach，不过，如果您有多核工作站（带有doMC/multicore包），则为您提供一个非常简单的并行块处理接口（查看DoMC 和 foreach 入门了解详细信息）学生。如果不是唯一的原因，plyr恕我直言是很好的解决方案。

于 2010-03-27T11:59:42.707 回答
2

就个人而言，我发现 plyr 不是很容易理解。我更喜欢 data.table 这也更快。例如，您想为每个 ID 计算 colum my_column 的标准差。

dt <- datab.table[df] # one time operation...changing format of df to table
result.sd <- dt[,sd(my_column),by="ID"] # result with each ID and SD in second column

三个这样的语句，最后是一个 cbind - 这就是你所需要的。您还可以使用 dt 仅对一个 ID 执行一些操作，而无需使用新语法中的子集命令：

result.sd.oneiD<- dt[ID="oneID",sd(my_column)]

第一条语句指的是行 (i)，第二条语句指的是列 (j)。

如果发现播放器更容易阅读并且更灵活，因为您还可以在“子集”中创建子域......文档描述它使用类似 SQL 的方法。例如，by 几乎是 SQL 中的“分组依据”。好吧，如果您了解 SQL，您可能可以做更多事情，但没有必要使用包。最后，它非常快，因为每个操作不仅是并行的，而且 data.table 会抓取计算所需的数据。然而，子集保持整个矩阵的水平并将其拖过内存。

于 2010-03-27T17:23:30.270 回答
2

您已经建议矢量化并避免对中间结果进行不必要的复制，因此您肯定走在正确的轨道上。让我提醒你不要做我所做的，只是假设向量化总是会给你带来性能提升（就像它在其他语言中所做的那样，例如 Python + NumPy、MATLAB）。

一个例子：

# small function to time the results:
time_this = function(...) {
  start.time = Sys.time(); eval(..., sys.frame(sys.parent(sys.parent()))); 
  end.time = Sys.time(); print(end.time - start.time)
}

# data for testing: a 10000 x 1000 matrix of random doubles
a = matrix(rnorm(1e7, mean=5, sd=2), nrow=10000)

# two versions doing the same thing: calculating the mean for each row
# in the matrix
x = time_this( for (i in 1:nrow(a)){ mean( a[i,] ) } )
y = time_this( apply(X=a, MARGIN=1, FUN=mean) )

print(x)    # returns => 0.5312099
print(y)    # returns => 0.661242

'apply' 版本实际上比 'for' 版本慢。（根据Inferno作者的说法，如果你这样做，你不是在矢量化，你是在“循环隐藏”。）

但是你可以通过使用built-ins来获得性能提升。下面，我使用内置函数“rowMeans”，对上面两个相同的操作进行了计时：

z = time_this(rowMeans(a))
print(z)    # returns => 0.03679609

与“for”循环（和矢量化版本）相比，改进了一个数量级。

apply 系列的其他成员不仅仅是原生“for”循环的包装器。

a = abs(floor(10*rnorm(1e6)))

time_this(sapply(a, sqrt))
# returns => 6.64 secs

time_this(for (i in 1:length(a)){ sqrt(a[i])})
# returns => 1.33 secs

与“for”循环相比，“sapply”大约慢5 倍。

最后，w / r / t矢量化与'for'循环，如果我可以使用矢量化函数，我认为我永远不会使用循环 - 后者通常是更少的击键并且它是一种更自然的方式（对我来说）编码，我想这是一种不同的性能提升。

于 2010-03-27T02:21:07.510 回答