两个矩阵的所有行之间的相关性/ p 值组合

Question

我希望能够计算每个物种(bac)与第二个数据框中的每个因子(fac)之间的相关性和p值。两者在相同数量的站点上进行了测量，但bac和fac的数量不匹配。

bac1 <- c(1,2,3,4,5)
bac2 <- c(2,3,4,5,1)
bac3 <- c(4,5,1,2,3)
bac4 <- c(5,1,2,3,4)
bac <- as.data.frame(cbind(bac1, bac2, bac3, bac4 ))
colnames(bac) <- c("station1", "station2", "station3", "station4")
rownames(bac) <- c("bac1", "bac2", "bac3", "bac4", "bac5")

fac1 <- c(1,2,3,4,5,6)
fac2 <- c(2,3,4,5,1,6)
fac3<- c(3,4,5,1,2,6)
fac4<- c(4,5,1,2,3, 6)
fac <- as.data.frame(cbind(fac1, fac2, fac3, fac4))
colnames(fac) <- c("station1", "station2", "station3", "station4")
rownames(fac) <- c("fac1", "fac2", "fac3", "fac4", "fac5", "fac6")

我想象结果看起来有点像这样，同时保留名称以知道呈现的组合是哪个：

bac1-fac1 cor1 p1
bac1-fac2 cor2 p2
bac1-fac3 cor3 p3

bac2-fac1 corx px...

我查看了Hmist的rcorr函数和psych的corr.test函数，但是找不到必要行重排的示例... 有什么想法吗？

- Helena

Answer 1

如果您重构数据，使其计算成对列之间的相关性，那将非常容易。

tbac <- data.frame(t(bac))
tfac <- data.frame(t(fac))

f <- function (x, y) cor(x, y)

tab <- outer(tfac, tbac, Vectorize(f))

as.data.frame.table(tab)

我使用同样的想法得出了一个答案：匹配数据并计算相同值的数量。

- Zheyuan Li

---

非常紧凑的方式来完成这个。一如既往的好答案！ - akrun

我不记得了，但感谢分享。 - akrun

谢谢，这似乎非常有效！我想知道为什么与fac6的任何相关性都会产生NA，但是我已经弄清楚了（所有值都是6）。 - Helena

Answer 2

您可以直接将完整的矩阵传递给cor函数（或psych::corr.test），它会负责找到相关列的相关性。例如：

cor(t(fac), t(bac))
#            bac1        bac2        bac3        bac4        bac5
# fac1  0.9899495 -0.07559289 -0.60000000 -0.60000000 -0.07559289
# fac2  0.9899495 -0.07559289 -0.60000000 -0.60000000 -0.07559289
# fac3 -0.3207135  0.94285714 -0.07559289 -0.07559289 -0.48571429
# fac4 -0.8000000 -0.32071349  0.98994949  0.98994949 -0.32071349
# fac5 -0.3207135 -0.48571429 -0.07559289 -0.07559289  0.94285714
# fac6         NA          NA          NA          NA          NA

您可以使用reshape2::melt将其转换为长格式。

reshape2::melt(cor(t(fac), t(bac)))
#    Var1 Var2       value
# 1  fac1 bac1  0.98994949
# 2  fac2 bac1  0.98994949
# 3  fac3 bac1 -0.32071349
# 4  fac4 bac1 -0.80000000
# ---
# ---

使用相同的方法获取p值。

test <- psych::corr.test(t(fac), t(bac), adjust="none")

然后像之前一样融化并连接在一起

merge(melt(test$r, value.name="cor"), melt(test$p, value.name="p-value"), by=c("Var1", "Var2"))
#   Var1 Var2         cor    p-value
# 1 fac1 bac1  0.98994949 0.01005051
# 2 fac1 bac2 -0.07559289 0.92440711
# 3 fac1 bac3 -0.60000000 0.40000000
# 4 fac1 bac4 -0.60000000 0.40000000
# 5 fac1 bac5 -0.07559289 0.92440711
# 6 fac2 bac1  0.98994949 0.01005051

- user20650

---

那是一个很好的选择。我错过了转置部分。 - akrun

Answer 3

我们可以使用expand.grid来获取'bac'和'fac'的rownames的组合，通过指定MARGIN为1，循环遍历行，并基于行名对'bac'和'fac'的行进行子集划分，进行corr.test并提取'p'值作为list。

library(psych)
do.call(c, apply(expand.grid(rownames(bac), rownames(fac)), 1, 
  function(x) list(corr.test(cbind(unlist(bac[1,]), unlist(fac[1,])))$p)))

- akrun

---

@李哲源ZheyuanLi 它还会将其他参数以list的形式返回。我认为Hmisc中的rcorr也会执行类似的操作。 - akrun

我最近发现了expand.grid，我非常喜欢它。但是当我尝试你的解决方案时，输出似乎不正确...而且我没有任何行/列名称？ - Helena

Answer 4

你只需要循环遍历expand.grid的行

pairs <- as.matrix(expand.grid(1:nrow(bac),1:nrow(fac)))
pairs <- cbind(pairs,NA,NA)
b <- as.matrix(bac)
f <- as.matrix(fac)
for(i in 1:nrow(pairs)){
    pairs[i,3] <- cor(b[pairs[i,1],], f[pairs[i,2],])
    pairs[i,4] <- cor.test(b[pairs[i,1],], f[pairs[i,2],])$p.value
}
colnames(pairs) <- c('bac','fac','corr','p')
pairs
##      bac fac        corr          p
## [1,]   1   1  0.98994949 0.01005051
## [2,]   2   1 -0.07559289 0.92440711
## [3,]   3   1 -0.60000000 0.40000000
## [4,]   4   1 -0.60000000 0.40000000
## [5,]   5   1 -0.07559289 0.92440711
## [6,]   1   2  0.98994949 0.01005051

如果你想要名称，那么可以这样做：

pairs <- as.data.frame(pairs)
pairs[,1] <- sapply(pairs[,1],function(x) rownames(bac)[x])
pairs[,2] <- sapply(pairs[,2],function(x) rownames(fac)[x])

虽然在那个时候，使用李哲源 Zheyuan Li 的解决方案可能更容易。

- Ryan

---

谢谢，也很有帮助，但是没有保留原始名称，在我的“真实”情况下会很有帮助！ - Helena