等級比較：兩個獨立樣本

Statsmodels 提供了統計量和檢定，用於計算 x1 的值大於 x2 的機率。這些衡量標準是基於使用等級的順序比較。

將 p 定義為第一個樣本母體中隨機抽取的值大於第二個樣本母體中隨機抽取的值的機率，具體來說

p = P(x1 > x2) + 0.5 * P(x1 = x2)

這是 Wilcoxon-Mann-Whitney 的 U 檢定、Fligner-Policello 檢定和 Brunner-Munzel 檢定的基礎衡量標準。推論是基於 Brunner-Munzel 檢定的漸近分佈。對於離散變數，連繫的半機率對應於中位數等級的使用，使其有效。

隨機相等性的虛無假設為 p = 0.5，這對應於 Brunner-Munzel 檢定。

此筆記本簡要概述了 statsmodels 中提供的統計量。

[1]:

import numpy as np

from statsmodels.stats.nonparametric import (
    rank_compare_2indep, rank_compare_2ordinal, prob_larger_continuous,
    cohensd2problarger)

範例

主要函式是 rank_compare_2indep，它計算 Brunner-Munzel 檢定並傳回具有其他方法的 RankCompareResult 實例。

範例的資料取自 Munzel 和 Hauschke 2003，並以頻率計數給出。我們需要將其擴展為觀察陣列，才能將其與 rank_compare_2indep 一起使用。請參閱下面的函式，該函式直接採用頻率計數。標籤或層級被視為順序的，只要它們定義順序 (>、=)，則特定值無關。

[2]:

levels = [-2, -1, 0, 1, 2]
new = [24, 37, 21, 19, 6]
active = [11, 51, 22, 21, 7]

x1 = np.repeat(levels, new)
x2 = np.repeat(levels, active)
np.bincount(x1 + 2), np.bincount(x2 + 2)

[2]:

(array([24, 37, 21, 19,  6]), array([11, 51, 22, 21,  7]))

[3]:

res = rank_compare_2indep(x1, x2) #, use_t=False)
res

[3]:

<class 'statsmodels.stats.nonparametric.RankCompareResult'>
statistic = np.float64(-1.1757561456581607)
pvalue = np.float64(0.24106066495471642)
s1 = np.float64(1164.3327014635863)
s2 = np.float64(701.9050836550837)
var1 = np.float64(0.09281989010392111)
var2 = np.float64(0.06130710836361985)
var = np.float64(0.3098544504968025)
var_prob = np.float64(0.0014148605045516095)
nobs1 = 107
nobs2 = 112
nobs = 219
mean1 = np.float64(105.04672897196262)
mean2 = np.float64(114.73214285714286)
prob1 = np.float64(0.4557743658210948)
prob2 = np.float64(0.5442256341789052)
somersd1 = np.float64(-0.08845126835781036)
somersd2 = np.float64(0.08845126835781048)
df = np.float64(204.29842398679557)
use_t = True
tuple = (np.float64(-1.1757561456581607), np.float64(0.24106066495471642))

結果實例的方法是

[4]:

[i for i in dir(res) if not i.startswith("_") and callable(getattr(res, i))]

[4]:

['conf_int',
 'confint_lintransf',
 'effectsize_normal',
 'summary',
 'test_prob_superior',
 'tost_prob_superior']

[5]:

print(res.summary())  # returns SimpleTable

                  Probability sample 1 is stochastically larger
==================================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
----------------------------------------------------------------------------------
prob(x1>x2) c0     0.4558      0.038     -1.176      0.241       0.382       0.530
==================================================================================

[6]:

ci = res.conf_int()
ci

[6]:

(np.float64(0.3816117144128266), np.float64(0.529937017229363))

[7]:

res.test_prob_superior()

[7]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(-1.1757561456581602)
pvalue = np.float64(0.24106066495471665)
df = np.float64(204.29842398679557)
distribution = 't'
tuple = (np.float64(-1.1757561456581602), np.float64(0.24106066495471665))

