如何优化我在 R 中的相关性问题？答案

【问题标题】：How to optimize my correlation problem in R?如何优化我在 R 中的相关性问题？
【发布时间】：2022-01-01 16:13:56
【问题描述】：

我在 R 中有三个数据框，我们称它们为 A、B 和 C。数据框 C 包含两列，第一列包含数据框 A 中的各种行名，第二列包含数据框 B 中的行名：

C <- data.frame(col1 = c("a12", "a9"), col2 = c("b6","b54"))

我想使用表 A 和 B 的行中的相应值计算表 C 的每一行的相关系数和 p 值（即将表 A 中 a12 行中的值与 b6 行中的值相关联来自表 B，来自表 A 的 a9 行和来自表 B 的 b54 行等）并将结果值放在表 C 的其他列中。这是我当前的幼稚且效率极低的代码：

for (i in 1:nrow(C)) { 
correlation <- cor.test(unlist(A[C[i,1],]), unlist(B[C[i,2],]), method = "spearman")
C[i,3] <-correlation$estimate
C[i,4] <- correlation$p.value
}

主要问题是，对于我目前的大型数据集，这种分析实际上可能需要几个月的时间。所以我正在寻找一种更有效的方法来完成这项任务。我还使用“Hmisc”包尝试了以下代码，但我正在处理的服务器无法处理大向量：

A <- t(A)
B <- t(B)
ind.A <- match(C[,1], colnames(A)) 
A<- A[,ind.A]
ind.B <- match(C[,2], colnames(B))
B<- B[,ind.B]
C[,3]<- diag(rcorr(as.matrix(A),as.matrix(B),type = "spearman")$r[c(1:ncol(A)),c(1:ncol(A))])
C[,4]<- diag(rcorr(as.matrix(A),as.matrix(B),type = "spearman")$P[c(1:ncol(A)),c(1:ncol(A))])

【问题讨论】：

加速它的一种方法是使用并行化。看看下面的帖子是否有帮助stackoverflow.com/questions/46532657/…
非常感谢您的指导。这种方法将我的分析速度加快了大约 4 倍

标签： r performance correlation

【解决方案1】：

根据@HYENA 的评论，我尝试了并行处理。这种方法将过程加速了大约 4 倍（使用 8 个内核）。代码：

library(foreach)
library(doParallel)
cl<- makeCluster(detectCores())
registerDoParallel(cl)
cor.res<- foreach (i=1:nrow(C)) %dopar% {
  a<- C[i,1]
  b<- C[i,2]
  correlation<- cor.test(unlist(A[a,]),unlist(B[b,]), method = "spearman")
  c(correlation$estimate,correlation$p.value)
}
cor.res<- data.frame(Reduce("rbind",cor.res))
C[,c(3,4)]<- cor.res

【讨论】：

你的机器上有哪个cpu？
@HYENA 我在自己的笔记本电脑上测试了这种方法，配备 12 Gb RAM 和以下处理器：Intel Corei7-10510U CPU @ 1.80GHz 2.30 GHz 但我正在使用的服务器有 30 Gb RAM以及以下 20 核处理器：Intel Xeon CPU E5-2699 v4 @ 2.20GHz 2.20 GHz
您可以采取其他措施来加快速度： 1、尝试使用 Intel oneAPI MKL 而不是默认库； 2、将线程数设置为略高于物理核心数，以使容量饱和； 3，去除不必要的参数检查代码，如@G。格洛腾迪克的回答者在这里。您可以直接使用Rcpp编写c++代码以更进一步。

【解决方案2】：

仅从 cor.test 中提取您需要的部分并给出 cor_test1 并使用它，或者另外，为给出 cor_test2 的 p 值创建一个查找表，这比 cor_test1 稍快。

基于具有 10 个向量的中值列，它们的运行速度比 cor.test 快约 3 倍。虽然 cor_test2 只比 cor_test1 稍微快一点，但我们将其包括在内，因为速度可能取决于我们没有的输入大小，但您可以自己尝试使用任何大小。

# given correlation and degrees of freedom output p value
r2pval <- function(r, dof) {
  tval <- sqrt(dof) * r/sqrt(1 - r^2)
  min(pt(tval, dof), pt(tval, dof, lower.tail = FALSE))
}

# faster version of cor.test
cor_test1 <- function(x, y) {
  r <- cor(x, y)
  dof <- length(x) - 2
  tval <- sqrt(dof) * r/sqrt(1 - r^2)
  pval <- min(pt(tval, dof), pt(tval, dof, lower.tail = FALSE))
  c(r, pval)
}

# even faster version of cor.test.
# Given x, y and the pvals table calculate a 2-vector of r and p value
cor_test2 <- function(x, y, pvals) {
 r <- cor(x, y)
 c(r, pvals[100 * round(r, 2) + 101])
}

# test
set.seed(123)
n <- 10
x <- rnorm(n); y <- rnorm(n)
dof <- n - 2
# pvals is the 201 p values for r = -1, -0.99, -0.98, ..., 1
pvals <- sapply(seq(-1, 1, 0.01), r2pval, dof = dof)

library(microbenchmark)
microbenchmark(cor.test(x, y), cor_test1(x, y), cor_test2(x, y, pvals))

给予：

Unit: microseconds
                   expr   min    lq    mean median     uq     max neval cld
         cor.test(x, y) 253.7 256.7 346.278 266.05 501.45   650.6   100   a
        cor_test1(x, y)  84.8  87.2 346.777  89.10 107.40 22974.4   100   a
 cor_test2(x, y, pvals)  72.4  75.0 272.030  79.45  91.25 17935.8   100   a

【讨论】：