排名比较：两个独立样本

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 中的大多数双样本假设检验都允许“非零”零值，即不要求两个样本中的效应相同的零值。一些假设检验，例如用于列联表的 Fisher 精确检验，以及大多数 k 样本方差分析型检验，只允许零假设所有样本中的效应相同。

样本之间没有差异的假设的零值 p0 取决于效应统计量，差值和相关性的零值为零，比率的零值为一，而随机优越性概率的零值为 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

在等效性检验中，我们想要证明治疗与对照具有近似相同的效果，或者两种治疗具有近似相同的效果。等效性由上下等效边际 (low 和 upp) 定义。如果效果测量位于此区间内，则两种治疗方法是等效的。

零假设和备择假设是

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 月 3 日