Hi,
I am using scipy.stats.chi2.ppf in python to convert GWAS p-values into chi-squared statistics. As I am interested in the upper tail, I insert 1-p in the method.
The method works generally well. However, for very small p-values (e.g. 3.760e-18) and 1 df, the method returns Inf:
>>> scipy.stats.chi2.ppf(1-3.760e-18, 1)
inf
On the other hand, R's qchisq function takes care of this problem with the lower.tail=False option:
> qchisq(3.760e-18, 1, lower.tail=F)
[1] 75.44357
...even though:
> qchisq(1-3.760e-18, 1)
[1] Inf
Is there a way to make the python method behave like the R function?
Thanks!
After searching, I found the answer myself :-)
One can use the inverse survival function on the original p-value:
>>> scipy.stats.chi2.isf(3.760e-18, 1)
75.443573728536265
It looks like you've reached a precision limit in Python as indicated by Jean-Karim.
Try the following:
from decimal import *
decimalP = Decimal(1) - Decimal(3.760e-18)
scipy.stats.chi2.ppf(decimalP, 1)
Thank you Samuel, I had tried that already, but the ppf method doesn't support the Decimal type. I found another solution: use the isf method by scipy.stats
Log in to answer this question.
Do you mean something equivalent to lower.tail = FALSE for python ? I can't answer that question but I'd like to shed some light on the issue you're facing.
The problem is that 3.76e-18 is below machine precision (defined as the smallest x such that 1+x!=1) so 1-3.76e-18 = 1. However, due to the way numbers are represented, 3.76e-18 is not 0 so 0+3.76e-18=3.76e-18 not 0. So with lower.tail = FALSE, you're looking for x such that P(X>x) = 3.76e-18 whereas otherwise, you're looking for x such that P(X≤x) = 1 (and P(X≤Inf) = 1).