是否有一個R函數將函數應用於每對列？

Question

我經常需要對數據框/矩陣中的每對列應用函數，並將結果以矩陣形式返回。現在我總是寫一個循環來做到這一點。例如，爲了使含我寫相關的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 

其作品，但對於非常大的矩陣相當緩慢。

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") 

但兩者都相當我可以在R（不與假設如上對稱的結果切削時間縮短了一半打擾）寫一個函數爲這個減緩甚至在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 

所以我的問題是：

1. 由於這些函數的簡單性，我認爲這已經在R的某個地方了。是否有應用程序或plyr函數執行此操作？我一直在尋找它，但一直沒能找到它。
2. 如果是這樣，它是否更快？

Answer 1

它不會更快，但您可以使用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) 

Answer 2

我不確定這是否以正確的方式解決您的問題，但看看William Revelle的psych包。 corr.test返回具有相關係數，obs數，t檢驗統計量和p值的矩陣列表。我知道我一直都在使用它（而AFAICS你也是一名心理學家，所以它也可以滿足你的需求）。編寫循環並不是這樣做的最優雅的方式。

library(psych) 
corr.test(mtcars) 
(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" 

Answer 3

的時間

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 ... 

Answer 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