是否有将函数应用于每对列的 R 函数？

Posted 2023-02-14

【中文标题】是否有将函数应用于每对列的 R 函数？

【英文标题】：Is there a R function that applies a function to each pair of columns?

【发布时间】：2011-07-11 03:17:59

【问题描述】：

我经常需要对数据框/矩阵中的每一对列应用一个函数，并以矩阵形式返回结果。现在我总是写一个循环来做到这一点。例如，要制作一个包含相关 p 值的矩阵，我会写：

df <- data.frame(x=rnorm(100),y=rnorm(100),z=rnorm(100))

n <- ncol(df)

foo <- matrix(0,n,n)

for ( i in 1:n)

    for (j in i:n)
    
        foo[i,j] <- cor.test(df[,i],df[,j])$p.value
    


foo[lower.tri(foo)] <- t(foo)[lower.tri(foo)]

foo
          [,1]      [,2]      [,3]
[1,] 0.0000000 0.7215071 0.5651266
[2,] 0.7215071 0.0000000 0.9019746
[3,] 0.5651266 0.9019746 0.0000000

这可行，但对于非常大的矩阵来说非常慢。我可以在 R 中为此编写一个函数（不必担心通过假设上述对称结果将时间减半）：

Papply <- function(x,fun)

n <- ncol(x)

foo <- matrix(0,n,n)
for ( i in 1:n)

    for (j in 1:n)
    
        foo[i,j] <- fun(x[,i],x[,j])
    

return(foo)

或带有 Rcpp 的函数：

library("Rcpp")
library("inline")

src <- 
'
NumericMatrix x(xR);
Function f(fun);
NumericMatrix y(x.ncol(),x.ncol());

for (int i = 0; i < x.ncol(); i++)

    for (int j = 0; j < x.ncol(); j++)
    
        y(i,j) = as<double>(f(wrap(x(_,i)),wrap(x(_,j))));
    

return wrap(y);
'

Papply2 <- cxxfunction(signature(xR="numeric",fun="function"),src,plugin="Rcpp")

但即使在包含 100 个变量的非常小的数据集上，两者都相当慢（我认为 Rcpp 函数会更快，但我猜 R 和 C++ 之间的转换一直都会造成损失）：

> system.time(Papply(matrix(rnorm(100*300),300,100),function(x,y)cor.test(x,y)$p.value))
   user  system elapsed 
   3.73    0.00    3.73 
> system.time(Papply2(matrix(rnorm(100*300),300,100),function(x,y)cor.test(x,y)$p.value))
   user  system elapsed 
   3.71    0.02    3.75 

所以我的问题是：

由于这些函数的简单性，我假设这已经在 R 中的某个地方。是否有执行此操作的 apply 或 plyr 函数？我已经找过了，但没找到。 如果是这样，是不是更快？

【参考方案1】：

我不确定这是否能以适当的方式解决您的问题，但请查看 William Revelle 的 psych 包。 corr.test 返回具有相关系数、# of obs、t 检验统计量和 p 值的矩阵列表。我知道我一直都在使用它（而且 AFAICS 你也是一名心理学家，所以它也可以满足你的需求）。编写循环并不是最优雅的方式。

> library(psych)
> ( k <- corr.test(mtcars[1:5]) )
Call:corr.test(x = mtcars[1:5])
Correlation matrix 
       mpg   cyl  disp    hp  drat
mpg   1.00 -0.85 -0.85 -0.78  0.68
cyl  -0.85  1.00  0.90  0.83 -0.70
disp -0.85  0.90  1.00  0.79 -0.71
hp   -0.78  0.83  0.79  1.00 -0.45
drat  0.68 -0.70 -0.71 -0.45  1.00
Sample Size 
     mpg cyl disp hp drat
mpg   32  32   32 32   32
cyl   32  32   32 32   32
disp  32  32   32 32   32
hp    32  32   32 32   32
drat  32  32   32 32   32
Probability value 
     mpg cyl disp   hp drat
mpg    0   0    0 0.00 0.00
cyl    0   0    0 0.00 0.00
disp   0   0    0 0.00 0.00
hp     0   0    0 0.00 0.01
drat   0   0    0 0.01 0.00

> str(k)
List of 5
 $ r   : num [1:5, 1:5] 1 -0.852 -0.848 -0.776 0.681 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
  .. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
 $ n   : num [1:5, 1:5] 32 32 32 32 32 32 32 32 32 32 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
  .. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
 $ t   : num [1:5, 1:5] Inf -8.92 -8.75 -6.74 5.1 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
  .. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
 $ p   : num [1:5, 1:5] 0.00 6.11e-10 9.38e-10 1.79e-07 1.78e-05 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
  .. ..$ : chr [1:5] "mpg" "cyl" "disp" "hp" ...
 $ Call: language corr.test(x = mtcars[1:5])
 - attr(*, "class")= chr [1:2] "psych" "corr.test"

【讨论】：

很好，谢谢！相关性 p 值只是我今天碰巧遇到的一个例子。

【参考方案2】：

92% 的时间都花在了cor.test.default 和它调用的例程上，因此通过简单地重写Papply 来获得更快的结果是毫无希望的（除了只计算对角线之上或之下的那些，假设你的函数在x 和y 中是对称的）。

> M <- matrix(rnorm(100*300),300,100)
> Rprof(); junk <- Papply(M,function(x,y) cor.test( x, y)$p.value); Rprof(NULL)
> summaryRprof()
$by.self
                 self.time self.pct total.time total.pct
cor.test.default      4.36    29.54      13.56     91.87
# ... snip ...

【参考方案3】：

不会更快，但您可以使用outer 来简化代码。它确实需要一个向量化函数，所以这里我使用Vectorize 来制作函数的向量化版本来获取两列之间的相关性。

df <- data.frame(x=rnorm(100),y=rnorm(100),z=rnorm(100))
n <- ncol(df)

corpij <- function(i,j,data) cor.test(data[,i],data[,j])$p.value
corp <- Vectorize(corpij, vectorize.args=list("i","j"))
outer(1:n,1:n,corp,data=df)

【讨论】：

如果我的函数有 data1 和 data2 并且只有 i，这是否正确？ corpij <- function(i,data1,data2) cor.test(data1[,i],data2[i])$p.value corp <- Vectorize(corpij, vectorize.args=list("i")) outer(1:n,1:n,corp,data1=df,data2=df2) 但是，这给了我错误：Error in FUN(X, Y, ...) : unused argument (Y)

【参考方案4】：

您可以使用mapply，但正如其他答案所述，它不太可能更快，因为大部分时间都被cor.test 用完了。

matrix(mapply(function(x,y) cor.test(df[,x],df[,y])$p.value,rep(1:3,3),sort(rep(1:3,3))),nrow=3,ncol=3)

您可以通过使用对称假设并注意零对角线来减少 mapply 所做的工作量，例如

v <- mapply(function(x,y) cor.test(df[,x],df[,y])$p.value,rep(1:2,2:1),rev(rep(3:2,2:1)))
m <- matrix(0,nrow=3,ncol=3)
m[lower.tri(m)] <- v
m[upper.tri(m)] <- v