python - Poisson 的 E-test 在 Python 中的实现

Question

是否有泊松电子检验的 Python 实现？对于二项式 scipy 有 Fisher's Exact test作为 stats.fisher_exact 和 Gaussians scipy.stats 有 Welch's T-test作为 ttest_ind。我似乎找不到任何用于比较两个泊松分布的 Python 实现。

For context look here

For the algorithm look here

For R implementation look here

Answer 1

更新 这些函数现在包含在 statsmodels 中 https://www.statsmodels.org/dev/generated/statsmodels.stats.rates.test_poisson_2indep.html
带有用于 TOST 等效测试的配套版本 https://www.statsmodels.org/dev/generated/statsmodels.stats.rates.tost_poisson_2indep.html

这是一个开始。这实现了 Gu 等人的三个测试。 al 2008 基于渐近正态分布，现在也基于精确分布的两个条件检验。

如果计数不是太小(比如大于 10 或 20)，并且曝光时间不是很不均等，则得分测试的效果相当好。对于较小的计数，结果可能有点保守或宽松，其他方法会更好。 “sqrt”版本在他们的模拟中表现非常好，但在该版本有效时可能比分数测试功能稍差。

'''Test for ratio of Poisson intensities in two independent samples

Author: Josef Perktold
License: BSD-3

destination statsmodels

'''


from __future__ import division
import numpy as np
from scipy import stats


# copied from statsmodels.stats.weightstats
def _zstat_generic2(value, std_diff, alternative):
    '''generic (normal) z-test to save typing

    can be used as ztest based on summary statistics
    '''
    zstat = value / std_diff
    if alternative in ['two-sided', '2-sided', '2s']:
        pvalue = stats.norm.sf(np.abs(zstat))*2
    elif alternative in ['larger', 'l']:
        pvalue = stats.norm.sf(zstat)
    elif alternative in ['smaller', 's']:
        pvalue = stats.norm.cdf(zstat)
    else:
        raise ValueError('invalid alternative')
    return zstat, pvalue


def poisson_twosample(count1, exposure1, count2, exposure2, ratio_null=1,
                      method='score', alternative='2-sided'):
    '''test for ratio of two sample Poisson intensities

    If the two Poisson rates are g1 and g2, then the Null hypothesis is

    H0: g1 / g2 = ratio_null

    against one of the following alternatives

    H1_2-sided: g1 / g2 != ratio_null
    H1_larger: g1 / g2 > ratio_null
    H1_smaller: g1 / g2 < ratio_null

    Parameters
    ----------
    count1: int
        Number of events in first sample
    exposure1: float
        Total exposure (time * subjects) in first sample
    count2: int
        Number of events in first sample
    exposure2: float
        Total exposure (time * subjects) in first sample
    ratio: float
        ratio of the two Poisson rates under the Null hypothesis. Default is 1.
    method: string
        Method for the test statistic and the p-value. Defaults to `'score'`.
        Current Methods are based on Gu et. al 2008
        Implemented are 'wald', 'score' and 'sqrt' based asymptotic normal
        distribution, and the exact conditional test 'exact-cond', and its mid-point
        version 'cond-midp', see Notes
    alternative : string
        The alternative hypothesis, H1, has to be one of the following

           'two-sided': H1: ratio of rates is not equal to ratio_null (default)
           'larger' :   H1: ratio of rates is larger than ratio_null
           'smaller' :  H1: ratio of rates is smaller than ratio_null

    Returns
    -------
    stat, pvalue two-sided

    not yet
    #results : Results instance
    #    The resulting test statistics and p-values are available as attributes.


    Notes
    -----
    'wald': method W1A, wald test, variance based on separate estimates
    'score': method W2A, score test, variance based on estimate under Null
    'wald-log': W3A
    'score-log' W4A
    'sqrt': W5A, based on variance stabilizing square root transformation
    'exact-cond': exact conditional test based on binomial distribution
    'cond-midp': midpoint-pvalue of exact conditional test

    The latter two are only verified for one-sided example.

    References
    ----------
    Gu, Ng, Tang, Schucany 2008: Testing the Ratio of Two Poisson Rates,
    Biometrical Journal 50 (2008) 2, 2008

    '''

    # shortcut names
    y1, n1, y2, n2 = count1, exposure1, count2, exposure2
    d = n2 / n1
    r = ratio_null
    r_d = r / d

    if method in ['score']:
        stat = (y1 - y2 * r_d) / np.sqrt((y1 + y2) * r_d)
        dist = 'normal'
    elif method in ['wald']:
        stat = (y1 - y2 * r_d) / np.sqrt(y1 + y2 * r_d**2)
        dist = 'normal'
    elif method in ['sqrt']:
        stat = 2 * (np.sqrt(y1 + 3 / 8.) - np.sqrt((y2 + 3 / 8.) * r_d))
        stat /= np.sqrt(1 + r_d)
        dist = 'normal'
    elif method in ['exact-cond', 'cond-midp']:
        from statsmodels.stats import proportion
        bp = r_d / (1 + r_d)
        y_total = y1 + y2
        stat = None
        pvalue = proportion.binom_test(y1, y_total, prop=bp, alternative=alternative)
        if method in ['cond-midp']:
            # not inplace in case we still want binom pvalue
            pvalue = pvalue - 0.5 * stats.binom.pmf(y1, y_total, bp)

        dist = 'binomial'

    if dist == 'normal':
        return _zstat_generic2(stat, 1, alternative)
    else:
        return stat, pvalue


