【问题标题】：R - fastest way to look at all unique pairs of columnsR - 查看所有独特的列对的最快方法
【发布时间】：2013-10-18 22:18:11
【问题描述】：

我有一个数据框 M，我想计算 M 的列之间的所有成对相关性。我可以使用应用函数轻松完成此操作，例如

pvals = laply(M, function(x) llply(M, function(y) cor.test(x, y)$p.value))

但是，由于 x 和 y 之间的相关性与 y 和 x 之间的相关性相同，因此该解决方案所做的工作量是所需工作的 2 倍。

我正在寻找一种快速、简单的方法来计算唯一列对之间的所有相关性。我希望结果是一个 NxN 矩阵，其中 N=ncol(M)。我在 Stack Overflow 上搜索了很长时间，但找不到任何这样做的东西。谢谢！

【问题讨论】：

我考虑过在 i=1:n 和 j=1:i 上使用 apply 语句，但是我无法将所有内容映射回 NxN 矩阵中的正确位置。如果有人知道如何做到这一点，它将解决我的问题。谢谢！
你能做一个小的示例数据集吗，M？
?expand.grid 和 ?combn 并且这些肯定已经被问过，所以你真的应该按照要求进行一些搜索。
我不知道相关性（对角线上的 1）和 cor.test 的输出（对角线上的 0）之间的区别。如果你想要前者，有一种简单的方法可以使用矩阵乘法。
是的！ x = data.frame(复制(10,rnorm(25)));行名（x）= 1:25； colnames(x) = 字母[1:10]。然后，我想得到 a,b,...,j 之间的所有相关性。

标签： r plyr apply correlation lapply

【解决方案1】：

对于虹膜数据，您可以：

data(iris)
r <- cor(iris[1:4])

得到相关矩阵。

您可以查看cor.test 对stats:::cor.test 的实际作用并找到它...

    df <- n - 2L
    ESTIMATE <- c(cor = r)
    PARAMETER <- c(df = df)
    STATISTIC <- c(t = sqrt(df) * r/sqrt(1 - r^2))
    p <- pt(STATISTIC, df)

这都是矢量化的，所以你可以运行它。

在维基百科上有一个关于不同测试的很好的讨论：http://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient

【讨论】：

对于双面测试，您还需要p <- 2 * pmin(p, 1 - p)。
抱歉，有点困惑...这是否计算了 A 与 B 的相关性，以及 B 与 A 的相关性？还是只计算一次相关性？
@ab - cor() 计算变量集的相关矩阵，这确实利用了对称性。其余的只是将其转换为 p 值。
upper.tri 和 lower.tri 对于只获取一组相关性并消除矩阵的另一半很有用。

【解决方案2】：

你可以使用combn:

#Some data:
DF <- USJudgeRatings
#transform to matrix for better subset performance:
m <- as.matrix(DF)
#use combn and its `FUN` argument: 
res <- matrix(nrow=ncol(DF), ncol=ncol(DF))
res[lower.tri(res)] <- combn(seq_along(DF), 2, function(ind) cor.test(m[, ind[[1]]], m[, ind[[2]]])$p.value)
res[upper.tri(res)] <- t(res)[upper.tri(res)]
diag(res) <- 0

基准测试：

corpRoland <- function(DF) {
  m <- as.matrix(DF)
  res <- matrix(nrow=ncol(DF), ncol=ncol(DF))
  res[lower.tri(res)] <- combn(seq_along(DF), 2, function(ind) cor.test(m[, ind[[1]]], m[, ind[[2]]])$p.value)
  res[upper.tri(res)] <- t(res)[upper.tri(res)]
  diag(res) <- 0
  res}

corpNeal <- function(DF) {
  cors <- cor(DF)
  df <- nrow(DF)-2
  STATISTIC <- c(t = sqrt(df) * cors/sqrt(1 - cors^2))
  p <- pt(STATISTIC, df)
  matrix(2 * pmin(p, 1 - p),nrow=ncol(DF))}


library(microbenchmark)
DF <- as.data.frame(matrix(rnorm(1e3), ncol=10))
microbenchmark(corpRoland(DF), corpNeal(DF))
#Unit: microseconds
#           expr       min         lq    median       uq       max neval
# corpRoland(DF) 14021.003 14228.2040 14950.212 15157.27 17013.574   100
#   corpNeal(DF)   342.631   351.6775   373.636   385.34   467.773   100

DF <- as.data.frame(matrix(rnorm(1e4), ncol=100))
microbenchmark(corpRoland(DF), corpNeal(DF), times=10)
# Unit: milliseconds
#           expr         min          lq      median          uq         max neval
# corpRoland(DF) 1595.878487 1601.221980 1615.391891 1633.746678 1637.373231    10
#   corpNeal(DF)    8.359662    8.751755    9.021532    9.509576    9.753154    10

所以，你应该使用@NealFultz 的答案。

【讨论】：

我不知道 diag 应该是一还是零，但我猜你可以将它初始化为初始 matrix(1,... 中的那个。
这取决于替代方案。双面测试为零，见cor.test(1:5,1:5)。