scipy.stats.

permutation_test#

scipy.stats.permutation_test(data, statistic, *, permutation_type='independent', vectorized=None, n_resamples=9999, batch=None, alternative='two-sided', axis=0, rng=None)[源代码]#

对所提供数据的给定统计量执行置换检验。

对于独立样本统计量,零假设是数据是从同一分布中随机抽样的。对于配对样本统计量,可以检验两个零假设:数据是随机配对的,或者数据是随机分配到样本中的。

参数:
dataarray-like 的可迭代对象

包含样本,每个样本都是一个观测值数组。样本数组的维度必须兼容广播,除非沿着 axis

statistic可调用对象

要计算假设检验的 p 值的统计量。statistic 必须是一个可调用对象,它接受样本作为单独的参数 (例如 statistic(*data)) 并返回结果统计量。如果 vectorized 设置为 True,则 statistic 还必须接受一个关键字参数 axis,并且被向量化以计算沿样本数组提供的 axis 的统计量。

permutation_type{‘independent’,‘samples’,‘pairings’},可选

要执行的置换类型,与零假设一致。前两种置换类型用于配对样本统计量,其中所有样本都包含相同数量的观测值,并且沿着 axis 的对应索引的观测值被认为是配对的;第三种用于独立样本统计量。

  • 'samples':观测值被分配到不同的样本,但仍与来自其他样本的相同观测值配对。此置换类型适用于配对样本假设检验,例如 Wilcoxon 符号秩检验和配对 t 检验。

  • 'pairings':观测值与不同的观测值配对,但它们仍保留在同一样本内。此置换类型适用于具有诸如 Spearman 的 \(\rho\),Kendall 的 \(\tau\) 和 Pearson 的 \(r\) 等统计量的关联/相关性测试。

  • 'independent'(默认):观测值被分配到不同的样本。样本可能包含不同数量的观测值。此置换类型适用于独立的样本假设检验,例如 Mann-Whitney \(U\) 检验和独立样本 t 检验。

    请参阅下面的“注释”部分,以获取有关置换类型的更详细说明。

vectorizedbool,可选

如果 vectorized 设置为 False,则不会将关键字参数 axis 传递给 statistic,并且预期仅计算 1D 样本的统计量。如果 True,则将关键字参数 axis 传递给 statistic,并且预期在传递 ND 样本数组时计算沿 axis 的统计量。如果 None(默认),如果 axisstatistic 的参数,则 vectorized 将设置为 True。使用向量化的统计量通常会减少计算时间。

n_resamplesint 或 np.inf,默认值:9999

用于近似零分布的随机置换(重采样)的数量。如果大于或等于不同的置换数,则将计算确切的零分布。请注意,不同的置换数随样本大小的增加而迅速增长,因此精确检验仅适用于非常小的数据集。

batchint,可选

在每次调用 statistic 中处理的置换数量。内存使用量为 O( batch * n ),其中 n 是所有样本的总大小,而与 vectorized 的值无关。默认值为 None,在这种情况下,batch 是置换的数量。

alternative{‘two-sided’,‘less’,‘greater’},可选

计算 p 值的备择假设。对于每个备择假设,p 值对于精确检验定义如下。

  • 'greater':大于或等于检验统计量的观测值的零分布的百分比。

  • 'less':小于或等于检验统计量的观测值的零分布的百分比。

  • 'two-sided'(默认):上述 p 值中较小者的两倍。

请注意,随机检验的 p 值是根据 [2][3] 中建议的保守(高估)近似值计算的,而不是 [4] 中建议的无偏估计量。也就是说,在计算与检验统计量的观测值一样极端的随机零分布的比例时,分子和分母中的值都增加 1。对此调整的一种解释是,检验统计量的观测值始终包含作为随机零分布的元素。双侧 p 值使用的约定并非普遍;返回观测的检验统计量和零分布,以防首选不同的定义。

axisint,默认值:0

(广播的)样本上计算统计量的轴。如果样本的维度数量不同,则在考虑 axis 之前,会将单例维度预加到维度较少的样本中。

rng{None, int, numpy.random.Generator},可选

