r - 计算多个变量的多个组之间的倍数变化的快速而优雅的方法？

Question

出于好奇，我一直在玩几种计算折叠变化的方法，我试图找到最快和最优雅的方法来做到这一点(希望这也是相同的解决方案)。

我感兴趣的矩阵类型如下所示:

# Some data
nvars <- 10000
nsamples <- 500
sample_groups <- 5
MAT <- replicate(nvars, runif(n=nsamples))

和一个看起来像这样的分组向量:

f <- rep_len(1:sample_groups, nsamples)
f <- LETTERS[f]

我最终希望从上述输入中得到一个 10 x 10,000 的矩阵，其中 f 中的每个级别组合都有一行。

为此，第一个任务是计算所有列的每个组的平均值。
我想出了 4 种可能的方法来做到这一点:

# Settings
aggr_FUN  <- mean
combi_FUN <- function(x,y) "/"(x,y) 

# helper function
pasteC <- function(x,y) paste(x,y,sep=" - ")

#A1. Loop
system.time({
f_un <- unique(f)
temp1 <- matrix(NA,nrow = length(f_un),ncol=ncol(MAT))
rownames(temp1) <- f_un

for(i in 1:length(f_un)){
  temp1[i,] <- apply(MAT[f_un[i] == f,,drop=FALSE],2,aggr_FUN)
}
})


 user  system elapsed 
   0.41    0.00    0.41 

#A2. aggregate
system.time({
temp2 <- aggregate(. ~ class, data = cbind.data.frame(class=f,MAT), aggr_FUN)
})


 user  system elapsed 
   7.76    0.05    7.81 

#A3. reshape2
library(reshape2)
system.time({
temp3 <- recast(data.frame(class=f,MAT),class ~ variable,id.var="class",aggr_FUN)
})


  user  system elapsed 
   1.82    0.30    2.12 

#A4. purrr
library(purrr)
system.time({
temp4 <- data.frame(class = f, MAT) %>%
  slice_rows("class") %>%
  by_slice(map, aggr_FUN)
})


   user  system elapsed 
   0.47    0.00    0.47 

