Qual è il modo più veloce per applicare t.test a ciascuna colonna di una matrice di grandi dimensioni?

Question

Supponiamo di avere una grande matrice:

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

Inoltre supponiamo che ciascuna colonna rappresenta un campione. Dire che vorrei applicare t.test() a ciascuna colonna, c'è un modo per farlo che è molto più veloce rispetto all'utilizzo di apply()?

apply(M, 2, t.test) 

Ci sono voluti un po 'meno di 2 minuti per eseguire l'analisi sul mio computer:

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

Answer 1

Se si dispone di una macchina multicore ci sono alcuni guadagni da utilizzare tutti i core, ad esempio utilizzando 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 

Questo mini esempio mostra che le cose vanno come pianificato. Ora scalare:

> 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 

Questo sta utilizzando 8 core virtuali. Il tuo chilometraggio può variare. Non è un guadagno enorme, ma viene da un piccolo sforzo.

EDIT

Se vi interessa soltanto la statistica t per sé, l'estrazione del campo corrispondente ($statistic) rende le cose un po 'più veloce, in particolare nel caso multicore:

> 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 

Or ancora più veloce, calcolare il valore t direttamente

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

Poi

> 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

È possibile fare meglio di questo con la funzione colttests dal pacchetto genefilter (su 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 

Rif: "Informatica migliaia di statistiche test contemporaneamente in R", SCGN, Vol 18 (1), 2007, http://stat-computing.org/newsletter/issues/scgn-18-1.pdf.