[ ]:

單邊檢定、優越性檢定和非劣性檢定

假設檢定函式具有 alternative 關鍵字，用於指定單邊檢定，以及 value 關鍵字，用於指定非零或不相等假設。這兩個關鍵字可以一起用於具有邊際的非劣性檢定或優越性檢定。

非劣性檢定指定一個邊際和替代方案，因此我們可以檢定新處理幾乎與參考處理一樣好或更好的假設。

一般單邊假設是

H0: p = p0
與
HA: p > p0 的替代方案是效果大於指定的虛無值 p0
HA: p < p0 的替代方案是效果小於指定的虛無值 p0

注意：單邊檢定的虛無假設通常指定為弱不等式。然而，p 值通常是在假設虛無假設位於邊界值的情況下得出的。在大多數情況下，邊界值是 p 值的最差情況，虛無假設內部的點通常具有較大的 p 值，並且在這些情況下檢定是保守的。

statsmodels 中的大多數雙樣本假設檢定都允許「非零」虛無值，也就是說，虛無值不要求兩個樣本中的效果相同。一些假設檢定，例如列聯表的費雪精確檢定，以及大多數 k 樣本 anova 型檢定，只允許虛無假設，即所有樣本中的效果相同。

樣本之間沒有差異的假設的虛無值 p0 取決於效果統計量。差異和相關性的虛無值為零，比率的虛無值為 1，而隨機優越機率的虛無值為 0.5。

非劣性和優越性檢定只是允許非零虛無假設的單邊假設檢定的特殊情況。TOST 等效檢定是兩個具有非零虛無值的單邊檢定的組合。

請注意，我們現在使用「優越」的兩種不同含義。檢定中的優越和非劣是指處理效果是否優於控制效果。rank_compare 中的效果衡量標準是基於樣本 1 相對於樣本 2 的隨機優越性的機率，但隨機優越性可能是好的或壞的。

非劣性：較小的值更好

如果具有較低的值更好，例如，如果感興趣的變數是疾病發生率，那麼非劣性檢定會指定一個大於相等性的閾值，其替代方案是參數小於該閾值。

在隨機優越性度量的情況下，兩個樣本的相等性表示機率等於 0.5。如果我們指定一個非劣性邊際（例如 0.6），則虛無假設和替代假設為

H0: p >= 0.6（基於 H0: p = 0.6 的推論的虛無假設）
HA: p < 0.6

如果我們拒絕虛無假設，那麼我們的資料支援該處理非劣於對照組。

非劣性：較大的值更好

在較大的值更好的情況下，例如技能或健康度量，非劣性表示該處理具有幾乎與控制組一樣高或更高的值。這定義了替代假設。假設 p0 是非劣性閾值，例如 p0 = 0.4，則虛無假設和替代假設為

H0: p <= p0 (p0 < 0.5)
HA: p > p0

如果拒絕虛無假設，那麼我們有證據表明該處理（樣本 1）非劣於參考值（樣本 2），也就是說，優越性機率大於 p0。

優越性檢定

優越性檢定通常定義為單邊替代方案，其中相等性為虛無假設，而較好的不等式為替代方案。但是，優越性也可以用邊際定義，在這種情況下，處理必須比邊際指定的值好出一個不可忽略的量。

這表示檢定與上述相同的單邊檢定，其中 p0 等於 0.5，或者如果較大值較好，則 p0 為大於 0.5 的值，如果較小值較好，則 p0 為小於 0.5 的值。

範例：非劣性較小的值更好

假設我們的非劣性閾值為 p0 = 0.55。替代方案為「較小」的單邊檢定之 p 值約為 0.0065，我們在 0.05 的 alpha 值下拒絕虛無假設。資料提供證據表明處理（樣本 1）非劣於控制組（樣本 2），也就是說，我們有證據表明處理最多比控制組差 5 個百分點。