如果通过关键字传递 rng,则将 numpy.random.Generator 以外的类型传递给 numpy.random.default_rng 以实例化一个 Generator。如果 rng 已经是 Generator 实例,则使用提供的实例。指定 rng 以实现可重复的函数行为。

如果此参数是通过位置传递的,或者 random_state 是通过关键字传递的,则适用于参数 random_state 的旧行为

  • 如果 random_state 为 None (或 numpy.random),则使用 numpy.random.RandomState 单例。

  • 如果 random_state 是一个整数,则会使用一个新的 RandomState 实例,并使用 random_state 作为种子。

  • 如果 random_state 已经是 GeneratorRandomState 实例,则直接使用该实例。

在 1.15.0 版本中更改: 作为从使用 numpy.random.RandomState 过渡到使用 numpy.random.GeneratorSPEC-007 过渡的一部分,此关键字已从 random_state 更改为 rng。在过渡期内,这两个关键字都将继续工作,但一次只能指定一个。过渡期结束后,使用 random_state 关键字的函数调用将发出警告。random_staterng 的行为如上所述,但在新代码中应仅使用 rng 关键字。

返回:
resPermutationTestResult

一个具有以下属性的对象

statisticfloat 或 ndarray

数据的观察检验统计量。

pvaluefloat 或 ndarray

给定备择假设的 p 值。

null_distributionndarray

在零假设下生成的检验统计量的值。

注释

此函数支持的三种排列检验类型如下所述。

非配对统计量 (permutation_type='independent')

与此排列类型相关的零假设是,所有观测值都是从相同的底层分布中采样的,并且它们已被随机分配到其中一个样本。

假设 data 包含两个样本;例如 a, b = data。当 1 < n_resamples < binom(n, k) 时,其中

  • ka 中的观测值数量,

  • nab 中的观测值总数,并且

  • binom(n, k) 是二项式系数(n 选择 k),

将数据池化(连接),随机分配给第一个或第二个样本,并计算统计量。此过程重复执行 permutation 次,生成零假设下统计量的分布。将原始数据的统计量与此分布进行比较以确定 p 值。

n_resamples >= binom(n, k) 时,执行精确检验:数据在每个不同的方式中被分割到样本之间,并且形成精确的零分布。请注意,对于给定样本之间的数据分割,仅考虑数据每个样本内的顺序/排列。对于不依赖于样本内数据顺序的统计量,这会大大降低计算成本,而不会影响零分布的形状(因为每个值的频率/计数受到相同因素的影响)。

对于 a = [a1, a2, a3, a4]b = [b1, b2, b3],此排列类型的示例是 x = [b3, a1, a2, b2]y = [a4, b1, a3]。由于在精确检验中只考虑数据每个样本内的顺序/排列,因此像 x = [b3, a1, b2, a2]y = [a4, a3, b1] 这样的重采样将被视为与上述示例不同。

permutation_type='independent' 不支持单样本统计量,但它可以应用于具有两个以上样本的统计量。在这种情况下,如果 n 是每个样本中观测值数量的数组,则不同分区的数量为

np.prod([binom(sum(n[i:]), sum(n[i+1:])) for i in range(len(n)-1)])

配对统计量,排列配对 (permutation_type='pairings')

与此排列类型相关的零假设是,每个样本中的观测值都是从相同的底层分布中抽取的,并且与其他样本元素的配对是随机分配的。

假设 data 仅包含一个样本;例如 a, = data,我们希望考虑 a 的元素与第二个样本 b 的元素的所有可能配对。设 na 中的观测值数量,它也必须等于 b 中的观测值数量。

1 < n_resamples < factorial(n) 时,a 的元素被随机排列。用户提供的统计量接受一个数据参数,比如 a_perm,并计算考虑 a_permb 的统计量。此过程重复执行 permutation 次,生成零假设下统计量的分布。将原始数据的统计量与此分布进行比较以确定 p 值。

n_resamples >= factorial(n) 时,执行精确检验:a 以每种不同的方式精确地排列一次。因此,对于 ab 之间样本的每个唯一配对,statistic 被精确地计算一次。

对于 a = [a1, a2, a3]b = [b1, b2, b3],此排列类型的示例是 a_perm = [a3, a1, a2],而 b 保持其原始顺序。

