一標本および二標本ポアソン率の統計と推論¶

著者: Josef Perktold

このノートブックは、一標本および二標本の場合におけるポアソン率の仮説検定、信頼区間、その他の統計の概要を簡単に説明します。詳細なオプションと追加情報は、docstringを参照してください。

statsmodels.stats.rates内のすべての関数は、データのサマリー統計を引数として受け取ります。それらは、イベントの数と観測数または総曝露数です。ポアソンに関するいくつかの関数には、過剰分散のオプションがあります。負の二項分布（NB2）に関する関数は、分散パラメータが必要です。過剰分散と分散パラメータはユーザーが提供する必要があり、それぞれGLM-ポアソンと離散負の二項モデルを使用して元のデータから推定できます。

一部はまだ実験段階であり、変更される可能性があります。一部の機能はまだ欠けており、今後のバージョンで追加される予定です。

一標本関数

二標本関数

[1]:

import numpy as np
from numpy.testing import assert_allclose
import statsmodels.stats.rates as smr
from statsmodels.stats.rates import (
    # functions for 1 sample
    test_poisson,
    confint_poisson,
    tolerance_int_poisson,
    confint_quantile_poisson,

    # functions for 2 sample
    test_poisson_2indep,
    etest_poisson_2indep,
    confint_poisson_2indep,
    tost_poisson_2indep,
    nonequivalence_poisson_2indep,

    # power functions
    power_poisson_ratio_2indep,
    power_poisson_diff_2indep,
    power_equivalence_poisson_2indep,
    power_negbin_ratio_2indep,
    power_equivalence_neginb_2indep,

    # list of statistical methods
    method_names_poisson_1samp,
    method_names_poisson_2indep,
    )

一標本関数¶

現在、一標本ポアソン率の主な関数はtest_poissonとconfint_poissonです。どちらもいくつかの利用可能なメソッドがあり、そのほとんどは仮説検定と信頼区間で一貫しています。許容区間と分位数の信頼区間用に2つの追加関数が利用可能です。

詳細はdocstringを参照してください。

[2]:

count1, n1 = 60, 514.775
count1 / n1

[2]:

0.11655577679568745

[3]:

test_poisson(count1, n1, value=0.1, method="midp-c")

[3]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = nan
pvalue = np.float64(0.23913820865664664)
distribution = 'Poisson'
method = 'midp-c'
alternative = 'two-sided'
rate = 0.11655577679568745
nobs = 514.775
tuple = (nan, np.float64(0.23913820865664664))

[4]:

confint_poisson(count1, n1, method="midp-c")

[4]:

(np.float64(0.0897357524941493), np.float64(0.1490015282355224))

仮説検定と信頼区間の利用可能なメソッドは、辞書method_names_poisson_1sampで確認できます。詳細はdocstringを参照してください。

[5]:

method_names_poisson_1samp

[5]:

{'test': ['wald',
  'score',
  'exact-c',
  'midp-c',
  'waldccv',
  'sqrt-a',
  'sqrt-v',
  'sqrt'],
 'confint': ['wald',
  'score',
  'exact-c',
  'midp-c',
  'jeff',
  'waldccv',
  'sqrt-a',
  'sqrt-v',
  'sqrt',
  'sqrt-cent',
  'sqrt-centcc']}

[6]:

for meth in method_names_poisson_1samp["test"]:
    tst = test_poisson(count1, n1, method=meth, value=0.1,
                       alternative='two-sided')
    print("%-12s" % meth, tst.pvalue)

wald         0.2712232025335152
score        0.23489608509894766
exact-c      0.2654698417416039
midp-c       0.23913820865664664
waldccv      0.27321266612309003
sqrt-a       0.25489746088635834
sqrt-v       0.2281700763432699
sqrt         0.2533006997208508

[7]:

for meth in method_names_poisson_1samp["confint"]:
    tst = confint_poisson(count1, n1, method=meth)
    print("%-12s" % meth, tst)

wald         (np.float64(0.08706363801159746), np.float64(0.14604791557977745))
score        (np.float64(0.0905597500576385), np.float64(0.15001420714831387))
exact-c      (np.float64(0.08894433674907924), np.float64(0.15003038882355074))
midp-c       (np.float64(0.0897357524941493), np.float64(0.1490015282355224))
jeff         (np.float64(0.08979284758964944), np.float64(0.14893677466593855))
waldccv      (np.float64(0.08694100904696915), np.float64(0.14617054454440576))
sqrt-a       (np.float64(0.08883721953786133), np.float64(0.14800553586080228))
sqrt-v       (np.float64(0.08975547672311084), np.float64(0.14897854470462502))
sqrt         (np.float64(0.08892923891524183), np.float64(0.14791351648342183))
sqrt-cent    (np.float64(0.08883721953786133), np.float64(0.1480055358608023))
sqrt-centcc  (np.float64(0.0879886777703761), np.float64(0.1490990831089978))

