单样本和两样本泊松率的统计与推断¶

作者: Josef Perktold

本笔记本简要概述了在单样本和双样本情况下泊松率的假设检验、置信区间和其他统计方法。更多选项和详细信息请参阅文档字符串。

在 statsmodels.stats.rates 中的所有函数都接受数据的汇总统计作为参数。这些参数包括事件的计数、观察次数或总暴露量。一些泊松分布的函数具有超分散选项。负二项分布（NB2）的函数需要分散参数。超分散和分散参数需要由用户提供，并且可以通过GLM-泊松和离散负二项分布模型分别从原始数据中估计。

注意，某些部分仍然是实验性的，可能会发生变化，某些功能仍然缺失，将在未来的版本中添加。

One sample functions

Two sample functions

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

单样本函数¶

The main functions for one sample Poisson rates currently are test_poisson and confint_poisson. Both have several methods available, most of them are consistent between hypothesis test and confidence interval. Two additional functions are available for tolerance intervals and for confidence intervals of quantiles.

See docstrings for details.

[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.23489608509894777
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 用于泊松分位数的置信区间。

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

A tolerance interval is specified by two probabilities, prob is the coverage of the Poisson interval, alpha is the confidence level for the confidence interval of the rate estimate.

Note, that probabilities cannot be exactly equal to the nominal probabilites because counts are discrete random variables. The properties of the intervals are specified in term of inequalities, coverage is at least prob, coverage of the confidence interval of the estimated rate is at least 1 - alpha. However, most methods will not guarantee that the coverage inequalities hold in small samples even if the distribution is correctly specified.

在以下示例中，如果总暴露量或观察次数为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。

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

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

字典 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.000712009328506111
   score        0.0006322188820470973
   score-log    0.000399251966184898
   wald-log     0.0008399438093390381
   exact-cond   0.0006751826586863219
   cond-midp    0.0005572624066190538
   sqrt         0.0005700355621795108
   etest-score  0.0005672617581627985
   etest-wald   0.000643144612489785
diff
   wald         0.0007120093285061094
   score        0.0006322188820470941
   waldccv      0.0007610462660136598
   etest-score  0.0005672617581627926
   etest-wald   0.0006431446124897783

以类似的方式，我们可以为所有当前可用方法计算率比和率差的置信区间。

[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 函数实现了等效性检验，其中备择假设指定两个比率在彼此的区间内。

非等价性测试实现了一种测试，其中备择假设指定两个比率至少相差一个给定的非零值。这也通常被称为最小效应测试。该测试使用两个单侧测试，类似于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.7754797240370243e-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.7754797240370243e-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.19557869693808483)
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.19557869693808483)
    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.19557869693808483))
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.19557869693808483))

函数 nonequivalence_poisson_2indep 用于检验两个率之间存在非微小差异的备择假设。

在下面的示例中，备择假设指定率比在区间 (0.95, 1/0.95) 之外。原假设是比率比在区间内。如果检验拒绝原假设，则提供证据表明率比差异超过区间限制指定的不重要数量。

关于在大样本中点假设检验和区间假设检验之间关系的注释。如果原假设不完全成立且样本量足够大，test_poisson_2indep的点原假设检验将拒绝任何小的偏离原假设的情况。如果比率差异不超过指定的可忽略量，非等价性或最小效应检验在大样本（样本趋近于无穷大）中不会拒绝原假设。

在示例中，点假设和区间假设都没有被拒绝。我们没有足够的证据表明这些比率在统计上是不同的。随后，我们将样本量增加20倍，同时保持观察到的比率不变。在这种情况下，点假设检验被拒绝，p值为0.01，而区间假设没有被拒绝，p值等于1。

注意：非等价性检验通常是保守的，其大小受alpha限制，但在大样本极限且固定非等价性边际的情况下，大小趋近于alpha / 2。如果非等价区间缩小到一个点，则非等价性检验与点假设检验相同。（参见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.011085166380602694)

功率¶

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

Last update: Oct 16, 2024