Existe-t-il une implémentation Python du E-Test pour Poissons? Pour binomials scipy a le Fisher's Exact test comme stats.fisher_exact et pour les gaussiens scipy.stats a Welch's T-test comme ttest_ind. Je n'arrive pas à trouver une implémentation Python pour la comparaison de deux Poissons.Implémentation d'E-test pour Poisson en Python
For R implementation look here
est un début ici. Ceci implémente trois tests de Gu et. al 2008 basé sur la distribution normale asymptotique, et maintenant aussi deux test conditionnel basé sur la distribution exacte. Le test de score fonctionne raisonnablement bien si les comptes ne sont pas trop petits (disons plus grands que 10 ou 20), et les temps d'exposition ne sont pas très inégaux. Pour les petits comptes, les résultats peuvent être soit un peu conservateurs ou libéraux, et d'autres méthodes seraient mieux. La version 'sqrt' s'est très bien comportée dans leurs simulations, mais pourrait avoir un peu moins de puissance que le test de score quand celui-ci fonctionne.
'''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
Imprime
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
Le test conditionnel exacte serait relativement facile à mettre en œuvre mais elle est très conservatrice et a une faible puissance. Les tests à peu près exacts nécessiteront un peu plus d'effort (pour lequel je n'aurai pas le temps pour le moment).
(Le plus souvent. Le calcul réel sont quelques lignes Décider sur l'interface, ajout de tests de documentation et de l'unité est plus de travail.)
EDIT
Le script ci-dessus inclut maintenant l'exacte conditionnelle test et sa version de valeur p à mi-chemin, vérifiée avec les deux exemples d'alternative unilatérale dans Gu et al.
Exemple 1: à sens unique
exact-cond None 0.00042805269405
cond-midp None 0.000310132441983
Exemple 2: à sens unique
exact-cond None 0.291453753765
cond-midp None 0.245173718501
actuellement il n'y a pas de statistiques de test pour les tests conditionnels retourné
Mes chiffres sont rarement supérieurs à 20 et le temps d'exposition est égal à la précision de la machine. Je pense que je vais aller avec la méthode sqrt pour l'instant. Je préférerais la méthode exacte parce que je veux être conservateur. Ce sera un bon point de départ dans tous les cas. – Keith
Je ne comprends pas ce que signifie avoir un temps d'exposition égal à la précision de la machine. Le taux de Poisson sera relatif à n'importe quelle unité d'exposition.Travailler à la précision de la machine pourrait faire exploser le bruit numérique et il pourrait être préférable de changer l'unité d'exposition. – user333700
Tests conditionnels exacts semblent faciles, je viens d'essayer et ajouterai conditionnel et conditionnel à mi-p bientôt. La puissance dans les petits échantillons (comptes) ne sera importante que pour les grandes différences dans les taux de Poisson. – user333700
Je ne l'ai jamais vu dans python. Statsmodels GLM-poisson ou Poisson fournit un test basé sur la distribution asymptotique. J'ai récemment regardé dans C, E tests et similaires pour les taux et les proportions. Le test E de Poisson dans la référence ressemble aux tests de score analogues. https://github.com/statsmodels/statsmodels/issues/2607 – user333700
Je veux faire un test unilatéral sur deux échantillons pour tester qu'un échantillon provient d'une population avec une moyenne plus grande. c'est-à-dire un test A/B standard basé sur les valeurs p. Je suis tout à fait sûr qu'un test E serait idéal mais pour l'instant je serais très content de trouver un test de Poisson qui soit analogue aux tests que je liste ci-dessus. Veuillez répondre avec un extrait de code. – Keith