在排列测试中使用秩和测试

Question

我对 scipy 提供的 permutation_test 函数的定义有点困惑。 以下是我写的用于计算 p 和 null 分布的内容

x = disease['Theta']
y = nodisease['Theta']

def ranksum_statistic(x, y, axis):
    stat, _ = ranksums(x, y, axis = axis, alternative='greater')
    return stat

permutation_test(x, y, ranksum_statistic, permutation_type='independent', 
vectorized=None, n_resamples=9999, batch=None, 
alternative='greater', axis=0, random_state=None)

这将返回一个与直接使用ranksum测试完全不同的p值

(-log10p = 3)

（另外，如果我多次运行ranksum测试，p值也会改变很多）如果我想比较来自相同分布（值不同）的两个不同大小的数组（非参数），我很困惑哪个更合适，因为无疾病应该已经是零分布。

Answer 1

在下文中，请注意

(-log10p = 150) 执行与

ranksums

基本相同的测试（Mann Whitney U 测试）。不同之处在于

mannwhitneyu

报告统计量的 z 分数，而不是积分

ranksums

统计量，并且它无法像

mannwhitneyu

那样执行连续性校正并计算精确的 p 值，具体取决于您的参数。所以我总是建议使用

mannwhitneyu

，但根据要求，我会用

mannwhitneyu

来回答，仅使用

ranksums

作为参考。 您的代码基本上是正确的。为了进行检查，我们将

mannwhitneyu

给出的 p 值与

permutation_test

提供的 p 值与

mannwhitneyu

进行比较。method='exact'

请注意，我对您的代码所做的唯一有意义的更改是将 

import numpy as np
from scipy import stats

rng = np.random.default_rng(32590345094235)
x = rng.random(6)
y = rng.random(7)

def ranksum_statistic(x, y, axis):
    stat, _ = stats.ranksums(x, y, axis = axis, alternative='greater')
    return stat

res = stats.permutation_test((x, y), ranksum_statistic, permutation_type='independent', 
                             n_resamples=9999, alternative='greater', axis=0,
                             random_state=rng)
ref = stats.mannwhitneyu(x, y, alternative='greater', method='exact')
np.allclose(res.pvalue, ref.pvalue, rtol=1e-14)

作为元组而不是作为单独的参数传递。 （这是必需的；

(x, y)

不接受单独的样本作为单独的参数。）对于大小为 6 和 7 的小样本，我们得到精确的 p 值，因为

permutation_test

足以执行所有

n_resamples=9999

排列。因此，在这种情况下我们可以获得的最小可能 p 值是

binom(6 + 7, 6) = binom(6 + 7, 7) = 1716

。如果

1/1716

不足以执行精确检验，则会执行随机检验，在这种情况下，当统计量的观测值比所有其他排列的观测值更极端时，您可以获得的最小 p 值是

n_resamples

。如果

1 / (n_resamples + 1)

，则意味着从

n_resamples=9999

可以获得的最小 p 值为

permutation_test

。是的，p 值是随机的。如果您需要重现性，请像我一样明确传递 NumPy

-log10(p) = 3

(

Generator

)。请注意，从长远来看（即许多单独的数据集），这些随机 p 值仍然几乎精确地控制假阳性率，即使它们只是给定数据的近似值。有关更多信息，请参阅

排列 p 值永远不应该为零

。

random_state=rng

始终使用零分布的分析近似，因此它可以为您提供更小的 p 值。这“可能”更准确，但前提是您的样本量“足够大”，使得渐近近似有效。要评估您获得的微小 p 值是否合理，请注意 Mann Whitney U 检验的最大可能统计量是

ranksums

，其中 m*n 和

m

是样本量。如果这是您观察到的数据的统计量，那么确切的 p 值将是

n

。如果

m * n / binom(m + n, m)

返回的 p 值小于该值，则它太小。