Python 泊松电子检验的实施
Implementation of E-test for Poisson in Python
是否有 Python 泊松电子检验的实施?对于二项式 scipy 有 Fisher's Exact test as stats.fisher_exact and for Gaussians scipy.stats has Welch's T-test 作为 ttest_ind。我似乎找不到任何 Python 用于比较两个泊松的实现。
更新 这些函数现在包含在 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
精确条件测试相对容易实现,但非常保守且功率低。大致精确的测试需要更多的努力(我现在没有时间)。
(通常:实际计算是几行。确定接口,添加文档和单元测试是更多工作。)
编辑
上面的脚本现在还包括精确的条件测试和它的中点 p 值版本,用 Gu 等人的两个单边替代示例进行了检查
例一:单方面
exact-cond None 0.00042805269405
cond-midp None 0.000310132441983
例2:单方面
exact-cond None 0.291453753765
cond-midp None 0.245173718501
当前返回的条件测试没有测试统计
我包裹了 Krishnamoorthy 的 fortran code, based on the published paper, using numpy bindings and packaged it up. Source code is on github。
通过
安装
pip install poisson-etest
用法
from poisson_etest import poisson_etest
sample1_k, sample1_n = 10, 20
sample2_k, sample2_n = 15, 20
prob = poisson_etest(sample1_k, sample2_k, sample1_n, sample2_n)
是否有 Python 泊松电子检验的实施?对于二项式 scipy 有 Fisher's Exact test as stats.fisher_exact and for Gaussians scipy.stats has Welch's T-test 作为 ttest_ind。我似乎找不到任何 Python 用于比较两个泊松的实现。
更新 这些函数现在包含在 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
精确条件测试相对容易实现,但非常保守且功率低。大致精确的测试需要更多的努力(我现在没有时间)。
(通常:实际计算是几行。确定接口,添加文档和单元测试是更多工作。)
编辑
上面的脚本现在还包括精确的条件测试和它的中点 p 值版本,用 Gu 等人的两个单边替代示例进行了检查
例一:单方面
exact-cond None 0.00042805269405
cond-midp None 0.000310132441983
例2:单方面
exact-cond None 0.291453753765
cond-midp None 0.245173718501
当前返回的条件测试没有测试统计
我包裹了 Krishnamoorthy 的 fortran code, based on the published paper, using numpy bindings and packaged it up. Source code is on github。
通过
安装pip install poisson-etest
用法
from poisson_etest import poisson_etest
sample1_k, sample1_n = 10, 20
sample2_k, sample2_n = 15, 20
prob = poisson_etest(sample1_k, sample2_k, sample1_n, sample2_n)