[8]:

res.test_prob_superior(value=0.55, alternative="smaller")

[8]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(-2.505026112598482)
pvalue = np.float64(0.006512753894336686)
df = np.float64(204.29842398679557)
distribution = 't'
tuple = (np.float64(-2.505026112598482), np.float64(0.006512753894336686))

範例：非劣性較大的值更好

現在考慮較大的值更好的情況，且非劣性閾值為 0.45。單邊檢定的 p 值為 0.44，我們無法拒絕虛無假設。因此，我們沒有證據證明該處理最多比控制組差 5 個百分點。

[9]:

res.test_prob_superior(value=0.45, alternative="larger")

[9]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(0.15351382128216023)
pvalue = np.float64(0.43907230923278956)
df = np.float64(204.29842398679557)
distribution = 't'
tuple = (np.float64(0.15351382128216023), np.float64(0.43907230923278956))

等效檢定 TOST

在等效檢定中，我們想要證明該處理與控制組具有大致相同的效果，或者兩個處理具有大致相同的效果。等效性由下限和上限等效邊際（下限和上限）定義。如果效果衡量標準在此區間內，則處理是等效的。

虛無假設和替代假設為

H0: p <= p_low 或 p >= p_upp
HA: p_low < p < p_upp

在此情況下，虛無假設是效果衡量標準在等效區間之外。如果我們拒絕虛無假設，則資料提供證據表明處理和控制組是等效的。

在 TOST 流程中，我們評估兩個單邊檢定，一個用於虛無假設，即效果等於或低於下限閾值，另一個用於虛無假設，即效果等於或高於上限閾值。如果我們拒絕這兩個檢定，則資料提供證據表明效果在等效區間內。tost 方法的 p 值將是兩個檢定的 p 值的最大值。

假設我們的等效邊際為 0.45 和 0.55，則等效檢定的 p 值為 0.43，且我們無法拒絕兩個處理不等效的虛無假設，也就是說，資料不提供對等效性的支持。

[10]:

res.tost_prob_superior(low=0.45, upp=0.55)

[10]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(0.15351382128216023)
pvalue = np.float64(0.43907230923278956)
results_larger = <class 'statsmodels.stats.base.HolderTuple'>
    statistic = np.float64(0.15351382128216023)
    pvalue = np.float64(0.43907230923278956)
    df = np.float64(204.29842398679557)
    distribution = 't'
    tuple = (np.float64(0.15351382128216023), np.float64(0.43907230923278956))
results_smaller = <class 'statsmodels.stats.base.HolderTuple'>
    statistic = np.float64(-2.505026112598482)
    pvalue = np.float64(0.006512753894336686)
    df = np.float64(204.29842398679557)
    distribution = 't'
    tuple = (np.float64(-2.505026112598482), np.float64(0.006512753894336686))
title = 'Equivalence test for Prob(x1 > x2) + 0.5 Prob(x1 = x2) '
tuple = (np.float64(0.15351382128216023), np.float64(0.43907230923278956))

[11]:

res.test_prob_superior(value=0.55, alternative="smaller")

[11]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(-2.505026112598482)
pvalue = np.float64(0.006512753894336686)
df = np.float64(204.29842398679557)
distribution = 't'
tuple = (np.float64(-2.505026112598482), np.float64(0.006512753894336686))

[12]:

res.test_prob_superior(value=0.6, alternative="larger")

[12]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(-3.834296079538801)
pvalue = np.float64(0.9999161566635063)
df = np.float64(204.29842398679557)
distribution = 't'
tuple = (np.float64(-3.834296079538801), np.float64(0.9999161566635063))

[13]:

res.test_prob_superior(value=0.529937, alternative="smaller")

[13]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(-1.9716432456640076)
pvalue = np.float64(0.02500002636314127)
df = np.float64(204.29842398679557)
distribution = 't'
tuple = (np.float64(-1.9716432456640076), np.float64(0.02500002636314127))

