Calculating adjusted p-values in Python

Question

So, I've been spending some time looking for a way to get adjusted p-values (aka corrected p-values, q-values, FDR) in Python, but I haven't really found anything. There's the R function p.adjust, but I would like to stick to Python coding, if possible. Is there anything similar for Python?

If this is somehow a bad question, sorry in advance! I did search for answers first, but found none (except a Matlab version)... Any help is appreciated!

Answer 1

It is available in statsmodels.

http://statsmodels.sourceforge.net/devel/stats.html#multiple-tests-and-multiple-comparison-procedures

http://statsmodels.sourceforge.net/devel/generated/statsmodels.sandbox.stats.multicomp.multipletests.html

and some explanations, examples and Monte Carlo http://jpktd.blogspot.com/2013/04/multiple-testing-p-value-corrections-in.html

Answer 2

According to the biostathandbook, the BH is easy to compute.

def fdr(p_vals):

    from scipy.stats import rankdata
    ranked_p_values = rankdata(p_vals)
    fdr = p_vals * len(p_vals) / ranked_p_values
    fdr[fdr > 1] = 1

    return fdr

Answer 3

You can try the module rpy2 that allows you to import R functions (b.t.w., a basic search returns How to implement R's p.adjust in Python).

Another possibility is to look at the maths an redo it yourself, because it is still relatively easy.

Apparently there is an ongoing implementation in statsmodels: http://statsmodels.sourceforge.net/ipdirective/_modules/scikits/statsmodels/sandbox/stats/multicomp.html . Maybe it is already usable.

Answer 4

You mentioned in your question q-values and no answer provided a link which addresses this. I believe this package (at least it seems so from the documentation) calculates q-values in python

https://puolival.github.io/multipy/

and also this one

https://github.com/nfusi/qvalue

Answer 5

We also added FDR in SciPy (will be in the next release, SciPy 1.11).

https://scipy.github.io/devdocs/reference/generated/scipy.stats.false_discovery_control.html

from scipy import stats

ps = [0.0001, 0.0004, 0.0019, 0.0095, 0.0201, 0.0278, 0.0298, 0.0344,
      0.0459, 0.3240, 0.4262, 0.5719, 0.6528, 0.7590, 1.000]
stats.false_discovery_control(ps)

# array([0.0015    , 0.003     , 0.0095    , 0.035625  , 0.0603    ,
#        0.06385714, 0.06385714, 0.0645    , 0.0765    , 0.486     ,
#        0.58118182, 0.714875  , 0.75323077, 0.81321429, 1.        ])