現在、一標本ポアソン率用に2つの追加関数が利用可能です。tolerance_int_poissonは許容区間用で、confint_quantile_poissonはポアソン分位数の信頼区間用です。

許容区間は予測区間と似ており、新しい観測値のランダム性と推定されたポアソン率の不確実性を組み合わせたものです。率が既知であれば、所定の率で逆累積分布関数を使用して、新しい観測値のポアソン区間を計算できます。許容区間は、率の信頼区間を使用して、率の不確実性を加えます。

許容区間は2つの確率で指定されます。probはポアソン区間のカバレッジで、alphaは率推定値の信頼区間に対する信頼水準です。

確率は、カウントが離散確率変数であるため、名目上の確率に正確に等しくなることはできません。区間の特性は不等式で指定されます。カバレッジは少なくともprobであり、推定率の信頼区間のカバレッジは少なくとも1 - alphaです。ただし、分布が正しく指定されている場合でも、ほとんどのメソッドは、小さなサンプルではカバレッジの不等式が成り立つことを保証しません。

次の例では、総曝露数または観測数が100の場合、所定のカバレッジprobと信頼水準alphaで、4〜23のイベントが発生すると予想できます。許容区間は、観測された率（5、19）におけるポアソン区間よりも大きくなります。これは、許容区間がパラメータ推定値の不確実性を考慮に入れるためです。

[8]:

exposure_new = 100
tolerance_int_poisson(count1, n1, prob=0.95, exposure_new=exposure_new, method="score", alpha=0.05, alternative='two-sided')

[8]:

(np.float64(4.0), np.float64(23.0))

[9]:

from scipy import stats
stats.poisson.interval(0.95, count1 / n1 * exposure_new)

[9]:

(np.float64(5.0), np.float64(19.0))

補足：alpha=1を指定することで、許容区間がパラメータの不確実性を無視するように強制できます。

[10]:

tolerance_int_poisson(count1, n1, prob=0.95, exposure_new=exposure_new, method="score", alpha=1, alternative='two-sided')

[10]:

(np.float64(5.0), np.float64(19.0))

最後の関数は、ポアソン分位数の信頼区間を返します。分位数は累積分布関数（CDF）の逆関数であり、scipy.stats分布ではppfと呼ばれています。

次の例は、CDF確率0.975におけるポアソン区間の上限の信頼区間を示しています。片側カバレッジ確率を使用した上限信頼限界は、許容区間の限界と同じです。

[11]:

confint_quantile_poisson(count1, n1, prob=0.975, exposure_new=100, method="score", alpha=0.05, alternative='two-sided')

[11]:

(np.float64(15.0), np.float64(23.0))

二標本関数¶

二標本の統計関数は、比率または差によって率を比較できます。デフォルトは率の比率を比較することです。

etest関数は、test_poisson_2indepを通じて直接アクセスできます。

[12]:

count1, n1, count2, n2 = 60, 514.775, 30, 543.087

[13]:

test_poisson_2indep(count1, n1, count2, n2, method='etest-score')

[13]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(3.4174018390002145)
pvalue = np.float64(0.0005672617581628009)
distribution = 'poisson'
compare = 'ratio'
method = 'etest-score'
alternative = 'two-sided'
rates = (np.float64(0.11655577679568745), np.float64(0.055239768213932575))
ratio = np.float64(2.10999757175465)
diff = np.float64(0.06131600858175487)
value = 1
rates_cmle = None
ratio_null = 1
tuple = (np.float64(3.4174018390002145), np.float64(0.0005672617581628009))

[14]:

confint_poisson_2indep(count1, n1, count2, n2, method='score',
                       compare="ratio")

[14]:

(np.float64(1.3659624311981189), np.float64(3.2593061483872257))

[15]:

confint_poisson_2indep(count1, n1, count2, n2, method='score',
                       compare="diff")

[15]:

(np.float64(0.026579645509259224), np.float64(0.0989192191413259))

二標本検定関数であるtest_poisson_2indepには、等式を指定しない帰無仮説を指定するvalueオプションがあります。これは、片側検定を用いた優越性検定と非劣性検定に役立ちます。

例として、次の検定は、率の比率が2であるという両側帰無仮説を検定します。この仮説のp値は0.81であり、最初の率が2番目の率の2倍であるという帰無仮説を棄却できません。

