¿Hay una función R que aplica una función a cada par de columnas?

Question

A menudo necesito aplicar una función a cada par de columnas en un marco de datos/matriz y devolver los resultados en una matriz. Ahora siempre escribo un ciclo para hacer esto. Por ejemplo, para hacer una matriz que contiene los valores de p de correlaciones que escribo:

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 

que funciona, pero es bastante lento para matrices muy grandes. Puedo escribir una función para este en R (sin molestarse con el corte de la mitad el tiempo asumiendo un resultado simétrico al anterior):

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

o una función con 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") 

Pero ambos son bastante lento incluso en un muy pequeño conjunto de datos de 100 variables (pensé que la función RCPP sería más rápido, pero supongo que la conversión entre R y C++ todo el tiempo se cobra su peaje):

> 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 

Así que mi pregunta es:

1. Debido a la simplicidad de estas funciones, supongo que esto ya está en algún lugar de R. ¿Existe una función apply o plyr que hace esto? Lo he buscado pero no he podido encontrarlo.
2. Si es así, ¿es más rápido?

Answer 1

No sería más rápido, pero se puede utilizar outer para simplificar el código. Requiere una función vectorizada, así que aquí he usado Vectorize para hacer una versión vectorizada de la función para obtener la correlación entre dos columnas.

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

No estoy seguro de si esto resuelve su problema de manera adecuada, pero eche un vistazo al paquete psych de William Revelle. corr.test devuelve una lista de matrices con coeficientes de correlación, # de obs, estadística de prueba t y valor p. Sé que lo uso todo el tiempo (y AFAICS también eres un psicólogo, por lo que puede satisfacer tus necesidades también). Escribir bucles no es la forma más elegante de hacerlo.

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% del tiempo que se gasta en cor.test.default y rutinas que llama por lo que su desesperada tratando de obtener resultados más rápidos simplemente reescribiendo Papply (aparte de los ahorros de la computación sólo aquellos por encima o por debajo de la diagonal asumiendo que su la función es simétrica en x y 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

Puede utilizar mapply, pero como dicen las otras respuestas es poco probable que sea mucho más rápido ya que la mayoría de las veces está siendo utilizado por 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) 

Se podría reducir la cantidad de trabajo mapply hace mediante el uso de la suposición de simetría y tomando nota de la diagonal cero, por ejemplo

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