順位比較：独立した2標本¶

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

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

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

これは、ウィルコクソン・マン・ホイットニーのU検定、フリグナー・ポリチェロ検定、ブルーナー・ムンツェル検定の基礎となる尺度です。推論は、ブルーナー・ムンツェル検定の漸近分布に基づいています。タイに対する半分の確率は、ミッドランクの使用に対応し、離散変数に対しても有効になります。

確率的等価性に対する帰無仮説は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]), 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のほとんどの2標本仮説検定は、「ゼロ以外」の帰無値、つまり2つの標本における効果が同じであることを要求しない帰無値を許容します。フィッシャーの正確確率検定（分割表用）やほとんどのk標本ANOVA型検定など、一部の仮説検定では、効果がすべての標本で同じであるという帰無仮説しか許容しません。

標本間に差がないという仮説に対する帰無値p0は、効果統計量に依存します。差と相関関係に対する帰無値はゼロであり、比率に対する帰無値は1であり、確率的優越確率に対する帰無値は0.5です。

非劣性検定と優越性検定は、ゼロ以外の帰無仮説を許容する片側仮説検定の特別なケースにすぎません。TOST等価性検定は、ゼロ以外の帰無値を持つ2つの片側検定の組み合わせです。

注：「優越」という言葉を2つの異なる意味で使用しています。検定における優越と非劣性は、治療効果がコントロールよりも優れているかどうかを示します。rank_compareにおける効果尺度は、標本1が標本2よりも確率的に優れていることに基づく確率ですが、確率的優越性は良い場合も悪い場合もあります。

非劣性：値が小さい方が良い

例えば、関心のある変数が疾患発生率である場合、値が小さい方が良い場合、非劣性検定は、パラメータが閾値より小さいという代替案を持つ等価性よりも大きい閾値を指定します。

確率的優越性尺度の場合、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に等しいか、値が大きい方が良い場合は0.5より大きく、値が小さい方が良い場合は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 = 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¶

等価性検定では、治療がコントロールとほぼ同じ効果を持っていること、または2つの治療がほぼ同じ効果を持っていることを示すことを目的としています。等価であることは、下限と上限の等価性マージン（lowとupp）によって定義されます。効果尺度がこの範囲内にある場合、治療は等価です。

帰無仮説と対立仮説は

H0: p <= p_low または p >= p_upp

HA: p_low < p < p_upp

この場合、帰無仮説は、効果尺度が等価性範囲の外にあることです。帰無仮説を棄却した場合、データは治療とコントロールが等価であるという証拠を提供します。

TOST手順では、効果が下限閾値以下であるという帰無仮説と、効果が上限閾値以上であるという帰無仮説の2つの片側検定を評価します。両方の検定を棄却した場合、データは効果が等価性範囲内にあるという証拠を提供します。tostメソッドのp値は、2つの検定のp値の最大値になります。

等価性マージンが0.45と0.55であると仮定すると、等価性検定のp値は0.43であり、2つの治療が等価ではないという帰無仮説を棄却できません。つまり、データは等価性を支持する証拠を提供しません。

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

補足：仮説検定における証明責任¶

仮説検定は一般的に次の2つの特性を持っています。

小標本は帰無仮説を支持します。小標本における不確実性は大きく、仮説検定の検出力は小さいです。この場合の研究は検出力不足であり、大きな不確実さのために帰無仮説を棄却できないことがよくあります。仮説検定には帰無仮説を棄却する能力があまりないため、棄却できないことを帰無仮説を支持する証拠とはみなせません。
大標本は対立仮説を支持します。大標本では不確実性が小さくなり、帰無仮説からのわずかなずれでも棄却につながります。この場合、統計的に有意だが実質的に無関係な効果を示す検定結果を得る可能性があります。

非劣性試験と同等性試験は、これら両方の問題を解決します。最初の問題である、小さなサンプルにおける帰無仮説の支持は、帰無仮説と対立仮説を逆にすることで回避できます。対立仮説は証明したい仮説であるべきであり、検出力があまりにも小さいために単に関心のある仮説を支持するだけにならないようにする必要があります。2番目の問題である、大きなサンプルにおける対立仮説の支持は、仮説検定において限界値を指定することで回避できます。これにより、些細なずれも依然として帰無仮説の一部となります。これにより、点推定の無関係なずれのために、大きなサンプルにおいて対立仮説を支持することはなくなります。非劣性試験は、ある治療法が対照群とほぼ同等に有効であることを示そうとします。同等性試験は、点推定のように完全に同じであると規定するのではなく、2つのサンプルの統計量がほぼ同じであることを示そうとします。

サンプルの反転¶

文献において、p≠0.5という意味での確率優位性は、statsmodelsで定義されている確率的に大きい代わりに、確率的に小さい尺度を用いて定義されることが多いです。

効果尺度は不等式を逆転させ、以下のようになります。

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

これは、statsmodelsの関数においてサンプルの順序を逆にすることで得られます。(x1, x2)の代わりに(x2, x1)を使用できます。2つの定義、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日