如何计算两个矩阵对应列之间的相关性并仅输出这些相关性，而不是其他相关性。

Question

我有这些数据

> a
     a    b    c
1    1   -1    4
2    2   -2    6
3    3   -3    9
4    4   -4   12
5    5   -5    6

> b
     d    e    f
1    6   -5    7
2    7   -4    4
3    8   -3    3
4    9   -2    3
5   10   -1    9

> cor(a,b)
           d            e             f
a  1.0000000    1.0000000     0.1767767
b -1.0000000    -1.000000    -0.1767767
c  0.5050763    0.5050763    -0.6964286

我想要的结果只是:

cor(a,d) = 1
cor(b,e) = -1
cor(c,f) = -0.6964286

- rder

Answer 1

上面的第一个答案计算了所有成对的相关性，如果矩阵很大，这是可以的，但第二个答案不起作用。据我所知，必须直接进行高效的计算，例如从arrayMagic Bioconductor包借鉴的此代码，适用于大型矩阵：

> colCors = function(x, y) { 
+   sqr = function(x) x*x
+   if(!is.matrix(x)||!is.matrix(y)||any(dim(x)!=dim(y)))
+     stop("Please supply two matrices of equal size.")
+   x   = sweep(x, 2, colMeans(x))
+   y   = sweep(y, 2, colMeans(y))
+   cor = colSums(x*y) /  sqrt(colSums(sqr(x))*colSums(sqr(y)))
+   return(cor)
+ }

> set.seed(1)
> a=matrix(rnorm(15),nrow=5)
> b=matrix(rnorm(15),nrow=5)
> diag(cor(a,b))
[1]  0.2491625 -0.5313192  0.5594564
> mapply(cor,a,b)
 [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
> colCors(a,b)
[1]  0.2491625 -0.5313192  0.5594564

- user1048410

---

是否可以添加p值和多重比较的调整后的p值？ - ferrelwill

Answer 2

我个人会使用diag：

> diag(cor(a,b))
[1]  1.0000000 -1.0000000 -0.6964286

但您也可以使用mapply：

> mapply(cor,a,b)
         a          b          c 
 1.0000000 -1.0000000 -0.6964286

- Joshua Ulrich

Answer 3

mapply 可以用于数据框，但不能用于矩阵。这是因为在数据框中，每一列都是一个元素，而在矩阵中，每个条目都是一个元素。

在上面的答案中，mapply(cor,as.data.frame(a),as.data.frame(b))可以正常工作。

set.seed(1)
a=matrix(rnorm(15),nrow=5)
b=matrix(rnorm(15),nrow=5)
diag(cor(a,b))
[1]  0.2491625 -0.5313192  0.5594564
mapply(cor,as.data.frame(a),as.data.frame(b))
    V1         V2         V3 
 0.2491625 -0.5313192  0.5594564 

这对于大矩阵来说更加高效。

- Cão

Answer 4

对于较大的矩阵，Rfast::corpairs 比 mapply 稍微快一些。

library(Rfast)

a <- matrix(rnorm(1e6), 1e3)
b <- matrix(rnorm(1e6), 1e3)

microbenchmark::microbenchmark(
  diag = diag(cor(a, b)),
  mapply = as.numeric(mapply(cor,as.data.frame(a),as.data.frame(b))),
  corpairs = corpairs(a, b),
  check = "equal"
)
#> Unit: milliseconds
#>      expr       min         lq       mean     median         uq       max neval
#>      diag 1294.1268 1303.72885 1313.41297 1306.84705 1314.44330 1426.1978   100
#>    mapply   49.4884   54.91380   58.52463   57.11995   60.49635  115.6425   100
#>  corpairs    8.6714   12.57265   15.29742   14.81435   16.18880   65.3170   100

- jblood94