尝试对每一行执行t.test，并计算p值小于0.0 5的所有行。

Question

在过去的四个小时里，我一直在绞尽脑汁，试图找到一个让我发疯的R问题的解决方案。我到处寻找一个像样的答案，但到目前为止，我一直在一堵又一堵的墙上。我现在恳请您对这一优秀社区的善意帮助。

考虑以下数据集：

set.seed(2112)
DataSample <- matrix(rnorm(24000),nrow=1000)
colnames(DataSample) <- c(paste("Trial",1:12,sep=""),paste("Control",13:24,sep=""))

我需要对DataSample中的每一行执行t检验，以确定组、试用和控制是否不同(适用相同的方差)。

然后，我需要计算p值等于或低于的行数。

下面是我尝试过的代码，我知道这是错误的：

set.seed(2112)
DataSample <- matrix(rnorm(24000),nrow=1000)
colnames(DataSample) <- c(paste("Trial",1:12,sep=""),paste("Control",13:24,sep=""))

pValResults <- apply(
  DataSample[,1:12],1,function(x) t.test(x,DataSample[,13:24], var.equal=T)$p.value
  )

sum(pValResults < 0.05) # Returns the wrong answer (so I was told)

我确实试着看了很多类似的关于堆栈溢出的问题，但是我经常会出现语法错误或维度不匹配。上面的代码是我所能得到的最好的，而不返回一个R错误--但是由于代码返回了错误的答案，我没有什么值得骄傲的。

任何建议都将不胜感激！提前谢谢你的时间。

Answer 1

一种选择是循环遍历计算每一行t测试的数据集，但它不那么优雅。

set.seed(2112)
DataSample <- matrix(rnorm(24000),nrow=1000)
colnames(DataSample) <- c(paste("Trial",1:12,sep=""),paste("Control",13:24,sep=""))

# initialize vector of stored p-values
pvalue <- rep(0,nrow(DataSample))

for (i in 1:nrow(DataSample)){
   pvalue[i] <- t.test(DataSample[i,1:12],DataSample[i,13:24])$p.value
}
# finding number that are significant
sum(pvalue < 0.05)

Answer 2

我转换成了data.table，我得到的答案是45：

DataSample.dt <- as.data.table(DataSample)
sum(sapply(seq_len(nrow(DataSample.dt)), function(x)
    t.test(DataSample.dt[x, paste0('Trial', 1:12), with=F],
           DataSample.dt[x, paste0('Control', 13:24), with=F],
           var.equal=T)$p.value) < 0.05)

Answer 3

要执行配对T检验，需要提供paired = TRUE参数。t.test函数不是矢量化的，但是一次测试整个矩阵非常简单。以下是三种方法(包括使用apply)：

library("genefilter")
library("matrixStats")
library("microbenchmark")
dd <- DataSample[, 1:12] - DataSample[, 13:24]
microbenchmark::microbenchmark(
  manual = {ps1 <- 2 * pt(-abs(rowMeans(dd) / sqrt(rowVars(dd) / ncol(dd))), ncol(dd) - 1)},
  apply = {ps2 <- apply(DataSample, 1, function(x) t.test(x[1:12], x[13:24], paired=TRUE)$p.value)},
  rowttests = {ps3 <- rowttests(dd)[, "p.value"]})
#Unit: milliseconds
#      expr        min         lq       mean     median         uq        max
#    manual   1.611808   1.641783   1.677010   1.663122   1.709401   1.852347
#     apply 390.869635 398.720930 404.391487 401.508382 405.715668 634.932675
# rowttests   2.368823   2.417837   2.639671   2.574320   2.757870   7.207135
# neval
#   100
#   100
#   100

您可以看到，手动方法比应用程序快200倍以上。

如果您实际上指的是未配对测试，下面是等效的比较：

microbenchmark::microbenchmark(
  manual = {x <- DataSample[, 1:12]; y <- DataSample[, 13:24]; ps1 <- 2 * pt(-abs((rowMeans(x) - rowMeans(y)) / sqrt((rowVars(x) + rowVars(y)) / ncol(x))), ncol(DataSample) - 2)},
  apply = { ps2 <- apply(DataSample, 1, function(x) t.test(x[1:12], x[13:24], var.equal = TRUE)$p.value)},
  rowttests = {ps3 <- rowttests(DataSample, factor(rep(1:2, each = 12)))[, "p.value"]})

注意，手动方法假定这两个组的大小相同。