[16]:

test_poisson_2indep(count1, n1, count2, n2, value=2, method='etest-score')

[16]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(0.23946504079843253)
pvalue = np.float64(0.813504857205675)
distribution = 'poisson'
compare = 'ratio'
method = 'etest-score'
alternative = 'two-sided'
rates = (np.float64(0.11655577679568745), np.float64(0.055239768213932575))
ratio = np.float64(2.10999757175465)
diff = np.float64(0.06131600858175487)
value = 2
rates_cmle = None
ratio_null = 2
tuple = (np.float64(0.23946504079843253), np.float64(0.813504857205675))

method_names_poisson_2indep辞書は、率の比率または率の差によって2つのサンプルを比較する場合に利用可能なメソッドを示しています。

この辞書を使用して、利用可能なすべてのメソッドでp値と信頼区間を計算できます。

[17]:

method_names_poisson_2indep

[17]:

{'test': {'ratio': ['wald',
   'score',
   'score-log',
   'wald-log',
   'exact-cond',
   'cond-midp',
   'sqrt',
   'etest-score',
   'etest-wald'],
  'diff': ['wald', 'score', 'waldccv', 'etest-score', 'etest-wald']},
 'confint': {'ratio': ['waldcc',
   'score',
   'score-log',
   'wald-log',
   'sqrtcc',
   'mover'],
  'diff': ['wald', 'score', 'waldccv', 'mover']}}

[18]:

for compare in ["ratio", "diff"]:
    print(compare)
    for meth in method_names_poisson_2indep["test"][compare]:
        tst = test_poisson_2indep(count1, n1, count2, n2, value=None,
                                  method=meth, compare=compare,
                                  alternative='two-sided')
        print("   %-12s" % meth, tst.pvalue)

ratio
   wald         0.0007120093285061108
   score        0.0006322188820470972
   score-log    0.0003992519661848979
   wald-log     0.0008399438093390379
   exact-cond   0.0006751826586863219
   cond-midp    0.0005572624066190538
   sqrt         0.0005700355621795108
   etest-score  0.0005672617581628009
   etest-wald   0.0006431446124897875
diff
   wald         0.0007120093285061094
   score        0.0006322188820470944
   waldccv      0.0007610462660136599
   etest-score  0.000567261758162795
   etest-wald   0.0006431446124897808

同様に、現在利用可能なすべてのメソッドについて、率の比率と率の差の信頼区間を計算できます。

[19]:

for compare in ["ratio", "diff"]:
    print(compare)
    for meth in method_names_poisson_2indep["confint"][compare]:
        ci = confint_poisson_2indep(count1, n1, count2, n2,
                                  method=meth, compare=compare)
        print("   %-12s" % meth, ci)

ratio
   waldcc       (np.float64(1.354190544703406), np.float64(3.233964238781885))
   score        (np.float64(1.3659624311981189), np.float64(3.2593061483872257))
   score-log    (np.float64(1.3903411228996467), np.float64(3.4348249508085043))
   wald-log     (np.float64(1.3612801263025065), np.float64(3.2705169691290763))
   sqrtcc       (np.float64(1.29635711135392), np.float64(3.132234781692197))
   mover        (np.float64(1.3614682485833316), np.float64(3.258622814678696))
diff
   wald         (np.float64(0.02581223514639487), np.float64(0.09681978201711487))
   score        (np.float64(0.026579645509259224), np.float64(0.0989192191413259))
   waldccv      (np.float64(0.025618973109117968), np.float64(0.09701304405439178))
   mover        (np.float64(0.026193641039269785), np.float64(0.09864127183950336))

区間仮説を指定する仮説検定のための2つの追加関数、tost_poisson_2indepとnonequivalence_poisson_2indepがあります。

TOST関数は、代替仮説が2つの率がお互いの区間内にあると指定する同等性検定を実装します。

nonequivalence検定は、代替仮説が2つの率が所定のゼロ以外の値以上だけ異なることを指定する検定を実装します。これは最小効果検定とも呼ばれます。この検定は、TOSTと同様に2つの片側検定を使用しますが、帰無仮説と代替仮説は同等性検定とは逆になります。

両方の関数はtest_poisson_2indepに委任するため、同じmethodオプションが利用可能です。

次の同等性検定は、率の比率が0.8と1/0.8の間にあるという代替仮説を指定します。観測された率の比率は0.89です。p値は0.107であり、所定のマージンで2つの率が同等であるという代替仮説を支持して帰無仮説を棄却できません。したがって、仮説検定は、2つの率が同等であるという証拠を提供しません。

