置换检验中的 P 值等于 0

机器算法验证 p 值 置换检验
2022-02-02 01:08:06

我有两个数据集,我想知道它们是否显着不同(这来自“两组显着不同?测试使用”)。

我决定使用置换测试,在 R 中执行以下操作:

permutation.test <- function(coding, lncrna) {
    coding <- coding[,1] # dataset1
    lncrna <- lncrna[,1] # dataset2

    ### Under null hyphotesis, both datasets would be the same. So:
    d <- c(coding, lncrna)

    # Observed difference
    diff.observed = mean(coding) - mean(lncrna)
    number_of_permutations = 5000
    diff.random = NULL

    for (i in 1:number_of_permutations) {
        # Sample from the combined dataset
        a.random = sample (d, length(coding), TRUE)
        b.random = sample (d, length(lncrna), TRUE)
        # Null (permuated) difference
        diff.random[i] = mean(b.random) - mean(a.random)
    }

    # P-value is the fraction of how many times the permuted difference is equal or more extreme than the observed difference
    pvalue = sum(abs(diff.random) >= abs(diff.observed)) / number_of_permutations
    pvalue
}

然而,根据本文,p 值不应为 0:http: //www.statsci.org/smyth/pubs/permp.pdf

你建议我做什么?这是计算p值的方法:

pvalue = sum(abs(diff.random) >= abs(diff.observed)) / number_of_permutations

一个好方法?还是做以下更好?

pvalue = sum(abs(diff.random) >= abs(diff.observed)) + 1 / number_of_permutations + 1
2个回答

讨论

置换测试生成数据集的所有相关置换,为每个这样的置换计算指定的测试统计量,并在得到的统计量置换分布的上下文中评估实际的测试统计量。评估它的一种常用方法是报告(在某种意义上)与实际统计数据“一样或更极端”的统计数据的比例。这通常被称为“p 值”。

因为实际的数据集是这些排列之一,所以它的统计数据必然在排列分布中发现的那些。因此,p 值永远不会为零。

除非数据集非常小(通常少于 20-30 个总数)或测试统计量具有特别好的数学形式,否则生成所有排列是不切实际的。(生成所有排列的示例出现在 R 中的 Permutation Test 中。)因此,排列测试的计算机实现通常从排列分布中采样。他们通过生成一些独立的随机排列来做到这一点,并希望结果是所有排列的代表性样本。

因此,从此类样本中得出的任何数字(例如“p 值”)仅是置换分布属性的估计值。估计的p 值为零是很可能的——并且经常发生在影响很大的时候。这并没有错,但它立即提出了迄今为止被忽视的问题,即估计的 p 值与正确的 p 值相差多少? 因为比例的抽样分布(例如估计的 p 值)是二项式的,所以可以使用二项式置信区间来解决这种不确定性。


建筑学

一个结构良好的实施将在所有方面密切关注讨论。 它将从计算检验统计量的例程开始,作为比较两组均值的例程:

diff.means <- function(control, treatment) mean(treatment) - mean(control)

编写另一个例程以生成数据集的随机排列并应用测试统计量。这个接口允许调用者提供测试统计数据作为参数。它将m数组的第一个元素(假定为参考组)与其余元素(“治疗”组)进行比较。

f <- function(..., sample, m, statistic) {
  s <- sample(sample)
  statistic(s[1:m], s[-(1:m)])
}

排列测试首先通过找到实际数据的统计量(这里假设存储在两个数组control和中treatment)然后找到其中许多独立随机排列的统计量来进行:

z <- stat(control, treatment) # Test statistic for the observed data
sim<- sapply(1:1e4, f, sample=c(control,treatment), m=length(control), statistic=diff.means)

现在计算 p 值的二项式估计及其置信区间。一种方法使用包中的内置binconf过程HMisc

require(Hmisc)                                    # Exports `binconf`
k <- sum(abs(sim) >= abs(z))                      # Two-tailed test
zapsmall(binconf(k, length(sim), method='exact')) # 95% CI by default

将结果与另一个测试进行比较并不是一个坏主意,即使已知这不太适用:至少您可能会对结果应该在哪里有一个数量级的感觉。在这个例子中(比较均值),Student t-test 通常给出一个好的结果:

t.test(treatment, control)

此架构在更复杂的情况下使用工作R代码进行了说明,测试变量是否遵循相同分布


例子

作为测试,我正态分布的“控制”值,从正态分布的“治疗”值100201.5

set.seed(17)
control <- rnorm(10)
treatment <- rnorm(20, 1.5)

在使用前面的代码运行排列测试后,我绘制了排列分布的样本以及一条垂直的红线来标记实际的统计数据:

h <- hist(c(z, sim), plot=FALSE)
hist(sim, breaks=h$breaks)
abline(v = stat(control, treatment), col="Red")

数字

二项式置信限计算导致

 PointEst Lower        Upper
        0     0 0.0003688199

换句话说,估计的p 值正好为零,置信区间从 (默认为 95%) 。学生 t 检验报告的 p 值为,这与此一致。这支持了我们更细致入微的理解,即在这种情况下估计的 p 值为零对应于一个非常小的 p 值,我们可以合理地认为它小于该信息虽然不确定,但通常足以对假设检验做出明确结论(因为远低于的常见阈值)。00.000373.16e-050.000370.000370.050.010.001


评论

排列分布样本中的N 个个被认为是“极端”时,都是对真实 p 值的合理估计。(其他估计也是合理的。)通常没有理由偏爱其中之一。如果它们导致不同的决定,那意味着太小了。取更大的置换分布样本,而不是捏造估计 p 值的方式。kN k/N(k+1)/(N+1)N

如果需要更高的估计精度,只需运行更长时间的排列测试。因为置信区间宽度通常与样本大小的平方根成反比,为了将置信区间提高倍,我运行了倍的排列。这次估计的 p 值为(其中五个排列结果至少与实际统计量一样远离零),置信区间为10102=1000.0000051.611.7百万分之几:比报告的学生 t 检验略小。尽管数据是使用正态随机数生成器生成的,这可以证明使用学生 t 检验是合理的,但排列检验结果与学生 t 检验结果不同,因为每组观察值内的分布并不完全正态。

由于使用估计的 p 值来决定是否拒绝原假设,因此重要的是要考虑估计量的选择如何影响错误拒绝的概率。Smyth & Phipson's 引用的论文指出,无偏估计器 ( ) 无法正确控制 I 类错误率。相比之下, ( ) 是一个有效的(但保守的)p 值估计量——它不会导致对空值的过度拒绝。BMB+1M+1

(B 是随机排列的数量,其中获得的统计量大于或等于观察到的统计量,M 是抽样的随机排列的总数)。

Smyth & Phipson 还证明了 ( ) 的无效性在多重比较设置中变得至关重要,其中导出非常小的 p 值估计值,然后通过乘以一个因子进行校正。在这些设置中,零 p 值的估计值尤其是灾难性的,因为无论应用何种校正,它都保持为零。BM