test - r correlation function

Una versión matricial de cor.test()

Question

Cor.test() toma vectores x e y como argumentos, pero tengo una matriz completa de datos que quiero probar, por pares. Cor() toma esta matriz como un argumento muy bien, y espero encontrar una manera de hacer lo mismo para cor.test() .

El consejo común de otras personas parece ser usar cor.prob() :

https://stat.ethz.ch/pipermail/r-help/2001-November/016201.html

Pero estos valores de p no son los mismos que los generados por cor.test() !!! Cor.test() también parece estar mejor equipado para manejar la eliminación por pares (tengo un poco de datos faltantes en mi conjunto de datos) que cor.prob() .

¿Alguien tiene alguna alternativa a cor.prob() ? Si la solución involucra anidados para bucles, que así sea (soy lo suficientemente nuevo para R incluso para que esto sea problemático para mí).

Answer 1

"La solución aceptada (función corr.test en el paquete psych) funciona, pero es extremadamente lenta para matrices grandes".

Si usas ci=FALSE , entonces la velocidad es mucho más rápida. Por defecto, se encuentran intervalos de confianza. Sin embargo, esto conduce a una ligera desaceleración de la velocidad. Entonces, solo para rs , ts y ps , establezca ci=FALSE .

Answer 2

La solución aceptada (función corr.test en el paquete psych) funciona, pero es extremadamente lenta para matrices grandes. Estaba trabajando con una matriz de expresión génica (~ 20,000 por ~ 1,000) correlacionada con una matriz de sensibilidad al fármaco (~ 1,000 por ~ 500) y tuve que detenerla porque tardaba una eternidad.

Tomé algo de código del paquete psych y usé la función cor () directamente en su lugar y obtuve resultados mucho mejores:

# find (pairwise complete) correlation matrix between two matrices x and y
# compare to corr.test(x, y, adjust = "none")
n <- t(!is.na(x)) %*% (!is.na(y)) # same as count.pairwise(x,y) from psych package
r <- cor(x, y, use = "pairwise.complete.obs") # MUCH MUCH faster than corr.test()
cor2pvalue = function(r, n) {
  t <- (r*sqrt(n-2))/sqrt(1-r^2)
  p <- 2*(1 - pt(abs(t),(n-2)))
  se <- sqrt((1-r*r)/(n-2))
  out <- list(r, n, t, p, se)
  names(out) <- c("r", "n", "t", "p", "se")
  return(out)
}
# get a list with matrices of correlation, pvalues, standard error, etc.
result = cor2pvalue(r,n)

Incluso con dos matrices de 100 x 200, la diferencia fue asombrosa. Un segundo o dos contra 45 segundos.

> system.time(test_func(x,y))
   user  system elapsed 
  0.308   2.452   0.130 
> system.time(corr.test(x,y, adjust = "none"))
   user  system elapsed 
 45.004   3.276  45.814 

Answer 3

Probablemente la forma más fácil es usar el rcorr() de Hmisc. Solo tomará una matriz, así que use rcorr(as.matrix(x)) si sus datos están en un data.frame. Le devolverá una lista con: 1) matriz de r por pares, 2) matriz de por pares n, 3) matriz de p valores para los r. Ignora automáticamente los datos faltantes.

Idealmente, una función de este tipo debería tomar también data.frames y también generar intervalos de confianza en línea con las '' Nuevas estadísticas ''.

Answer 4

Si está siguiendo estrictamente los valores de p en un formato de matriz de cor.test aquí hay una solución descaradamente robada de Vincent ( LINK ):

cor.test.p <- function(x){
    FUN <- function(x, y) cor.test(x, y)[["p.value"]]
    z <- outer(
      colnames(x), 
      colnames(x), 
      Vectorize(function(i,j) FUN(x[,i], x[,j]))
    )
    dimnames(z) <- list(colnames(x), colnames(x))
    z
}

cor.test.p(mtcars)

Nota: Tommy también proporciona una solución más rápida, aunque menos fácil de implementar. Ah y no para loops :)

Editar Tengo una función v_outer en mi paquete qdapTools que hace esta tarea bastante fácil:

library(qdapTools)
(out <- v_outer(mtcars, function(x, y) cor.test(x, y)[["p.value"]]))
print(out, digits=4)  # for more digits

Answer 5

corr.test en el paquete psych está diseñado para hacer esto:

library("psych")
data(sat.act)
corr.test(sat.act)

Como se señaló en los comentarios, para replicar los valores p de la función cor.test() base en toda la matriz, debe desactivar el ajuste de los valores p para comparaciones múltiples (el valor predeterminado es usar el método de Holm para ajuste):

 corr.test(sat.act, adjust = "none")

[¡Pero ten cuidado al interpretar esos resultados!]