正如你所看到的，循环实际上是最快的解决方案， purrr 包只是稍微慢了一点。 recast 慢 5 倍， aggregate 明显更宽松。
我也尝试过 dplyr 包，但由于某种原因( https://github.com/hadley/dplyr/issues/1395 )，结果非常慢。
purrr 既快速又优雅，因此对于这部分来说，这是一个寻找更好方法的学术练习。

此时我们的输出是:

> temp1[,1:6]
       [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
A 0.4804964 0.4779168 0.5292458 0.4401357 0.4728515 0.5009800
B 0.4819612 0.5260592 0.5291887 0.5095620 0.4840777 0.4792213
C 0.4661714 0.4886010 0.5006018 0.5061170 0.5058892 0.5432819
D 0.4566942 0.4519988 0.5334207 0.4912822 0.4542889 0.4898384
E 0.4967948 0.5630683 0.4941777 0.5239327 0.5045152 0.5227140

因此，如果您仍在阅读这里，那么更具挑战性的部分。我们需要计算所有组/行组合之间的倍数变化。

我找到了两种方法来做到这一点:

#B1. by loop
combs <- t(combn(as.character(f_un),2))
combi_FUN_vec <- Vectorize(combi_FUN)

out <- matrix(NA,nrow = nrow(combs),ncol=ncol(temp1))
rownames(out) <- pasteC(combs[,1],combs[,2])
colnames(out) <- 1:ncol(temp1)

system.time({

for( i in 1:nrow(combs)){
  out[i,] <- combi_FUN_vec(    temp1[combs[i,1],]      ,         temp1[combs[i,2],]   )
}

})


  user  system elapsed 
   0.13    0.00    0.13

#B2. by apply
class_computed <- as.character(temp2[,1])
temp2 <- as.matrix(temp2[,-1])
combs <- t(combn(class_computed,2))
rownames(temp2) <- class_computed

combi_FUN_vec <- Vectorize(combi_FUN)

system.time({

out <- apply(temp2,2,function(x){
  v <- combi_FUN_vec(    x[combs[,1]]      ,        x[combs[,2]]    )
  names(v) <- pasteC(combs[,1],combs[,2])
  return(v)
})

})


   user  system elapsed 
   0.91    0.00    0.91 

毫不奇怪，循环显然是赢家，输出是这样的:

> out[,1:5]
              1         2         3         4         5
A - B 1.2128952 1.0161608 0.9313115 0.9767619 1.0132362
A - C 1.0946079 1.0524154 0.9882857 0.9546686 0.9604382
A - D 1.1872958 0.9113349 0.9941437 0.8751611 0.9863873
A - E 1.1457396 0.9669100 0.9146375 0.8760513 1.0604971
B - C 0.9024753 1.0356780 1.0611763 0.9773810 0.9478918
B - D 0.9788940 0.8968413 1.0674664 0.8959820 0.9735018
B - E 0.9446320 0.9515325 0.9820962 0.8968933 1.0466435
C - D 1.0846768 0.8659461 1.0059275 0.9167172 1.0270179
C - E 1.0467123 0.9187532 0.9254788 0.9176496 1.1041804
D - E 0.9649993 1.0609821 0.9200254 1.0010171 1.0751326

现在这就是我的问题......最后两个方法非常难看。
有没有更好/更清洁/更快的方法？最好使用 dplyr/ purrr 风格的语法？也许甚至不必通过 combn ？

任何提示表示赞赏。

编辑:

我设法以 dplyr 风格制作了一个更紧凑的版本:

f_un <- unique(f)
  combs <- t(combn(as.character(f_un),2))

  out3 <- data.frame(class = f, MAT) %>%  slice_rows("class") %>% by_slice(map, aggr_FUN) %>% 
          do(combi_FUN( slice(.,match(combs[,1], class))[,-1]  ,slice(.,match(combs[,2], class))[,-1]      )) %>% 
          as.data.frame(row.names = pasteC(combs[,1],combs[,2]))

有没有办法简化并加快速度？它比上面最快的慢 10 倍。

编辑2:
根据迄今为止的建议，最快和最干净的是以下功能。

fold.change <- function(MAT,f,aggr_FUN=mean,combi_FUN=function(x,y) "/"(x,y)   ){

  # mean using purrr
  x <- data.frame(class = f, MAT) %>%  slice_rows("class") %>% by_slice(map, aggr_FUN)
  rownames <- as.character(as.data.frame(x[,1])[,1])
  x <- as.matrix(x[,-1])
  rownames(x) <- rownames

  # calculate changes between all rows
  i <- combn(unique(f), 2)
  ret <- combi_FUN(x[i[1,],] , x[i[2,],])
  rownames(ret) <- pasteC(i[1,], i[2,])

  # Put original colnames
  colnames(ret) <- colnames(MAT)

  return(ret)
}

Answer 1

矩阵运算和子集运算速度很快:

fold <- function(x, f, aggr_FUN = colMeans, combi_FUN = '/'){
  f <- as.factor(f)
  i <- split(1:nrow(x), f)
  x <- sapply(i, function(i){ aggr_FUN(x[i,])})
  x <- t(x)
  j <- combn(levels(f), 2)
  ret <- combi_FUN(x[j[1,],], x[j[2,],])
  rownames(ret) <- paste(j[1,], j[2,], sep = '-')
  ret
}

> system.time(ret <- fold(MAT, f))
   user  system elapsed 
   0.13    0.00    0.12 
> all.equal(ret, out, check.attributes = F)
[1] TRUE
> if(require(matrixStats))  
+   system.time(fold(MAT, f, aggr_FUN = colMedians))
   user  system elapsed 
   0.27    0.00    0.27 
> if(require(matrixStats))
+   system.time(fold(MAT, f, aggr_FUN = colSds))
   user  system elapsed 
   0.17    0.02    0.18 

除非，我真的误解了你想要做什么。