ポアソンの E-検定の Python 実装はありますか? Binomials の場合、scipy にはstats.fisher_exact としてFisher の Exact 検定があり、Gaussians の場合、scipy.stats には ttest_ind としてWelch の T 検定があります。2 つのポアソンを比較するための Python 実装が見つからないようです。
質問する
3760 次
2 に答える
10
更新これらの関数は現在、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らの3つのテストを実装しています。al 2008 は漸近正規分布に基づいており、現在は正確な分布に基づいた 2 つの条件付き検定も行われています。
スコア テストは、カウントが小さすぎず (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 et al. の片側代替の 2 つの例でチェックされています。
例 1: 片側
exact-cond None 0.00042805269405
cond-midp None 0.000310132441983
例 2: 片側
exact-cond None 0.291453753765
cond-midp None 0.245173718501
現在、返された条件付きテストのテスト統計はありません
于 2015-12-01T05:28:32.123 に答える
4
公開された論文に基づいて、クリシュナムーシーのfortran コードを numpy バインディングを使用してラップし、パッケージ化しました。ソースコードは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)
于 2019-03-24T23:00:06.887 に答える