2番目の例では、同等性がマージン（-0.04、0.04）で定義されている率の差の同等性を検定します。p値は約0.2であり、この検定は2つの率が同等であることを支持しません。

[20]:

low = 0.8
upp = 1 / low

count1, n1, count2, n2 = 200, 1000, 450, 2000

tost_poisson_2indep(count1, n1, count2, n2, low, upp, method='score', compare='ratio')

[20]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(1.2403473458920846)
pvalue = np.float64(0.10742347370282446)
method = 'score'
compare = 'ratio'
equiv_limits = (0.8, 1.25)
results_larger = <class 'statsmodels.stats.base.HolderTuple'>
    statistic = np.float64(1.2403473458920846)
    pvalue = np.float64(0.10742347370282446)
    distribution = 'normal'
    compare = 'ratio'
    method = 'score'
    alternative = 'larger'
    rates = (np.float64(0.2), np.float64(0.225))
    ratio = np.float64(0.888888888888889)
    diff = np.float64(-0.024999999999999994)
    value = 0.8
    rates_cmle = None
    ratio_null = 0.8
    tuple = (np.float64(1.2403473458920846), np.float64(0.10742347370282446))
results_smaller = <class 'statsmodels.stats.base.HolderTuple'>
    statistic = np.float64(-4.0311288741492755)
    pvalue = np.float64(2.7754797240370253e-05)
    distribution = 'normal'
    compare = 'ratio'
    method = 'score'
    alternative = 'smaller'
    rates = (np.float64(0.2), np.float64(0.225))
    ratio = np.float64(0.888888888888889)
    diff = np.float64(-0.024999999999999994)
    value = 1.25
    rates_cmle = None
    ratio_null = 1.25
    tuple = (np.float64(-4.0311288741492755), np.float64(2.7754797240370253e-05))
title = 'Equivalence test for 2 independent Poisson rates'
tuple = (np.float64(1.2403473458920846), np.float64(0.10742347370282446))

[21]:

upp = 0.04
low = -upp
tost_poisson_2indep(count1, n1, count2, n2, low, upp, method='score', compare='diff')

[21]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(0.8575203124598336)
pvalue = np.float64(0.19557869693808477)
method = 'score'
compare = 'diff'
equiv_limits = (-0.04, 0.04)
results_larger = <class 'statsmodels.stats.base.HolderTuple'>
    statistic = np.float64(0.8575203124598336)
    pvalue = np.float64(0.19557869693808477)
    distribution = 'normal'
    compare = 'diff'
    method = 'score'
    alternative = 'larger'
    rates = (np.float64(0.2), np.float64(0.225))
    ratio = np.float64(0.888888888888889)
    diff = np.float64(-0.024999999999999994)
    value = -0.04
    rates_cmle = (np.float64(0.19065363652113884), np.float64(0.23065363652113885))
    ratio_null = None
    tuple = (np.float64(0.8575203124598336), np.float64(0.19557869693808477))
results_smaller = <class 'statsmodels.stats.base.HolderTuple'>
    statistic = np.float64(-3.4807277010355238)
    pvalue = np.float64(0.00025002679047994814)
    distribution = 'normal'
    compare = 'diff'
    method = 'score'
    alternative = 'smaller'
    rates = (np.float64(0.2), np.float64(0.225))
    ratio = np.float64(0.888888888888889)
    diff = np.float64(-0.024999999999999994)
    value = 0.04
    rates_cmle = (np.float64(0.24581855699051405), np.float64(0.20581855699051405))
    ratio_null = None
    tuple = (np.float64(-3.4807277010355238), np.float64(0.00025002679047994814))
title = 'Equivalence test for 2 independent Poisson rates'
tuple = (np.float64(0.8575203124598336), np.float64(0.19557869693808477))

関数nonequivalence_poisson_2indepは、2つの率が無視できない量だけ異なるという代替仮説を検定します。

次の例では、代替仮説は、率の比率が区間（0.95、1/0.95）の外側にあることを指定します。帰無仮説は、比率の比率が区間内にあるということです。検定が帰無仮説を棄却した場合、それは率の比率が区間限界で指定された無視できる量よりも大きく異なるという証拠を提供します。

大サンプルにおける点帰無仮説検定と区間帰無仮説検定の関係に関する注記。帰無仮説が正確に成り立たず、サンプルサイズが十分に大きい場合、test_poisson_2indepの点帰無仮説は、帰無仮説からわずかなずれを棄却します。率が指定された無視できる量以下しか異ならない場合、非同等性検定または最小効果検定は大サンプル（サンプルが無限に近づく）では帰無仮説を棄却しません。

