converting very small p-values into chi-squared statistics with python
2
0
Entering edit mode
3.8 years ago

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


...even though:

> qchisq(1-3.760e-18, 1)
 Inf


Is there a way to make the python method behave like the R function?

Thanks!

python chi-squared p-value upper tail gwas • 2.3k views
ADD COMMENT
1
Entering edit mode

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

ADD REPLY
1
Entering edit mode
3.8 years ago

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

ADD COMMENT
0
Entering edit mode
3.8 years ago
Samuel Brady ▴ 320

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)

ADD COMMENT
0
Entering edit mode

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

ADD REPLY

Login before adding your answer.

Traffic: 1682 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.