permutation_type='pairings' 支持 data 包含任意数量的样本,每个样本必须包含相同数量的观测值。所有在 data 中提供的样本都独立排列。因此,如果 m 是样本数量,n 是每个样本中观测值的数量,则精确检验中的排列次数为

factorial(n)**m

请注意,如果一个双样本统计量,例如,并非固有地依赖于观测值的提供顺序,而仅依赖于观测值的配对,则只应在 data 中提供两个样本中的一个。这会大大降低计算成本,而不会影响零分布的形状(因为每个值的频率/计数受到相同因素的影响)。

配对统计量,排列样本 (permutation_type='samples')

与这种置换类型相关的零假设是,每对中的观测值来自相同的底层分布,并且它们被分配到的样本是随机的。

假设 data 包含两个样本;例如 a, b = data。设 na 中的观测值数量,它也必须等于 b 中的观测值数量。

1 < n_resamples < 2**n 时,a 的元素和 b 的元素在样本之间随机交换(保持它们的配对),并计算统计量。这个过程重复执行,执行 permutation 次,生成零假设下统计量的分布。原始数据的统计量与此分布进行比较,以确定 p 值。

n_resamples >= 2**n 时,执行精确检验:观测值以每种不同的方式(同时保持配对)恰好分配到两个样本一次。

对于 a = [a1, a2, a3]b = [b1, b2, b3],这种置换类型的一个例子是 x = [b1, a2, b3]y = [a1, b2, a3]

permutation_type='samples' 支持 data 包含任意数量的样本,每个样本必须包含相同数量的观测值。如果 data 包含多个样本,则 data 内的配对观测值在样本之间独立交换。因此,如果 m 是样本数量,n 是每个样本内的观测值数量,那么精确检验中的置换次数是

factorial(m)**n

一些配对样本统计检验,例如 Wilcoxon 符号秩检验和配对样本 t 检验,可以仅考虑两个配对元素之间的差异来执行。因此,如果 data 仅包含一个样本,则零分布是通过独立更改每个观测值的符号来形成的。

警告

p 值是通过计算零分布中与统计量的观测值一样极端或更极端的元素来计算的。由于使用了有限精度的算术,一些统计函数在理论值完全相等时会返回数值上不同的值。在某些情况下,这可能会导致计算出的 p 值出现较大误差。permutation_test 通过考虑零分布中与检验统计量的观测值“接近”(在非精确数据类型的浮点数 epsilon 的 100 倍相对容差内)的元素,将其视为等于检验统计量的观测值来防止这种情况。但是,建议用户检查零分布以评估这种比较方法是否合适,如果不合适,则手动计算 p 值。请参见下面的示例。

参考文献

[1]
    1. Fisher. 实验设计,第 6 版 (1951)。

[2]

B. Phipson 和 G. K. Smyth。“置换 p 值永远不应为零:当随机抽取置换时计算精确 p 值。”《遗传学和分子生物学统计应用》9.1 (2010)。

[3]

M. D. Ernst。“置换方法:精确推断的基础”。统计科学 (2004)。

[4]

B. Efron 和 R. J. Tibshirani. 《自举导论》 (1993)。

示例

假设我们希望检验两个样本是否来自相同的分布。假设我们不知道底层分布,并且在观察数据之前,我们假设第一个样本的均值将小于第二个样本的均值。我们决定使用样本均值之间的差异作为检验统计量,并且我们将认为 0.05 的 p 值具有统计学意义。

为了提高效率,我们以矢量化的方式编写定义检验统计量的函数:样本 xy 可以是 ND 数组,并且将为沿 axis 的每个轴切片计算统计量。

>>> import numpy as np
>>> def statistic(x, y, axis):
...     return np.mean(x, axis=axis) - np.mean(y, axis=axis)

收集数据后,我们计算检验统计量的观测值。

>>> from scipy.stats import norm
>>> rng = np.random.default_rng()
>>> x = norm.rvs(size=5, random_state=rng)
>>> y = norm.rvs(size=6, loc = 3, random_state=rng)
>>> statistic(x, y, 0)
-3.5411688580987266

实际上,检验统计量为负,这表明 x 底层分布的真实均值小于 y 底层分布的真实均值。为了确定如果两个样本来自同一分布,这种情况偶然发生的概率,我们执行置换检验。

