Python SciPy stats.permutation_test用法及代码示例

本文简要介绍 python 语言中 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： 类似数组的可迭代

包含样本，每个样本都是一组观察值。样本数组的维度必须与广播兼容，除了沿轴。

statistic： 可调用的

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

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

根据原假设，要执行的排列类型。前两种排列类型用于配对样本统计，其中所有样本包含相同数量的观察值，并且沿轴具有相应索引的观察值被认为是配对的；三是独立样本统计。

• 'samples' ：观察值分配给不同的样本，但仍与来自其他样本的相同观察值配对。此排列类型适用于配对样本假设检验，例如 Wilcoxon signed-rank 检验和配对 t-test。

• 'pairings' ：观察与不同的观察配对，但它们保留在同一个样本中。此排列类型适用于具有统计信息的关联/相关测试，例如 Spearman 的 \(\rho\) 、Kendall 的 \(\tau\) 和 Pearson 的 \(r\) 。

• 'independent'(默认)：观察被分配给不同的样本。样本可能包含不同数量的观察值。此排列类型适用于独立样本假设检验，例如 Mann-Whitney \(U\) 检验和独立样本 t-test。

  有关排列类型的更详细说明，请参阅下面的注释部分。

vectorized： 布尔型，可选

如果矢量化已设置False,统计不会传递关键字参数轴并且预计仅计算一维样本的统计量。如果True,统计将传递关键字参数轴并预计将计算统计数据轴当传递 ND 样本数组时。如果None(默认)，矢量化将被设置True如果axis是一个参数统计。使用矢量化统计通常可以减少计算时间。

n_resamples： int 或 np.inf，默认值：9999

用于近似零分布的随机排列(重采样)的数量。如果大于或等于不同排列的数量，将计算精确的零分布。请注意，不同排列的数量随着样本大小的增加而快速增长，因此精确测试仅适用于非常小的数据集。

batch： 整数，可选

每次调用中要处理的排列数统计.内存使用量为 O(批*n)， 在哪里n是所有样本的总大小，与矢量化.默认为None， 在这种情况下batch是排列的数量。

alternative： {‘双面’，‘less’, ‘greater’}，可选

计算 p 值的备择假设。对于每个替代方案，精确检验的 p 值定义如下。

• 'greater' ：大于或等于检验统计量观察值的空分布百分比。

• 'less' ：小于或等于检验统计量观察值的空分布百分比。

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

请注意，随机检验的 p 值是根据 [2] 和 [3] 中建议的保守 (over-estimated) 近似计算的，而不是根据 [4] 中建议的无偏估计量计算的。也就是说，当计算与检验统计量的观测值一样极端的随机零分布的比例时，分子和分母中的值都加一。这种调整的解释是检验统计量的观测值始终作为随机零分布的元素包含在内。用于两侧 p 值的约定并不通用；如果首选不同的定义，则返回观察到的测试统计量和零分布。

axis： 整数，默认值：0

计算统计数据的(广播)样本的轴。如果样本具有不同数量的维度，则在考虑轴之前，会将单一维度添加到具有较少维度的样本之前。

random_state： {无，整数， numpy.random.Generator ，

numpy.random.RandomState}, optional

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

如果random_state是None(默认)，numpy.random.RandomState使用单例。如果random_state是一个 int，一个新的RandomState使用实例，播种random_state.如果random_state已经是一个Generator或者RandomState实例然后使用该实例。

返回 ：：

res： PermutationTestResult

具有属性的对象：

统计 浮点数或 ndarray

观察到的数据检验统计量。

p值 浮点数或 ndarray

给定替代方案的 p 值。

null_distribution ndarray

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

注意：

该函数支持的三种置换测试如下所述。

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

与这种排列类型相关的零假设是所有观测值都是从相同的基础分布中采样的，并且它们已被随机分配给其中一个样本。

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

• k 是 a 中的观察数，

• n 是 a 和 b 中的观察总数，并且

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

数据被汇集(连接)，随机分配给第一个或第二个样本，并计算统计量。该过程重复执行多次，生成零假设下的统计分布。将原始数据的统计量与该分布进行比较以确定 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 的元素的所有可能配对。令 n 为 a 中的观察数，它也必须等于 b 中的观察数。

什么时候1 < n_resamples < factorial(n), 的元素a是随机排列的。用户提供的统计数据接受一个数据参数，比如a_perm, 并计算统计量a_perm和b.这个过程反复进行，排列次，生成零假设下的统计分布。将原始数据的统计量与该分布进行比较以确定 p 值。

什么时候n_resamples >= factorial(n), 进行精确测试：a以每种不同的方式恰好排列一次。因此，统计为每个唯一的样本对计算a和b恰好一次。

对于 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 。令 n 为 a 中的观察数，它也必须等于 b 中的观察数。

什么时候1 < n_resamples < 2**n, 的元素a是b在样本之间随机交换(保持它们的配对)并计算统计量。这个过程反复进行，排列次，生成零假设下的统计分布。将原始数据的统计量与该分布进行比较以确定 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-test，可以只考虑区别在两个配对元素之间。因此，如果data只包含一个样本，则通过独立改变标志每个观察。

警告

p 值是通过对零分布中与统计观察值一样极端或更极端的元素进行计数来计算的。由于使用有限精度算术，当理论值完全相等时，某些统计函数会返回数值上不同的值。在某些情况下，这可能会导致计算的 p 值出现较大误差。 permutation_test 通过将空分布中 “close” (在 1+1e-14 因子内)与检验统计量观测值之间的元素视为等于检验统计量观测值来防止这种情况。但是，建议用户检查零分布以评估这种比较方法是否合适，如果不合适，请手动计算 p 值。请参阅下面的示例。

参考：

[1]

18. 1. 费舍尔。实验设计，第 6 版(1951 年)。

[2]

B. Phipson 和 G. K. Smyth。 “排列 P 值不应该为零：随机抽取排列时计算精确的 P 值。”遗传学和分子生物学中的统计应用 9.1 (2010)。

[3]

M.D.恩斯特。 “排列方法：精确推断的基础”。统计科学(2004 年)。

[4]

B. Efron 和 R. J. Tibshirani。 Bootstrap 简介(1993 年)。

例子：

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

为了提高效率，我们以向量化的方式编写定义测试统计量的函数：样本x和y可以是 ND 数组，并且会为每个 axis-slice 计算统计量轴.

>>> 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.3, random_state=rng)
>>> res = permutation_test((x, y), statistic, n_resamples=100000,
...                        vectorized=True, alternative='less',
...                        random_state=rng)
>>> print(res.statistic)
-0.5230459671240913
>>> print(res.pvalue)
0.00016999830001699983

在原假设下获得小于或等于观测值的检验统计量的近似概率为 0.0225%。这再次低于我们选择的 5% 阈值，因此我们再次有重要证据可以拒绝零假设而支持替代假设。

对于大样本和排列数量，结果与相应的渐近检验的结果相当，独立样本t-test。

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

提供检验统计量的排列分布以供进一步研究。

>>> 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()

如果由于机器精度有限导致统计结果不准确，则必须检查零分布。考虑以下情况：

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

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

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

如果 permutation_test 单纯地执行比较，则值为 0.7999999999999999 的空分布的元素不会被视为与统计观测值一样极端或更极端，因此计算出的 p 值会太小。

>>> incorrect_pvalue = np.count_nonzero(null >= r) / len(null)
>>> incorrect_pvalue
0.1111111111111111  # 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

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