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

這種比較方法在大多數實際情況下預計是準確的，但建議用戶通過檢查接近統計觀察值的零分布元素來評估這一點。另外，考慮使用可以使用精確算術計算的統計數據(例如整數統計數據)。