>>> from scipy.stats import permutation_test
>>> # because our statistic is vectorized, we pass `vectorized=True`
>>> # `n_resamples=np.inf` indicates that an exact test is to be performed
>>> res = permutation_test((x, y), statistic, vectorized=True,
...                        n_resamples=np.inf, alternative='less')
>>> print(res.statistic)
-3.5411688580987266
>>> print(res.pvalue)
0.004329004329004329

在零假设下,获得小于或等于观测值的检验统计量的概率为 0.4329%。这小于我们选择的 5% 阈值,因此我们认为这是反对零假设而支持备择假设的显著证据。

由于上述样本的大小较小,permutation_test 可以执行精确检验。对于较大的样本,我们采用随机置换检验。

>>> x = norm.rvs(size=100, random_state=rng)
>>> y = norm.rvs(size=120, loc=0.2, random_state=rng)
>>> res = permutation_test((x, y), statistic, n_resamples=9999,
...                        vectorized=True, alternative='less',
...                        rng=rng)
>>> print(res.statistic)
-0.4230459671240913
>>> print(res.pvalue)
0.0015

在零假设下,获得小于或等于观测值的检验统计量的近似概率为 0.0225%。这再次小于我们选择的 5% 阈值,因此我们再次有显著证据拒绝零假设而支持备择假设。

对于较大的样本和置换次数,结果与相应的渐近检验(独立样本 t 检验)的结果相当。

>>> from scipy.stats import ttest_ind
>>> res_asymptotic = ttest_ind(x, y, alternative='less')
>>> print(res_asymptotic.pvalue)
0.0014669545224902675

提供了检验统计量的置换分布,以供进一步研究。

>>> import matplotlib.pyplot as plt
>>> plt.hist(res.null_distribution, bins=50)
>>> plt.title("Permutation distribution of test statistic")
>>> plt.xlabel("Value of Statistic")
>>> plt.ylabel("Frequency")
>>> plt.show()
../../_images/scipy-stats-permutation_test-1_00_00.png

如果统计量由于机器精度有限而存在不准确性,则检查零分布至关重要。考虑以下情况

>>> from scipy.stats import pearsonr
>>> x = [1, 2, 4, 3]
>>> y = [2, 4, 6, 8]
>>> def statistic(x, y, axis=-1):
...     return pearsonr(x, y, axis=axis).statistic
>>> res = permutation_test((x, y), statistic, vectorized=True,
...                        permutation_type='pairings',
...                        alternative='greater')
>>> r, pvalue, null = res.statistic, res.pvalue, res.null_distribution

在这种情况下,由于数值噪声,零分布的某些元素与相关系数 r 的观测值不同。我们手动检查零分布中与检验统计量的观测值几乎相同的元素。

>>> r
0.7999999999999999
>>> unique = np.unique(null)
>>> unique
array([-1. , -1. , -0.8, -0.8, -0.8, -0.6, -0.4, -0.4, -0.2, -0.2, -0.2,
    0. ,  0.2,  0.2,  0.2,  0.4,  0.4,  0.6,  0.8,  0.8,  0.8,  1. ,
    1. ])  # may vary
>>> unique[np.isclose(r, unique)].tolist()
[0.7999999999999998, 0.7999999999999999, 0.8]  # may vary

如果 permutation_test 要天真地执行比较,则零分布中值为 0.7999999999999998 的元素不会被认为与统计量的观测值一样极端或更极端,因此计算出的 p 值会太小。

>>> incorrect_pvalue = np.count_nonzero(null >= r) / len(null)
>>> incorrect_pvalue
0.14583333333333334  # may vary

相反,permutation_test 将零分布中与统计量观测值 rmax(1e-14, abs(r)*1e-14) 范围内的元素视为等于 r

>>> correct_pvalue = np.count_nonzero(null >= r - 1e-14) / len(null)
>>> correct_pvalue
0.16666666666666666
>>> res.pvalue == correct_pvalue
True

这种比较方法有望在大多数实际情况下都是准确的,但建议用户通过检查零分布中接近统计量观测值的元素来评估这一点。另外,请考虑使用可以使用精确算术计算的统计量(例如,整数统计量)。