¿Cuál es la forma más rápida de aplicar t.test a cada columna de una matriz grande?

Question

Supongamos que tengo una matriz grande:

M <- matrix(rnorm(1e7),nrow=20) 

Supongamos además que cada columna representa una muestra. Digamos que me gustaría aplicar t.test() a cada columna, ¿hay alguna manera de hacerlo que sea mucho más rápido que usar apply()?

apply(M, 2, t.test) 

Se tomó un poco menos de 2 minutos para ejecutar el análisis en mi equipo:

> system.time(invisible(apply(M, 2, t.test))) 
user system elapsed 
113.513 0.663 113.519 

Answer 1

Si usted tiene una máquina de múltiples núcleos que hay algunos beneficios del uso de todos los núcleos, por ejemplo usando mclapply.

> library(multicore) 
> M <- matrix(rnorm(40),nrow=20) 
> x1 <- apply(M, 2, t.test) 
> x2 <- mclapply(1:dim(M)[2], function(i) t.test(M[,i])) 
> all.equal(x1, x2) 
[1] "Component 1: Component 9: 1 string mismatch" "Component 2: Component 9: 1 string mismatch" 
# str(x1) and str(x2) show that the difference is immaterial 

Este mini ejemplo muestra que las cosas van según lo planeado. Ahora ampliar:

> M <- matrix(rnorm(1e7), nrow=20) 
> system.time(invisible(apply(M, 2, t.test))) 
    user system elapsed 
101.346 0.626 101.859 
> system.time(invisible(mclapply(1:dim(M)[2], function(i) t.test(M[,i])))) 
    user system elapsed 
55.049 2.527 43.668 

Esto está utilizando 8 núcleos virtuales. Su experiencia puede ser diferente. No es una gran ganancia, pero viene de muy poco esfuerzo.

EDITAR

Si sólo se preocupan por la propia estadística t, extrayendo el campo correspondiente ($statistic) hace las cosas un poco más rápido, en particular en el caso de varios núcleos:

> system.time(invisible(apply(M, 2, function(c) t.test(c)$statistic))) 
    user system elapsed 
80.920 0.437 82.109 
> system.time(invisible(mclapply(1:dim(M)[2], function(i) t.test(M[,i])$statistic))) 
    user system elapsed 
21.246 1.367 24.107 

O aún más rápido, calcular el valor de t directamente

my.t.test <- function(c){ 
    n <- sqrt(length(c)) 
    mean(c)*n/sd(c) 
} 

Entonces

> system.time(invisible(apply(M, 2, function(c) my.t.test(c)))) 
    user system elapsed 
21.371 0.247 21.532 
> system.time(invisible(mclapply(1:dim(M)[2], function(i) my.t.test(M[,i])))) 
    user system elapsed 
144.161 8.658 6.313 

Answer 2

Usted puede hacer mejor que esto con la función de la colttestsgenefilter paquete (en Bioconductor).

> library(genefilter) 
> M <- matrix(rnorm(40),nrow=20) 
> my.t.test <- function(c){ 
+ n <- sqrt(length(c)) 
+ mean(c)*n/sd(c) 
+ } 
> x1 <- apply(M, 2, function(c) my.t.test(c)) 
> x2 <- colttests(M, gl(1, nrow(M)))[,"statistic"] 
> all.equal(x1, x2) 
[1] TRUE 
> M <- matrix(rnorm(1e7), nrow=20) 
> system.time(invisible(apply(M, 2, function(c) my.t.test(c)))) 
    user system elapsed 
27.386 0.004 27.445 
> system.time(invisible(colttests(M, gl(1, nrow(M)))[,"statistic"])) 
    user system elapsed 
    0.412 0.000 0.414 

Ref: "computación miles de estadísticas de prueba simultáneamente en R", SCGN, Vol 18 (1), 2007, http://stat-computing.org/newsletter/issues/scgn-18-1.pdf.