from numpy.testing import assert_allclose

# testing against two examples in Gu et al

print('\ntwo-sided')
# example 1
count1, n1, count2, n2 = 60, 51477.5, 30, 54308.7

s1, pv1 = poisson_twosample(count1, n1, count2, n2, method='wald')
pv1r = 0.000356
assert_allclose(pv1, pv1r*2, rtol=0, atol=5e-6)
print('wald', s1, pv1 / 2)   # one sided in the "right" direction

s2, pv2 = poisson_twosample(count1, n1, count2, n2, method='score')
pv2r = 0.000316
assert_allclose(pv2, pv2r*2, rtol=0, atol=5e-6)
print('score', s2, pv2 / 2)   # one sided in the "right" direction

s2, pv2 = poisson_twosample(count1, n1, count2, n2, method='sqrt')
pv2r = 0.000285
assert_allclose(pv2, pv2r*2, rtol=0, atol=5e-6)
print('sqrt', s2, pv2 / 2)   # one sided in the "right" direction

print('\ntwo-sided')
# example2
# I don't know why it's only 2.5 decimal agreement, rounding?
count1, n1, count2, n2 = 41, 28010, 15, 19017
s1, pv1 = poisson_twosample(count1, n1, count2, n2, method='wald', ratio_null=1.5)
pv1r = 0.2309
assert_allclose(pv1, pv1r*2, rtol=0, atol=5e-3)
print('wald', s1, pv1 / 2)   # one sided in the "right" direction

s2, pv2 = poisson_twosample(count1, n1, count2, n2, method='score', ratio_null=1.5)
pv2r = 0.2398
assert_allclose(pv2, pv2r*2, rtol=0, atol=5e-3)
print('score', s2, pv2 / 2)   # one sided in the "right" direction

s2, pv2 = poisson_twosample(count1, n1, count2, n2, method='sqrt', ratio_null=1.5)
pv2r = 0.2499
assert_allclose(pv2, pv2r*2, rtol=0, atol=5e-3)
print('sqrt', s2, pv2 / 2)   # one sided in the "right" direction

print('\none-sided')
# example 1 onesided
count1, n1, count2, n2 = 60, 51477.5, 30, 54308.7

s1, pv1 = poisson_twosample(count1, n1, count2, n2, method='wald', alternative='larger')
pv1r = 0.000356
assert_allclose(pv1, pv1r, rtol=0, atol=5e-6)
print('wald', s1, pv1)   # one sided in the "right" direction

s2, pv2 = poisson_twosample(count1, n1, count2, n2, method='score', alternative='larger')
pv2r = 0.000316
assert_allclose(pv2, pv2r, rtol=0, atol=5e-6)
print('score', s2, pv2)   # one sided in the "right" direction