この例では、点帰無仮説と区間帰無仮説のどちらも棄却されません。率が統計的に異なるという十分な証拠はありません。それに続いて、観測された率を一定に保ちながら、サンプルサイズを20倍に増やします。この場合、点帰無仮説検定は棄却され、p値は0.01ですが、区間帰無仮説は棄却されず、p値は1になります。

注：非同等性検定は一般的に保守的であり、そのサイズはαによって制限されますが、非同等性マージンを固定した大サンプルの限界では、そのサイズはα/2に近づきます。非同等性区間が1点に縮小すると、非同等性検定は点仮説検定と同じになります。（docstringを参照）

[22]:

count1, n1, count2, n2 = 200, 1000, 420, 2000
low = 0.95
upp = 1 / low
nf = 1
nonequivalence_poisson_2indep(count1 * nf, n1 * nf, count2 * nf, n2 * nf, low, upp, method='score', compare='ratio')

[22]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(-1.1654330934961301)
pvalue = np.float64(1.0232437381644721)
method = 'score'
results_larger = <class 'statsmodels.stats.base.HolderTuple'>
    statistic = np.float64(0.02913582733740325)
    pvalue = np.float64(0.5116218690822361)
    distribution = 'normal'
    compare = 'ratio'
    method = 'score'
    alternative = 'smaller'
    rates = (np.float64(0.2), np.float64(0.21))
    ratio = np.float64(0.9523809523809524)
    diff = np.float64(-0.009999999999999981)
    value = 0.95
    rates_cmle = None
    ratio_null = 0.95
    tuple = (np.float64(0.02913582733740325), np.float64(0.5116218690822361))
results_smaller = <class 'statsmodels.stats.base.HolderTuple'>
    statistic = np.float64(-1.1654330934961301)
    pvalue = np.float64(0.8780781359377093)
    distribution = 'normal'
    compare = 'ratio'
    method = 'score'
    alternative = 'larger'
    rates = (np.float64(0.2), np.float64(0.21))
    ratio = np.float64(0.9523809523809524)
    diff = np.float64(-0.009999999999999981)
    value = 1.0526315789473684
    rates_cmle = None
    ratio_null = 1.0526315789473684
    tuple = (np.float64(-1.1654330934961301), np.float64(0.8780781359377093))
title = 'Equivalence test for 2 independent Poisson rates'
tuple = (np.float64(-1.1654330934961301), np.float64(1.0232437381644721))

[23]:

test_poisson_2indep(count1 * nf, n1 * nf, count2 * nf, n2 * nf, method='score', compare='ratio')

[23]:

<class 'statsmodels.stats.base.HolderTuple'>
statistic = np.float64(-0.5679618342470648)
pvalue = np.float64(0.5700608835629815)
distribution = 'normal'
compare = 'ratio'
method = 'score'
alternative = 'two-sided'
rates = (np.float64(0.2), np.float64(0.21))
ratio = np.float64(0.9523809523809524)
diff = np.float64(-0.009999999999999981)
value = 1
rates_cmle = None
ratio_null = 1
tuple = (np.float64(-0.5679618342470648), np.float64(0.5700608835629815))

[24]:

nf = 20
nonequivalence_poisson_2indep(count1 * nf, n1 * nf, count2 * nf, n2 * nf, low, upp, method='score', compare='ratio').pvalue

[24]:

np.float64(1.1036704302254083)

[25]:

test_poisson_2indep(count1 * nf, n1 * nf, count2 * nf, n2 * nf, method='score', compare='ratio').pvalue

[25]:

np.float64(0.01108516638060269)

検出力¶

Statsmodelsは、2標本ポアソン率と負の二項率の比較のための統計的検出力を計算するためのサポートが限られています。それらは、ZhuとLakkis、およびZhuによる両方の分布の比率比較に基づいており、ポアソン率の差の基本的な正規分布に基づく比較に基づいています。特にGuなど、仮説検定関数で利用可能なメソッドにより密接に対応する他のメソッドはまだ利用できません。

利用可能な関数は次のとおりです。

[26]:

power_poisson_ratio_2indep
power_equivalence_poisson_2indep
power_negbin_ratio_2indep
power_equivalence_neginb_2indep

power_poisson_diff_2indep

[26]:

<function statsmodels.stats.rates.power_poisson_diff_2indep(rate1, rate2, nobs1, nobs_ratio=1, alpha=0.05, value=0, method_var='score', alternative='two-sided', return_results=True)>

最終更新日：2024年10月3日