[14]:

res.test_prob_superior(value=0.529937017, alternative="smaller")

[14]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(-1.9716436976157954)
pvalue = np.float64(0.025000000351072277)
df = np.float64(204.29842398679557)
distribution = 't'
tuple = (np.float64(-1.9716436976157954), np.float64(0.025000000351072277))

[15]:

res.conf_int()

[15]:

(np.float64(0.3816117144128266), np.float64(0.529937017229363))

旁注：假設檢定中的舉證責任

假設檢定通常具有以下兩個屬性

  • 小樣本有利於虛無假設。小樣本的不確定性很大，且假設檢定的功效很小。在這種情況下，研究的效力不足，且我們通常因為不確定性大而無法拒絕虛無假設。不拒絕不能被視為支持虛無的證據，因為假設檢定沒有太多功效來拒絕虛無假設。

  • 大樣本傾向於支持對立假設。在大樣本中，不確定性會變小，即使與虛無假設的微小偏差也會導致拒絕虛無假設。在這種情況下，我們可能會得到一個顯示統計上顯著但實質上無關的結果。

非劣性檢定和等效性檢定解決了這兩個問題。第一個問題，即在小樣本中傾向於支持虛無假設，可以透過反轉虛無假設和對立假設來避免。對立假設應該是我們想要證明的假設，因此檢定不應僅僅因為統計檢定力太小而支持感興趣的假設。第二個問題，即在大樣本中傾向於支持對立假設，可以透過在假設檢定中指定一個邊際值來避免，如此一來，微不足道的偏差仍然屬於虛無假設的一部分。透過使用這種方法，我們不會因為與點虛無假設無關的偏差而在大樣本中傾向於支持對立假設。非劣性檢定想要證明某種治療方法幾乎與對照組一樣好。等效性檢定試圖證明兩個樣本中的統計數據大致相同，這與指定它們完全相同的點假設相反。

反轉樣本

在文獻中，以機率 p 不等於 0.5 的意義下的隨機優勢，通常使用隨機較小的測量來定義，而不是像 statsmodels 中定義的隨機較大的測量。

效果測量會反轉不等式，因此變為

p2 = P(x1 < x2) + 0.5 * P(x1 = x2)

這可以在 statsmodels 函數中透過反轉樣本的順序來獲得，也就是說，我們可以 使用 (x2, x1) 而不是 (x1, x2)。這兩種定義 p 和 p2 的關係是 p2 = 1 - p。

rank_compare 的結果實例顯示了這兩個統計量，但假設檢定和信賴區間僅針對隨機較大的測量提供。

[16]:

res = rank_compare_2indep(x2, x1) #, use_t=False)
res

[16]:

<class 'statsmodels.stats.nonparametric.RankCompareResult'>
statistic = np.float64(1.1757561456581607)
pvalue = np.float64(0.24106066495471642)
s1 = np.float64(701.9050836550837)
s2 = np.float64(1164.3327014635863)
var1 = np.float64(0.06130710836361985)
var2 = np.float64(0.09281989010392111)
var = np.float64(0.3098544504968025)
var_prob = np.float64(0.0014148605045516095)
nobs1 = 112
nobs2 = 107
nobs = 219
mean1 = np.float64(114.73214285714286)
mean2 = np.float64(105.04672897196262)
prob1 = np.float64(0.5442256341789052)
prob2 = np.float64(0.4557743658210948)
somersd1 = np.float64(0.08845126835781048)
somersd2 = np.float64(-0.08845126835781036)
df = np.float64(204.29842398679557)
use_t = True
tuple = (np.float64(1.1757561456581607), np.float64(0.24106066495471642))

[17]:

res.conf_int()

[17]:

(np.float64(0.470062982770637), np.float64(0.6183882855871734))

[18]:

st = res.summary()
st