s2, pv2 = poisson_twosample(count1, n1, count2, n2, method='sqrt', alternative='larger')
pv2r = 0.000285
assert_allclose(pv2, pv2r, rtol=0, atol=5e-6)
print('sqrt', s2, pv2)   # one sided in the "right" direction

# 'exact-cond', 'cond-midp'

s2, pv2 = poisson_twosample(count1, n1, count2, n2, method='exact-cond',
                            ratio_null=1, alternative='larger')
pv2r = 0.000428  # typo in Gu et al, switched pvalues between C and M
assert_allclose(pv2, pv2r, rtol=0, atol=5e-4)
print('exact-cond', s2, pv2)   # one sided in the "right" direction

s2, pv2 = poisson_twosample(count1, n1, count2, n2, method='cond-midp',
                            ratio_null=1, alternative='larger')
pv2r = 0.000310
assert_allclose(pv2, pv2r, rtol=0, atol=5e-4)
print('cond-midp', s2, pv2)   # one sided in the "right" direction


print('\none-sided')
# example2 onesided
# I don't know why it's only 2.5 decimal agreement, rounding?
count1, n1, count2, n2 = 41, 28010, 15, 19017
s1, pv1 = poisson_twosample(count1, n1, count2, n2, method='wald',
                            ratio_null=1.5, alternative='larger')
pv1r = 0.2309
assert_allclose(pv1, pv1r, rtol=0, atol=5e-4)
print('wald', s1, pv1)   # one sided in the "right" direction

s2, pv2 = poisson_twosample(count1, n1, count2, n2, method='score',
                            ratio_null=1.5, alternative='larger')
pv2r = 0.2398
assert_allclose(pv2, pv2r, rtol=0, atol=5e-4)
print('score', s2, pv2)   # one sided in the "right" direction

s2, pv2 = poisson_twosample(count1, n1, count2, n2, method='sqrt',
                            ratio_null=1.5, alternative='larger')
pv2r = 0.2499
assert_allclose(pv2, pv2r, rtol=0, atol=5e-4)
print('score', s2, pv2)   # one sided in the "right" direction

# 'exact-cond', 'cond-midp'

s2, pv2 = poisson_twosample(count1, n1, count2, n2, method='exact-cond',
                            ratio_null=1.5, alternative='larger')
pv2r = 0.2913
assert_allclose(pv2, pv2r, rtol=0, atol=5e-4)
print('exact-cond', s2, pv2)   # one sided in the "right" direction

s2, pv2 = poisson_twosample(count1, n1, count2, n2, method='cond-midp',
                            ratio_null=1.5, alternative='larger')
pv2r = 0.2450
assert_allclose(pv2, pv2r, rtol=0, atol=5e-4)
print('cond-midp', s2, pv2)   # one sided in the "right" direction

这打印

two-sided /2
wald 3.38491255626 0.000356004664253
score 3.417401839 0.000316109441024
sqrt 3.44548501956 0.00028501778109

two-sided /2
wald 0.73544663636 0.231033764105
score 0.706630933035 0.239897930348
sqrt 0.674401392575 0.250028078819

one-sided
wald 3.38491255626 0.000356004664253
score 3.417401839 0.000316109441024
sqrt 3.44548501956 0.00028501778109

one-sided
wald 0.73544663636 0.231033764105
score 0.706630933035 0.239897930348
score 0.674401392575 0.250028078819

精确条件测试相对容易实现，但非常保守且功耗低。大致精确的测试需要更多的努力(我目前没有时间)。

(通常:实际计算是几行。确定接口(interface)、添加文档和单元测试是更多工作。)

编辑

上面的脚本现在还包括精确的条件测试和它的中点 p 值版本，用 Gu 等人的单边替代的两个例子进行了检查。

例子一:单边

exact-cond None 0.00042805269405
cond-midp None 0.000310132441983

例子2:单边

exact-cond None 0.291453753765
cond-midp None 0.245173718501

目前返回的条件测试没有测试统计