cor.test() 的矩阵版本答案

【问题标题】：A matrix version of cor.test()cor.test() 的矩阵版本
【发布时间】：2012-10-18 05:07:02
【问题描述】：

Cor.test() 将向量 x 和 y 作为参数，但我有一个完整的数据矩阵要成对测试。 Cor() 将此矩阵作为参数就好了，我希望找到一种方法来为 cor.test() 做同样的事情。

其他人的常见建议似乎是使用cor.prob()：

https://stat.ethz.ch/pipermail/r-help/2001-November/016201.html

但是这些p值和cor.test()生成的p值不一样！！！ Cor.test() 似乎也比 cor.prob() 更好地处理成对删除（我的数据集中有相当多的缺失数据）。

有没有人可以替代cor.prob()？如果解决方案涉及嵌套 for 循环，那就这样吧（我对 R 足够新，即使这对我来说也是个问题）。

【问题讨论】：

您可以将lapply 与cor.test 一起使用，或者将函数矢量化并将其提供给outer，如以下链接所示：stackoverflow.com/questions/9917242/…

标签： r correlation

【解决方案1】：

如果您严格追求来自cor.test 的矩阵格式的 pvalue，这里有一个从文森特 (LINK) 处无耻窃取的解决方案：

cor.test.p <- function(x){
    FUN <- function(x, y) cor.test(x, y)[["p.value"]]
    z <- outer(
      colnames(x), 
      colnames(x), 
      Vectorize(function(i,j) FUN(x[,i], x[,j]))
    )
    dimnames(z) <- list(colnames(x), colnames(x))
    z
}

cor.test.p(mtcars)

注意：Tommy 还提供了一个更快的解决方案，但不太容易实现。哦，没有 for 循环 :)

编辑我的 qdapTools 包中有一个函数 v_outer，它使这项任务变得非常简单：

library(qdapTools)
(out <- v_outer(mtcars, function(x, y) cor.test(x, y)[["p.value"]]))
print(out, digits=4)  # for more digits

【讨论】：

编辑和[[3]] 索引cor.test 输出的列表。该列表的第三个元素是 p.value。
@TylerRinker 我发现如果使用列表输出的命名版本，代码会更清晰。如果不是 cor.test(x, y)[[3]] 你有 cor.test(x, y)[["p.value"]] 你从测试中提取 p 值，那就更清楚了。
@Dason 我同意我只是懒惰，因为我猜测索引是基于输出的，并且太懒了，在 cor.test 的输出上使用 str 或 names找出。我真的为此责怪机器人。他们已经自动化了我们的生活，以至于我们太懒了。
你是说你的提案可以达到和p.mat.all <- psych:::cor.test(M.cor, alternative = "two.sided", method = c("pearson", "kendall", "spearman"), adjust = "none", ci = F)一样的结果？ - - 我想你只是在这里使用 Pearson cor。
我喜欢这种方法，谢谢！我需要计算多个成对相关的 p-val，并且 rcorr 没有在我的数据中运行，因为它是由非常大的向量组成的。这成功了！谢谢！！

【解决方案2】：

psych 包中的corr.test 旨在执行此操作：

library("psych")
data(sat.act)
corr.test(sat.act)

如 cmets 中所述，要在整个矩阵中复制基础 cor.test() 函数中的 p 值，则需要关闭 p-values 进行多重比较（默认是使用 Holm 的调整方法）：

 corr.test(sat.act, adjust = "none")

[但在解释这些结果时要小心！]

【讨论】：

漂亮，为什么要重新发明轮子。 +1g
如果您希望结果与统计数据匹配，请备注cor.test 使用corr.test(mtcars, adjust="none")
泰勒，我注意到了。谢谢！你们俩都很棒而且非常乐于助人！
如果你有一个大矩阵，这将非常非常慢！要加快速度，请设置参数 ci=F - 运行时间大约是 cor() 的两倍，而使用 ci=T（默认值）可能需要 100 倍的时间。
当我尝试执行 "ci = F ”。我在下面写了一个答案，它从函数中获取重要代码，然后运行 cor() 并给出 pvalues。

【解决方案3】：

可能最简单的方法是使用来自 Hmisc 的rcorr()。它只需要一个矩阵，所以如果您的数据在 data.frame 中，请使用rcorr(as.matrix(x))。它将返回一个列表，其中包含：1) r 成对矩阵，2) 成对 n 矩阵，3) r 的 p 值矩阵。它会自动忽略丢失的数据。

理想情况下，此类函数也应采用 data.frames 并输出符合“New Statistics”的置信区间。

【讨论】：

这是理想的，但它没有在我的大型数据集（50 个变量（我正在评估它们的相似性）x 46,000,000 个观察值）上运行。出现内存错误。
试试 weights 包中的wtd.cors()。我认为它使用某种快速的近似值。如果您需要 p 值等，您可以在每个成对变量上使用 wtd.cor()。如果您仍然想要更快的速度，您可以考虑一次执行一个变量并保存计算之间的 z 分数，因为这样可以省去多次重新计算它们的操作。

【解决方案4】：

公认的解决方案（psych 包中的 corr.test 函数）有效，但对于大型矩阵来说非常慢。我正在使用与药物敏感性矩阵（~1,000 x ~500）相关的基因表达矩阵（~20,000 x ~1,000），我不得不停止它，因为它需要很长时间。

我从 psych 包中提取了一些代码，直接使用了 cor() 函数，得到了更好的结果：

# find (pairwise complete) correlation matrix between two matrices x and y
# compare to corr.test(x, y, adjust = "none")
n <- t(!is.na(x)) %*% (!is.na(y)) # same as count.pairwise(x,y) from psych package
r <- cor(x, y, use = "pairwise.complete.obs") # MUCH MUCH faster than corr.test()
cor2pvalue = function(r, n) {
  t <- (r*sqrt(n-2))/sqrt(1-r^2)
  p <- 2*(1 - pt(abs(t),(n-2)))
  se <- sqrt((1-r*r)/(n-2))
  out <- list(r, n, t, p, se)
  names(out) <- c("r", "n", "t", "p", "se")
  return(out)
}
# get a list with matrices of correlation, pvalues, standard error, etc.
result = cor2pvalue(r,n)

即使有两个 100 x 200 矩阵，差异也是惊人的。一两秒对 45 秒。

> system.time(test_func(x,y))
   user  system elapsed 
  0.308   2.452   0.130 
> system.time(corr.test(x,y, adjust = "none"))
   user  system elapsed 
 45.004   3.276  45.814

【讨论】：

注意：我刚刚在上面看到您可以使用带有“ci = F”选项的 corr.test() 来加快速度。但是，当我尝试它时它给了我一个错误。
看起来代码中有一个小错误。在这里查看我的修复（我知道它是只读的）：github.com/cran/psych/pull/2/commits/… 我通过电子邮件向包维护者发送了邮件。

【解决方案5】：

“公认的解决方案（心理包中的corr.test 函数）有效，但对于大型矩阵来说非常慢。”

如果你使用ci=FALSE，那么速度会快很多。默认情况下，会找到置信区间。但是，这会导致速度略有下降。所以，对于rs、ts 和ps，设置ci=FALSE。

【讨论】：