假设我有一个大矩阵。
M <- matrix(rnorm(1e7),nrow=20)
再假设每一列代表一个样本 假设我想应用 t.test()
到每一列,是否有一种方法比使用 apply()
?
apply(M, 2, t.test)
在我的电脑上运行分析的时间略低于2分钟。
> system.time(invisible( apply(M, 2, t.test)))
user system elapsed
113.513 0.663 113.519
如果你有一台多核机器,使用所有的核心会有一些收益,例如使用 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
这个小例子表明,事情按我们的计划进行。现在扩大规模。
> 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
这是用8个虚拟核心。你的里程数可能会有所不同。不是一个巨大的收益,但它来自于很少的努力。
编辑
如果你只关心t统计量本身,提取相应的字段($statistic
)让事情变得更快一些,尤其是在多核情况下。
> 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
或者更快,直接计算t值。
my.t.test <- function(c){
n <- sqrt(length(c))
mean(c)*n/sd(c)
}
然后
> 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
你可以做得比这更好的 colttests
函数从 genefilter
包(在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
参考:"在R中同时计算数千个测试统计数据",SCGN,第18(1)卷,2007年。http:/stat-computing.orgnewsletterissuesscgn-18-1.pdf。.