一个和两个样本泊松率的统计和推断

作者: Josef Perktold

此笔记本简要概述了在一个和两个样本情况下泊松率的假设检验、置信区间和其他统计数据。有关更多选项和详细信息，请参阅文档字符串。

所有函数 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]:

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

目前还有两个用于单样本泊松率的额外函数，tolerance_int_poisson 用于容忍区间和 confint_quantile_poisson 用于泊松分位数的置信区间。

容忍区间类似于预测区间，它将新观测值的随机性和对估计的泊松率的不确定性结合在一起。如果速率已知，那么我们可以使用给定速率下逆累积分布函数来计算新观测值的泊松区间。容忍区间通过使用速率估计的置信区间来增加速率的不确定性。

容忍区间由两个概率指定，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))

最后一个函数返回泊松分位数的置信区间。分位数是累积分布函数的逆，在 scipy.stats 分布中名为 ppf。

以下示例显示了累积分布函数概率为 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))

双样本函数

用于两个样本的统计函数可以通过比率或差值来比较速率。默认情况下是比较速率比率。

可以通过 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 对计算比较 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 日