converting very small p-values into chi-squared statistics with python

Question

[yorgos.athanasiadis]

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!

[Jean-Karim Heriche]

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).

Answer 1

[yorgos.athanasiadis]

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

Answer 2

[Samuel Brady]

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)

[yorgos.athanasiadis]

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