單樣本與雙樣本 Poisson 率的統計與推論¶

作者：Josef Perktold

此筆記本簡要概述了單樣本和雙樣本情況下 Poisson 率的假設檢定、信賴區間和其他統計資訊。有關更多選項和其他詳細資訊，請參閱文件字串。

statsmodels.stats.rates 中的所有函數都將資料的摘要統計資訊作為引數。這些是事件計數和觀察次數或總曝露量。一些用於 Poisson 的函數具有過度離散的選項。用於負二項式（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,
    )

單樣本函數¶

目前單樣本 Poisson 率的主要函數是 test_poisson 和 confint_poisson。兩者都有多種方法可供使用，其中大多數方法在假設檢定和信賴區間之間保持一致。另有兩個函數可用於容差區間和分位數的信賴區間。

有關詳細資訊，請參閱文件字串。

[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 中。有關詳細資訊，請參閱文件字串。

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

目前另有兩個函數可用於單樣本 Poisson 率：用於容差區間的 tolerance_int_poisson 和用於 Poisson 分位數信賴區間的 confint_quantile_poisson。

容差區間類似於預測區間，它結合了新觀測值的隨機性和估計 Poisson 率的不確定性。如果已知速率，則可以使用給定速率下的反 cdf 計算新觀測值的 Poisson 區間。容差區間透過使用速率估計的信賴區間來增加速率的不確定性。

容差區間由兩個機率指定，prob 是 Poisson 區間的涵蓋範圍，alpha 是速率估計信賴區間的信賴水準。

請注意，機率不能完全等於名義機率，因為計數是離散隨機變數。區間的屬性以不等式的形式指定，涵蓋範圍至少為 prob，估計速率的信賴區間的涵蓋範圍至少為 1 - alpha。但是，即使正確指定了分佈，大多數方法也不能保證涵蓋不等式在小樣本中成立。

在以下範例中，如果總曝露量或觀察次數為 100，在給定的涵蓋範圍 prob 和信賴水準 alpha 下，我們可以預期觀察到 4 到 23 個事件。容差區間大於觀察速率下的 Poisson 區間 (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))

最後一個函數傳回 Poisson 分位數的信賴區間。分位數是 cdf 函數的反函數，在 scipy.stats 分佈中命名為 ppf。

以下範例顯示了 cdf 機率 0.975 下 Poisson 區間上限的信賴區間。使用單尾涵蓋機率的上限信賴限與容差區間的上限相同。

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

雙樣本函數¶

雙樣本的統計函數可以透過比率或差值比較速率。預設是比較速率比。

可以透過 test_poisson_2indep 直接存取 etest 函數。

[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，我們不能拒絕第一個速率是第二個速率的兩倍。

[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 字典顯示了在透過速率比或速率差比較兩個樣本時可用的方法。

我們可以使用字典來計算使用所有可用方法的 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))

我們有兩個額外的函數用於指定區間假設的假設檢定：tost_poisson_2indep 和 nonequivalence_poisson_2indep。

TOST 函數實作了等效性檢定，其中替代假設指定兩個速率彼此在一個區間內。

nonequivalence 檢定實作了一項檢定，其中替代假設指定兩個速率的差異至少為給定的非零值。這也常稱為最小效果檢定。此檢定使用兩個單邊檢定，類似於 TOST，但與等效性檢定相比，虛無假設和替代假設是相反的。

兩個函數都委派給 test_poisson_2indep，因此，可以使用相同的 method 選項。

以下等效性檢定指定替代假設，即速率比在 0.8 和 1/0.8 之間。觀察到的速率比為 0.89。 p 值為 0.107，我們無法拒絕虛無假設，轉而支持兩個速率在給定邊際下等效的替代假設。因此，假設檢定並未提供證據表明兩個速率是等效的。

在第二個範例中，我們檢定速率差中的等效性，其中等效性由邊際 (-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 = 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 檢定替代假設，即兩個速率的差異非小到可以忽略。

在以下範例中，替代假設指定速率比在區間 (0.95, 1/0.95) 之外。虛無假設是速率比在區間內。如果檢定拒絕虛無假設，則它會提供證據表明速率比的差異大於區間限制指定的不重要量。

關於大樣本中點假設檢定和區間假設檢定之間關係的注意事項。如果虛無假設不完全成立且樣本大小足夠大，則 test_poisson_2indep 的點虛無假設將拒絕與虛無假設的任何小偏差。如果速率的差異不超過指定的可忽略量，則在大型樣本中（樣本趨近無窮大），非等效性或最小效果檢定將不會拒絕虛無假設。

在範例中，點虛無假設和區間虛無假設均未被拒絕。我們沒有足夠的證據說明這些速率在統計上不同。接下來，我們將樣本大小增加 20 倍，同時保持觀察到的速率不變。在這種情況下，點虛無假設檢定被拒絕，p 值為 0.01，而區間虛無假設未被拒絕，p 值等於 1。

注意：非等效性檢定通常是保守的，其大小受 alpha 限制，但在具有固定非等效性邊際的大樣本限制中，大小接近 alpha / 2。如果非等效性區間縮小到一個點，則非等效性檢定與點假設檢定相同。（請參閱文件字串）

[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 對於計算雙樣本 Poisson 和負二項式速率比較的統計檢定力支援有限。這些檢定力以 Zhu 和 Lakkis 以及 Zhu 為基礎，用於兩種分佈的比率比較，以及基於基本常態的 Poisson 速率差比較。其他與假設檢定函數中的可用方法更密切對應的方法，尤其是 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 月 03 日