尝试对每一行执行 t.test 并计算 p 值小于 0.05 的所有行答案

【问题标题】：trying to perform a t.test for each row and count all rows where p-value is less than 0.05尝试对每一行执行 t.test 并计算 p 值小于 0.05 的所有行
【发布时间】：2015-10-14 07:23:24
【问题描述】：

在过去的四个小时里，我一直在努力寻找 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 检验，以确定 TRIAL 和 CONTROL 组是否不同（适用等方差）。

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

所以这是我试过的代码，我知道这是错误的：

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)

我确实尝试查看有关 stackoverflow 的许多类似问题，但我经常会遇到语法错误或尺寸不匹配。上面的代码是在不返回 R 错误的情况下我能得到的最好的代码——但由于代码返回了错误的答案，我没有什么好骄傲的。

任何建议将不胜感激！提前感谢您的宝贵时间。

【问题讨论】：

你只需要修正你的申请语句。 pvalresults
哈！您的代码返回 56（这是使用时的正确答案，请参阅 879）。在研究答案时，我偶然发现了这段代码，但不幸的是它返回了 55。我想知道这两者有什么区别？我确实对在 apply 中使用 FUNCTION 感到困惑（并且 apply 对我来说也是新闻。这是我被告知要使用的代码：pValResults <- apply(ME,1, function(x) t.test(x[1:12],x[13:24])$p.value) #this is what I used
唯一的区别是 t.test 默认使用 var.equal 为 False。我的代码（以及上面的示例代码）使用 var.equal = T 意味着 t.test 以不同的方式计算 p 值。检查 t.test 文档，它解释了它。
您是说配对吗？到目前为止，你得到的答案都没有给你。还值得注意的是，var.equal 参数对配对 t 检验没有任何作用。
这能回答你的问题吗？ doing t.test for columns for each row in data set

标签： r statistics hypothesis-test

【解决方案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)

【讨论】：

有趣。您的代码返回的值与我在 IRC 聊天中建议的另一个代码相同。考虑使用 set.seed(879) 而不是 2112（我必须使用的真正的），您的代码（以及我在 IRC 上获得的代码）返回 55。但我被告知正确答案是 56。
好吧，我肯定出了问题。我关闭了 R 然后重新打开它，现在你的代码突然返回了正确的值。您的代码中缺少的一件事是 var.equal=T，加上总和意味着 pvalue

【解决方案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)

【讨论】：

是的，45 是使用种子 2112 时的值。使用种子 879 时，我本来应该得到 56，而我得到了 55，除了上面的 ARobertson 代码实际上返回 56。现在不管原因我后来重新启动了 R 并重试了所有代码，所有代码都给出了 56。但我可以发誓我没有对数据集进行任何更改以导致这种差异发生......所以我想我会衰老。跨度>

【解决方案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"]})

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

【讨论】：

天哪！我真的很抱歉我在标题上犯了一个严重的错误。你是对的，但这并不意味着是配对测试。我的意思是说 var.equal=T 和我的电线交叉了。

【解决方案4】：

使用外部库添加替代方案。

执行测试：

library(matrixTests)
res <- row_t_equalvar(DataSample[,1:12], DataSample[,13:24])

结果格式：

res

   obs.x obs.y obs.tot      mean.x       mean.y    mean.diff     var.x     var.y var.pooled    stderr df    statistic     pvalue   conf.low   conf.high alternative mean.null conf.level
1     12    12      24  0.30569721  0.160622830  0.145074376 0.5034806 1.0769678  0.7902242 0.3629105 22  0.399752487 0.69319351 -0.6075559  0.89770469   two.sided         0       0.95
2     12    12      24 -0.27463354 -0.206396781 -0.068236762 0.8133311 0.2807800  0.5470556 0.3019535 22 -0.225984324 0.82329990 -0.6944500  0.55797651   two.sided         0       0.95
3     12    12      24 -0.19805092 -0.023207888 -0.174843032 0.4278359 0.5604078  0.4941219 0.2869733 22 -0.609265949 0.54858909 -0.7699891  0.42030307   two.sided         0       0.95

p <= 0.05 的行数：

> sum(res$pvalue <= 0.05)
[1] 4

【讨论】：