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)
[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!

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.

Powered by the version 2.3.6