【r<-方案】R检验中的“数据是恆量”问题
这是一般做基因差异表达分析在使用t检验或者其他统计检验中常出现的一个问题。之前我学习和自己分析时就遇到过,尝试使用判断的方式事先检查它是不是数据存在问题(这类数据明显不服从正态分布),可以使用正态性检验,或者直接判断是不是样本组内的数据是完全一样的,如果一样就不要这个了。
今天一个师弟问到这个问题,我才有想起来,在简书上找了下自己的笔记发现没有把之前的方案记录下来~~考虑到这个问题的普遍性,我在这篇文章里提供下方案。
所遇到的问题:
分析两个样本之间是否存在差异,每个样本三个重复。现在用的是t.test,但有些样本三个重复的值一样(比如有0,0,0或者2,2,2之类的),想问下像这种数据应该用什么检验方法呢?
举个例子:
> t.test(c(0,0,0), c(2,2,2))
Error in t.test.default(c(0, 0, 0), c(2, 2, 2)) : 数据是恆量
这就是最简单的一个重复例子了,我们需要解决的就是这个问题。
为什么出现这问题?如果解决?以下是我的回答:
数据是恒量是无法做t检验的,因为计算公式分母为0(不懂的看下统计量t的计算公式,一般标准差/标准误为分母,所以恒量是不能算的)。因为你要用t检验,我给你一个处理思路, 先不分组别,按基因名检查所有样本的基因表达值(循环)是否一样,如果一样就丢掉,如果不一样,则按组别判断样本(每组3个)基因表达是否一样,如果不一样进行t检验寻找一批差异基因,如果一样,则输出原始的结果,再筛选其中差异大的基因 。
假设有两万个基因的表达,我手头没数据,所以写个伪代码:
下面用geneExpr1与geneExpr2表示两组数据:
for循环1(geneExpr1, geneExpr2):
组合某基因表达 - c(geneExpr1, geneExpr2) -> mergedExpr
if (mergedExpr是恒量):
进行下一个循环,计算下一个基因表达差异,这个基因不算了
else:
if (geneExpr1与geneExpr2都是恒量):
输出该结果进行人为检查,可以赋给一个列表什么的。虽然两者都是恒量,但两者可能有差异,却不能用统计检验算。
else:
统计检验
在使用t
检验前尽量使用方差分析检验方差同质性。
最后提供两个参考函数:
1是判断恒量:
zero_range <- function(x, tol = .Machine$double.eps ^ 0.5) {
if (length(x) == 1) return(TRUE)
x <- range(x) / mean(x)
isTRUE(all.equal(x[1], x[2], tolerance = tol))
}
2是修正的t检验,其他检验可以参考该方式:
my.t.test.p.value <- function(...) {
obj<-try(t.test(...), silent=TRUE)
if (is(obj, "try-error")) return(NA) else return(obj$p.value)
}
这个函数可以帮助顺利的执行循环,如果出问题,返回相应的NA
,这样我们可以算完后再检查数据。
参考:
https://stackoverflow.com/questions/9508518/why-are-these-numbers-not-equal
https://stackoverflow.com/questions/23093095/t-test-failed-in-r