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!
It is available in statsmodels.
and some explanations, examples and Monte Carlo http://jpktd.blogspot.com/2013/04/multiple-testing-p-value-corrections-in.html
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
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.
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. ])
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
