公认的解决方案(心理包中的 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