1つおよび2つの標本におけるポアソン率の統計と推論¶

著者: Josef Perktold

このノートブックでは、1つおよび2つの標本の場合におけるポアソン率に関する仮説検定、信頼区間、その他の統計の概要を簡単に説明します。詳細や追加のオプションについては、ドキュメンテーション文字列を参照ください。

statsmodels.stats.rates 内のすべての関数は、データの要約統計を引数として受け取ります。それらは、イベントのカウント数と観測数、または合計暴露量です。ポアソンに関する一部の関数は、過剰分散のオプションを持っています。ネガティブ二項分布（NB2）の関数は、分散パラメータが必要です。過剰分散や分散パラメータはユーザーが指定する必要があり、それぞれ GLM-Poisson モデルや離散負二項モデルを使用して元データから推定することができます。

なお、一部の機能はまだ実験的であり、変更される可能性があります。また、いくつかの機能はまだ未実装であり、将来のバージョンで追加される予定です。

一標本関数

二標本関数

[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つの追加関数も利用可能です。

詳細についてはドキュメンテーション文字列（docstrings）を参照してください。

[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 = 0.23913820865664664
distribution = 'Poisson'
method = 'midp-c'
alternative = 'two-sided'
rate = 0.11655577679568745
nobs = 514.775
tuple = (nan, 0.23913820865664664)

[4]:

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

[4]:

(0.0897357524941493, 0.1490015282355224)

仮説検定および信頼区間の利用可能な方法は、辞書method_names_poisson_1sampに含まれています。詳細については、ドックストリングを参照してください。

[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         (0.08706363801159746, 0.14604791557977745)
score        (0.0905597500576385, 0.15001420714831387)
exact-c      (0.08894433674907924, 0.15003038882355074)
midp-c       (0.0897357524941493, 0.1490015282355224)
jeff         (0.08979284758964944, 0.14893677466593855)
waldccv      (0.08694100904696915, 0.14617054454440576)
sqrt-a       (0.08883721953786133, 0.14800553586080228)
sqrt-v       (0.08975547672311084, 0.14897854470462502)
sqrt         (0.08892923891524183, 0.14791351648342183)
sqrt-cent    (0.08883721953786133, 0.1480055358608023)
sqrt-centcc  (0.0879886777703761, 0.1490990831089978)

現在、一標本のポアソン率に対して利用可能な2つの追加機能として、許容区間を求める tolerance_int_poisson と、ポアソン分位数の信頼区間を求める confint_quantile_poisson があります。

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

許容区間は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]:

(4.0, 23.0)

[9]:

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

[9]:

(5.0, 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]:

(5.0, 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]:

(15.0, 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 = 3.4174018390002145
pvalue = 0.0005672617581628009
distribution = 'poisson'
compare = 'ratio'
method = 'etest-score'
alternative = 'two-sided'
rates = (0.11655577679568745, 0.055239768213932575)
ratio = 2.10999757175465
diff = 0.06131600858175487
value = 1
rates_cmle = None
ratio_null = 1
tuple = (3.4174018390002145, 0.0005672617581628009)

[14]:

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

[14]:

(1.3659624311981189, 3.2593061483872257)

[15]:

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

[15]:

(0.026579645509259224, 0.0989192191413259)

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

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

[16]:

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

[16]:

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

method_names_poisson_2indep辞書は、二つの標本をレート比またはレート差で比較する際に使用できる方法を示しています。

この辞書を使用して、利用可能なすべての方法を用いて 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.0006751826586863236
   cond-midp    0.0005572624066190556
   sqrt         0.0005700355621795109
   etest-score  0.0005672617581628009
   etest-wald   0.0006431446124897875
diff
   wald         0.0007120093285061096
   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       (1.354190544703406, 3.233964238781885)
   score        (1.3659624311981189, 3.2593061483872257)
   score-log    (1.3903411228996467, 3.4348249508085043)
   wald-log     (1.3612801263025065, 3.2705169691290763)
   sqrtcc       (1.29635711135392, 3.132234781692197)
   mover        (1.3614682485833316, 3.258622814678696)
diff
   wald         (0.02581223514639487, 0.09681978201711487)
   score        (0.026579645509259224, 0.0989192191413259)
   waldccv      (0.025618973109117968, 0.09701304405439178)
   mover        (0.026193641039269785, 0.09864127183950336)

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

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

nonequivalence 検定は、代替仮説が二つの率が少なくとも指定された非ゼロの値だけ異なることを指定する検定です。これは最小効果検定とも呼ばれます。この検定は、 TOST と似た二つの片側検定を使用しますが、等価性検定と比較して帰無仮説と代替仮説が逆転しています。

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

次の同等性検定では、代替仮説として比率が 0.8 と 1/0.8 の間であることを指定しています。観測された比率は 0.89 です。 p 値は 0.107 であり、与えられた限界において二つの比率が同等であるという代替仮説を支持するために帰無仮説を棄却することはできません。したがって、この仮説検定は、二つの比率が同等であるという証拠を提供しません。

2番目の例では、率の差の同等性を検定しています。同等性は限界（-0.04, 0.04）によって定義されます。 p 値は約 0.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 = 1.2403473458920846
pvalue = 0.10742347370282446
method = 'score'
compare = 'ratio'
equiv_limits = (0.8, 1.25)
results_larger = <class 'statsmodels.stats.base.HolderTuple'>
    statistic = 1.2403473458920846
    pvalue = 0.10742347370282446
    distribution = 'normal'
    compare = 'ratio'
    method = 'score'
    alternative = 'larger'
    rates = (0.2, 0.225)
    ratio = 0.888888888888889
    diff = -0.024999999999999994
    value = 0.8
    rates_cmle = None
    ratio_null = 0.8
    tuple = (1.2403473458920846, 0.10742347370282446)
results_smaller = <class 'statsmodels.stats.base.HolderTuple'>
    statistic = -4.0311288741492755
    pvalue = 2.7754797240370253e-05
    distribution = 'normal'
    compare = 'ratio'
    method = 'score'
    alternative = 'smaller'
    rates = (0.2, 0.225)
    ratio = 0.888888888888889
    diff = -0.024999999999999994
    value = 1.25
    rates_cmle = None
    ratio_null = 1.25
    tuple = (-4.0311288741492755, 2.7754797240370253e-05)
title = 'Equivalence test for 2 independent Poisson rates'
tuple = (1.2403473458920846, 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 = 0.8575203124598336
pvalue = 0.19557869693808477
method = 'score'
compare = 'diff'
equiv_limits = (-0.04, 0.04)
results_larger = <class 'statsmodels.stats.base.HolderTuple'>
    statistic = 0.8575203124598336
    pvalue = 0.19557869693808477
    distribution = 'normal'
    compare = 'diff'
    method = 'score'
    alternative = 'larger'
    rates = (0.2, 0.225)
    ratio = 0.888888888888889
    diff = -0.024999999999999994
    value = -0.04
    rates_cmle = (0.19065363652113884, 0.23065363652113885)
    ratio_null = None
    tuple = (0.8575203124598336, 0.19557869693808477)
results_smaller = <class 'statsmodels.stats.base.HolderTuple'>
    statistic = -3.4807277010355238
    pvalue = 0.00025002679047994814
    distribution = 'normal'
    compare = 'diff'
    method = 'score'
    alternative = 'smaller'
    rates = (0.2, 0.225)
    ratio = 0.888888888888889
    diff = -0.024999999999999994
    value = 0.04
    rates_cmle = (0.24581855699051405, 0.20581855699051405)
    ratio_null = None
    tuple = (-3.4807277010355238, 0.00025002679047994814)
title = 'Equivalence test for 2 independent Poisson rates'
tuple = (0.8575203124598336, 0.19557869693808477)

関数nonequivalence_poisson_2indepは、二つの比率が無視できない程度に異なるという対立仮説を検定します。

次の例では、対立仮説は比率が区間 (0.95, 1/0.95) の外にあることを指定しています。帰無仮説は比率がその区間内にあることです。もし帰無仮説が棄却されれば、比率が区間の限界で指定された無視できない量以上に異なるという証拠が得られたことになります。

大きな標本における点仮説検定と区間仮説検定の関係についての注記です。 test_poisson_2indep の帰無仮説は、もし帰無仮説が正確に成り立たない場合、標本サイズが十分大きければ、帰無仮説からの小さな偏差でも棄却します。無同等性または最小効果検定は、標本が無限大に近づくと、比率が指定された無視できない量を超えて異ならない場合、帰無仮説を棄却しません。

この例では、点仮説と区間仮説の両方が棄却されません。比率が統計的に異なるとは言えません。

その後、標本サイズを20倍に増加させ、観測された比率は一定に保ちます。この場合、点帰無仮説は棄却され、 p 値は 0.01 となり、区間帰無仮説は棄却されず、 p 値は 1 です。

注: 無同等性検定は一般に保守的であり、そのサイズは α で制約されますが、標本サイズが大きくなると（標本が無限大に近づくと）、無同等性の限界が固定されている場合、サイズは α / 2 に近づきます。無同等性区間が 1 点に収縮すると、無同等性検定は点仮説検定と同じになります。（ドキュメント参照）

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

[24]:

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

[24]:

1.1036704302254083

[25]:

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

[25]:

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

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