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, random_state=None)[源代码]#

对给定数据上特定统计量的置换检验。

对于独立样例统计量,原假设是数据是随机从同一分布中抽取的。对于成对样例统计量,可以检验两个原假设:一是数据是按随机配对的,二是数据是按随机分配到样例的。

参数:
data类似于数组的可迭代对象

包含样本,其中每个样本是观测值数组。样本数组的维度必须兼容遍历广播,但沿着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 检验。

    有关排列类型更详细的描述,请参阅下面的注释部分。

vectorized布尔值,可选

如果将 vectorized 设置为 False,则 statistic 将不会传递关键字参数 axis,并且仅会计算一维样本的统计值。如果将 Truestatistic 将传递关键字参数 axis,并且在传入 ND 采样数组时会按 axis 计算统计值。如果将 None (默认值),则如果 axisstatistic 的参数,则会将 vectorized 设置为 True。使用矢量化统计值通常能缩短计算时间。

n_resamplesint 或 np.inf,默认值:9999

用于近似零分布的随机置换(重新采样)的数量。如果随机置换的数量大于或等于不同的置换数量,则将计算确切零分布。请注意,不同的置换数量会随着样本数量而急剧增加,因此只有对于非常小的数据集才能执行确切检验。

batchint,可选

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

alternative{‘双边’、‘小于’、‘大于’},可选

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

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

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

  • '双边' (默认值):上面较小 p 值的两倍。

请注意:随机化测试的 p 值是基于在 [2][3] 中建议的保守(高估)近似值计算的,而不是基于在 [4] 中建议的无偏估计量,即,在计算与观察到的检验统计量值一样极端的随机化零分布部分时,分子和分母中的值都会增加 1。此调整的一项释义是,观察到的检验统计量值始终包含在随机化零分布中。双侧 p 值所用的惯例并不通用;观察到的检验统计量和零分布将根据是偏好的不同定义返回。

axisint,默认值:0

用于计算统计量的(广播)样本的轴。如果样本具有不同的维度数,那么在考虑axis之前,维度少的样本将被填充单例维度。

random_state{None、int、numpy.random.Generator,

用于生成排列的伪随机数生成器状态。

如果random_stateNone(默认值),则使用numpy.random.RandomState单例。如果random_state是 int,则使用新的RandomState实例,并用random_state进行种子填充。如果random_state已经是GeneratorRandomState实例,那么使用该实例。

返回:
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 的统计量。这一过程反复执行,排列次,根据无效假设生成统计量的分布。原始数据的统计量与该分布进行比较以确定 p 值。

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

对于 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',
...                        random_state=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 把零分布中与统计量 r 的观测值相差在 max(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

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