Cómo implementar R's p.adjust en Python

Question

Tengo una lista de valores p y me gustaría calcular los valores p de ajuste para comparaciones múltiples para FDR. En R, puedo usar:

pval <- read.csv("my_file.txt",header=F,sep="\t") 
pval <- pval[,1] 
FDR <- p.adjust(pval, method= "BH") 
print(length(pval[FDR<0.1])) 
write.table(cbind(pval, FDR),"pval_FDR.txt",row.names=F,sep="\t",quote=F) 

¿Cómo puedo implementar este código en Python? Aquí fue mi intento Feable en Python con la ayuda de Google:

pvalue_list [2.26717873145e-10, 1.36209234286e-11 , 0.684342083821...] # my pvalues 
pvalue_lst = [v.r['p.value'] for v in pvalue_list] 
p_adjust = R.r['p.adjust'](R.FloatVector(pvalue_lst),method='BH') 
for v in p_adjust: 
    print v 

El código anterior genera un error AttributeError: 'float' object has no attribute 'r'. ¿Alguien puede ayudarme a señalar mi problema? Gracias de antemano por la ayuda!

Answer 1

Si desea estar seguro de lo que está recibiendo de R, también puede indicar que desea utilizar la función en el paquete R 'alta':

from rpy2.robjects.packages import importr 
from rpy2.robjects.vectors import FloatVector 

stats = importr('stats') 

p_adjust = stats.p_adjust(FloatVector(pvalue_list), method = 'BH') 

Answer 2

Bueno, para obtener su código de trabajo, me imagino algo como esto funcionaría:

import rpy2.robjects as R 

pvalue_list = [2.26717873145e-10, 1.36209234286e-11 , 0.684342083821...] # my pvalues 
p_adjust = R['p.adjust'](R.FloatVector(pvalue_list),method='BH') 
for v in p_adjust: 
    print v 

Si p.adjust es bastante simple, se podría escribir en Python para que evitar la necesidad de poner en R. Y si quieres usarlo mucho, se puede hacer un simple envoltorio de Python:

def adjust_pvalues(pvalues, method='BH'): 
    return R['p.adjust'](R.FloatVector(pvalues), method=method) 

Answer 3

(sé que esto no es la respuesta ... tratando de ser útil.) el código de BH en I de p.adjust is just:

BH = { 
     i <- lp:1L # lp is the number of p-values 
     o <- order(p, decreasing = TRUE) # "o" will reverse sort the p-values 
     ro <- order(o) 
     pmin(1, cummin(n/i * p[o]))[ro] # n is also the number of p-values 
     } 

Answer 4

esta pregunta es un poco viejo , pero hay correcciones de comparación múltiples disponibles en los modelos de estadísticas para Python. Tenemos

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

Answer 5

vieja pregunta, pero aquí es una traducción del código R FDR en Python (que es probablemente bastante ineficiente):

def FDR(x): 
    """ 
    Assumes a list or numpy array x which contains p-values for multiple tests 
    Copied from p.adjust function from R 
    """ 
    o = [i[0] for i in sorted(enumerate(x), key=lambda v:v[1],reverse=True)] 
    ro = [i[0] for i in sorted(enumerate(o), key=lambda v:v[1])] 
    q = sum([1.0/i for i in xrange(1,len(x)+1)]) 
    l = [q*len(x)/i*x[j] for i,j in zip(reversed(xrange(1,len(x)+1)),o)] 
    l = [l[k] if l[k] < 1.0 else 1.0 for k in ro] 
    return l 

Answer 6

Aquí es una función de la casa que utilizo:

def correct_pvalues_for_multiple_testing(pvalues, correction_type = "Benjamini-Hochberg"):     
    """                         
    consistent with R - print correct_pvalues_for_multiple_testing([0.0, 0.01, 0.029, 0.03, 0.031, 0.05, 0.069, 0.07, 0.071, 0.09, 0.1]) 
    """ 
    from numpy import array, empty                   
    pvalues = array(pvalues) 
    n = float(pvalues.shape[0])                   
    new_pvalues = empty(n) 
    if correction_type == "Bonferroni":                 
     new_pvalues = n * pvalues 
    elif correction_type == "Bonferroni-Holm":                
     values = [ (pvalue, i) for i, pvalue in enumerate(pvalues) ]          
     values.sort() 
     for rank, vals in enumerate(values):                
      pvalue, i = vals 
      new_pvalues[i] = (n-rank) * pvalue                
    elif correction_type == "Benjamini-Hochberg":               
     values = [ (pvalue, i) for i, pvalue in enumerate(pvalues) ]          
     values.sort() 
     values.reverse()                     
     new_values = [] 
     for i, vals in enumerate(values):                 
      rank = n - i 
      pvalue, index = vals                   
      new_values.append((n/rank) * pvalue)               
     for i in xrange(0, int(n)-1): 
      if new_values[i] < new_values[i+1]:               
       new_values[i+1] = new_values[i]               
     for i, vals in enumerate(values): 
      pvalue, index = vals 
      new_pvalues[index] = new_values[i]                             
    return new_pvalues 

Answer 7

el uso de la biblioteca de Python numpy, sin necesidad de llamar a R en absoluto, he aquí una aplicación razonablemente eficiente del método BH:

import numpy as np 

def p_adjust_bh(p): 
    """Benjamini-Hochberg p-value correction for multiple hypothesis testing.""" 
    p = np.asfarray(p) 
    by_descend = p.argsort()[::-1] 
    by_orig = by_descend.argsort() 
    steps = float(len(p))/np.arange(len(p), 0, -1) 
    q = np.minimum(1, np.minimum.accumulate(steps * p[by_descend])) 
    return q[by_orig] 

(Basado en el código R BondedDust publicado)