[18]:
樣本 1 在機率上較大
係數 標準誤 t P>|t| [0.025 0.975]
prob(x1>x2) c0 0.5442 0.038 1.176 0.241 0.470 0.618

順序數據

基於排序的分析僅假設數據是順序性的，並且僅使用順序信息。此外，根據中位數排名所下的定義，允許在類似於 Brunner-Munzel 檢定的分析中使用離散數據。

當數據是離散的並且只有有限數量的支持點，即屬於有序類別時，Statsmodels 提供了一些額外的函數。

上面的數據是以順序數據的形式給出的，但我們必須擴展它才能使用 rank_compare_2indep。相反，我們可以將 rank_compare_2ordinal 直接用於頻率計數。後者函數與前者主要的不同之處在於，它針對特殊情況使用了更有效的計算。結果統計量將會相同。

在上述範例中，依據底層數值遞增排序的治療組（「新」）和對照組（「活性」）的計數分別為

[19]:

new, active

[19]:

([24, 37, 21, 19, 6], [11, 51, 22, 21, 7])

[20]:

res_o = rank_compare_2ordinal(new, active)
res_o.summary()

[20]:
樣本 1 在機率上較大
係數 標準誤 t P>|t| [0.025 0.975]
prob(x1>x2) c0 0.4558 0.038 -1.176 0.241 0.382 0.530 
[21]:

res_o

[21]:

<class 'statsmodels.stats.nonparametric.RankCompareResult'>
statistic = None
pvalue = None
s1 = None
s2 = None
var1 = np.float64(0.0919524144954732)
var2 = np.float64(0.06075972346751607)
var = np.float64(0.3098544504968023)
var_prob = np.float64(0.0014148605045516088)
nobs1 = np.int64(107)
nobs2 = np.int64(112)
nobs = np.int64(219)
mean1 = None
mean2 = None
prob1 = np.float64(0.4557743658210948)
prob2 = np.float64(0.5442256341789052)
somersd1 = np.float64(-0.08845126835781036)
somersd2 = np.float64(0.08845126835781048)
df = np.float64(204.29842398679557)
use_t = True
tuple = (None, None)

[22]:

res = rank_compare_2indep(x1, x2)
res

[22]:

<class 'statsmodels.stats.nonparametric.RankCompareResult'>
statistic = np.float64(-1.1757561456581607)
pvalue = np.float64(0.24106066495471642)
s1 = np.float64(1164.3327014635863)
s2 = np.float64(701.9050836550837)
var1 = np.float64(0.09281989010392111)
var2 = np.float64(0.06130710836361985)
var = np.float64(0.3098544504968025)
var_prob = np.float64(0.0014148605045516095)
nobs1 = 107
nobs2 = 112
nobs = 219
mean1 = np.float64(105.04672897196262)
mean2 = np.float64(114.73214285714286)
prob1 = np.float64(0.4557743658210948)
prob2 = np.float64(0.5442256341789052)
somersd1 = np.float64(-0.08845126835781036)
somersd2 = np.float64(0.08845126835781048)
df = np.float64(204.29842398679557)
use_t = True
tuple = (np.float64(-1.1757561456581607), np.float64(0.24106066495471642))

[23]:

res_o.test_prob_superior()

[23]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(-1.1757561456581607)
pvalue = np.float64(0.24106066495471642)
df = np.float64(204.29842398679557)
distribution = 't'
tuple = (np.float64(-1.1757561456581607), np.float64(0.24106066495471642))

[24]:

res.test_prob_superior()

[24]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(-1.1757561456581602)
pvalue = np.float64(0.24106066495471665)
df = np.float64(204.29842398679557)
distribution = 't'
tuple = (np.float64(-1.1757561456581602), np.float64(0.24106066495471665))

[25]:

res_o.conf_int(), res.conf_int()

[25]:

((np.float64(0.38161171441282665), np.float64(0.529937017229363)),
 (np.float64(0.3816117144128266), np.float64(0.529937017229363)))
上次更新：2024年10月03日