配对 t 检验崩溃应用循环（已编辑）

Posted 2023-03-12

【中文标题】配对 t 检验崩溃应用循环（已编辑）

【英文标题】：paired t-test crashes apply-loop (edited)

【发布时间】：2012-03-17 20:39:28

【问题描述】：

为了回应有用的 cmets，我编辑了原始问题（我曾假设 for-loop 和 apply-loop 给出不同的结果）。

我正在使用 R 运行大量 2 组 t 检验，使用来自分隔表的输入。根据这里和其他地方的建议，我尝试了“for-loops”和“apply”来实现这一点。对于“正常”的 t.test，两者都可以很好地工作并给出相同的结果。但是，对于配对 t 检验，for-look 似乎有效，而 apply-loop 则无效。后来，我发现两个循环都遇到了同样的问题（见下文），但 for 循环更优雅地处理了这种情况（循环只有一个循环返回无效结果），而 apply 循环完全失败。

我的输入文件如下所示：（第一行是标题行，数据行有一个名称，第 1 组有 4 个数据点，第 2 组有 4 个数据点）：

header g1.1 g1.2 g1.3 g1.4 g2.1 g2.2 g2.3 g2.4
name1  0    0.5  -0.2 -0.2 -0.1 0.4 -0.3 -0.3
name2  23.2 24.4 24.5 27.2 15.5 16.5 17.7 20.0
name3  .....

等等（总共约 50000 行）。第一个数据行（以name19开头的竟然是罪魁祸首。

这是效果更好的 for 循环版本（在有问题的行上失败，但可以正确处理所有其他行）：

table <- read.table('ttest_in.txt',head=1,sep='\t')
for(i in 1:nrow(table)) 
   g1<-as.numeric((table)[i,2:5])
   g2<-as.numeric((table)[i,6:9])
   pv <- t.test(g1,g2,paired=TRUE)$p.value

这是导致问题的“应用”版本

table <- read.table('ttest_in.txt',head=1,sep='\t')
pv.list <- apply(table[,2:9],1,function(x)t.test(x[1:4],x[5:8],paired=TRUE)$p.value)

大约 50000 条数据行中的一条是有问题的，因为所有成对比较的差异都是相同的，这在成对 t 检验中会导致未定义的 p 值（基本上为零）。应用循环因错误“数据基本不变”而崩溃。对我（作为一个 R 新手）来说，仅仅因为 t.test 不喜欢一个数据而使整个脚本崩溃似乎不是一个好主意。在 for 循环中，此数据行也会导致错误消息，但循环继续，所有其他 t 检验都给出正确的结果。

我做错了什么吗？这种行为本质上禁止使用应用循环进行这种批处理分析。或者有没有一种标准的方法来规避这个问题。为什么 t 检验不只是返回对该特定 p 值无效的东西而不是退出？

【问题讨论】：

您能否提供一个重现错误的最小示例？我们可能需要比您上面提供的更多的帮助。但是，不要在apply() 版本中使用as.matrix() 行。如果header 确实在文件中，那么您可能会遇到麻烦。 apply() 在数据框上工作得很好，所以你可以使用pv.list <- apply(table[, 2:9], 1 ,function(x) t.test(x[1:4], x[5:8], paired=TRUE)$p.value) 而不是在矩阵mat 上工作。这能解决问题吗？

添加到 Gavin 的优秀问题：尝试在打开 debug 的情况下运行您的 apply 代码。可能会出现一些半明显的索引错误。

感谢您的 cmets。 Gavin 建议跳过 as.matrix 步骤简化了脚本，但没有改变任何内容。我不熟悉“调试”，以后可能会尝试。接下来，我追查了罪魁祸首，现在我（部分）不明白这个问题。我已经更改了上面的示例文件，它是导致问题的第一行。我可以看到为什么配对 t 检验不喜欢这些数据，但是这个原因足以使整个“应用”循环崩溃吗？对此我能做些什么吗？

我应该补充一点，在理解了问题之后，我看到 for 循环也拒绝了这个特定的数据行。我没有注意到这一点，因为循环的其余部分工作正常，输出文件看起来不显眼（只缺少一个输出行）

后审工作做得很好。您可以编辑您的问题或将其删除以反映您更新的问题吗？

【参考方案1】：

更新既然你说for循环也给出错误并且你希望apply版本更健壮，为什么不简单地添加一个tryCatch？

pv.list <- apply(table[,2:9],1, function(x) tryCatch( 
  t.test(x[1:4],x[5:8],paired=TRUE)$p.value, error=function(x) NA ))

如果 p.value 无法计算，这应该返回 NA。您可以通过编辑错误处理函数更改为另一个值（例如NULL、0 或Inf）。

旧帖

我注意到t.test（有点）在某些值为Inf（这似乎是一个错误）时给出了您发现的错误：

> t.test(1:10, c(rep(1,9), Inf), paired=TRUE)
Error in if (stderr < 10 * .Machine$double.eps * abs(mx)) stop("data are essentially constant") : 
missing value where TRUE/FALSE needed

所以你真的得到这个还是它真的说：

Error in t.test: data are essentially constant

仍然不太清楚为什么 for 循环有效。但请注意，在您的 for 循环中您执行 as.numeric，而在 apply 案例中您不会执行此操作...

【讨论】：

我同意这种行为也是不可取的（但在数据中包含 Inf 比导致我的问题的看起来无害的数据线更容易发现）。 t.test 似乎不适合批处理（或者必须采取更多预防措施）。顺便说一句，for-loop 对于有问题的行也失败了，但至少会继续并处理所有其他数据行

@user900889 - 啊哈！请参阅我的更新答案以获取可能的解决方案。

为什么我没有简单地添加 TryCatch？因为我不知道它存在:-)。谢谢！

【参考方案2】：

在这种情况下，我会捕获所有警告和错误并在之后进行调查，如下所示：How do I save warnings and errors as output from a function?

您还可以在这里找到一些好主意：How to tell lapply to ignore an error and process the next thing in the list?

【参考方案3】：

t.test(rep(1,4),rep(0,4))会报同样的错误。而https://stat.ethz.ch/pipermail/r-help/2008-February/154167.html是处理1）零方差2）观察数不≥1的答案

【讨论】：

对不起，我的意思是可能是导致问题发生的零方差，我只是举一个简单的例子来简单通知。我会完整回答

【参考方案4】：

我正好遇到了这个问题，使用 tapply 运行大量配对 t 检验。我使用样本量作为指导来删除我可能不应该在任何情况下运行 t-tests 的数据。