ランク比較：二つの独立した標本¶

Statsmodelsは、x1がx2より大きい値を持つ確率に関する統計量と検定を提供します。これらの測定値は、ランクを使用した順序比較に基づいています。

pを、最初の標本の母集団からのランダムな引き出しが、2番目の標本の母集団からのランダムな引き出しよりも大きな値を持つ確率として定義します。具体的には、

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

これは、ウィルコクソン-マン-ホイットニー（Wilcoxon-Mann-Whitney）の U 検定、フリグナー-ポリチェロ（Fligner-Policello）検定、ブランナー-ムンゼル（Brunner-Munzel）検定に基づく測定値です。推論は、ブランナー-ムンゼル検定の漸近分布に基づいています。結びつきのための半確率は、中間ランクの使用に対応しており、離散変数に対して有効にします。

確率的同等性に対する帰無仮説はp = 0.5であり、これはブランナー-ムンゼル検定に対応します。

このノートブックでは、statsmodelsが提供する統計量の簡単な概要を紹介します。

[1]:

import numpy as np

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

例¶

主な関数は rank_compare_2indep で、これはブランナー-ムンゼル検定を計算し、追加のメソッドを持つ 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], dtype=int64),
 array([11, 51, 22, 21,  7], dtype=int64))

[3]:

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

[3]:

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

(0.3816117144128266, 0.529937017229363)

[7]:

res.test_prob_superior()

[7]:

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

[ ]:

片側検定、優越性検定および非劣性検定¶

仮説検定関数には、片側検定を指定するためのalternativeキーワードと、ゼロでないまたは不等式の仮説を指定するためのvalueキーワードがあります。これらの二つのキーワードを組み合わせることで、非劣性検定や優越性検定をマージン付きで行うことができます。

非劣性検定では、マージンと代替仮説が指定されており、新しい処置法が基準となる処置法とほぼ同等かそれ以上であるという仮説を検定できます。

一般的な片側仮説は次のようになります。

H0: p = p0
対して
HA: p > p0 (代替仮説は効果が指定された帰無仮説の値p0より大きい)
HA: p < p0 (代替仮説は効果が指定された帰無仮説の値p0より小さい)

注意：片側検定の帰無仮説は通常、弱い不等式として指定されます。しかし、p 値は通常、帰無仮説が境界値にあると仮定して導かれます。境界値は多くの場合、p 値にとって最悪のケースであり、帰無仮説の内部の点は通常 p 値が大きくなり、このような場合、検定は保守的です。

統計モデルの多くの2標本仮説検定では「ゼロでない」帰無値を許可しています。つまり、二つの標本で効果が同じである必要はないという帰無値です。一部の仮説検定（例えば、Fisherの正確検定や多くのk標本の分散分析型検定）では、すべての標本で効果が同じであるという帰無仮説のみが許可されます。

標本間に差がないという仮説の帰無値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 で帰無仮説を棄却します。このデータは、処置（標本1）が対照（標本2）よりも非劣性であることを示す証拠を提供しており、すなわち処置が対照より最大で 5 パーセントポイント悪いだけであるという証拠があります。

[8]:

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

[8]:

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

[11]:

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

[11]:

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

[12]:

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

[12]:

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

[13]:

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

[13]:

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

[14]:

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

[14]:

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

[15]:

res.conf_int()

[15]:

(0.3816117144128266, 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 = 1.1757561456581607
pvalue = 0.24106066495471642
s1 = 701.9050836550837
s2 = 1164.3327014635863
var1 = 0.06130710836361985
var2 = 0.09281989010392111
var = 0.3098544504968025
var_prob = 0.0014148605045516095
nobs1 = 112
nobs2 = 107
nobs = 219
mean1 = 114.73214285714286
mean2 = 105.04672897196262
prob1 = 0.5442256341789052
prob2 = 0.4557743658210948
somersd1 = 0.08845126835781048
somersd2 = -0.08845126835781036
df = 204.29842398679557
use_t = True
tuple = (1.1757561456581607, 0.24106066495471642)

[17]:

res.conf_int()

[17]:

(0.470062982770637, 0.6183882855871734)

[18]:

st = res.summary()
st

[18]:

Probability sample 1 is stochastically larger
	coef	std err	t	P>\|t\|	[0.025	0.975]
prob(x1>x2) c0	0.5442	0.038	1.176	0.241	0.470	0.618

順序データ¶

ランクベースの分析は、データが順序尺度であることのみを前提とし、順序情報のみを使用します。さらに、中央値を基にした定義は、ブランナー・ムンゼル検定に似た分析で離散データを扱うことを可能にします。

Statsmodelsは、データが離散的でサポート点が有限の順序カテゴリカルデータの場合に使用できる追加の関数を提供しています。

上記のデータは順序データとして与えられましたが、 rank_compare_2indep を使用するために展開する必要がありました。代わりに、頻度カウントを使って直接 rank_compare_2ordinal を使用できます。後者の関数は、特殊なケースに対してより効率的な計算を使用する点で前者と主に異なります。結果の統計は同じになります。

上記の例から、処置（」new」）と対照群（」active」）のカウントは、基礎となる値の増加順に並べると次のようになります。

[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]:

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

[21]:

res_o

[21]:

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

[22]:

res = rank_compare_2indep(x1, x2)
res

[22]:

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

[23]:

res_o.test_prob_superior()

[23]:

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

[24]:

res.test_prob_superior()

[24]:

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

[25]:

res_o.conf_int(), res.conf_int()

[25]:

((0.38161171441282665, 0.529937017229363),
 (0.3816117144128266, 0.529937017229363))

最終更新日: 2025年01月28日