如何在Python中实现R的p.adjust

Question

我有一组p值，现在想要计算调整后的p值，以便进行多重比较，主要是为了控制假发现率（FDR）。在R语言中，我可以使用：

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 )

我想知道如何在Python中实现这个代码。我在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

但是上面的代码出现了一个错误，提示AttributeError: 'float' object has no attribute 'r'。有没有人能帮我指出问题出在哪里？非常感谢大家的帮助！

Answer 1

这是我自己写的一个函数：

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 2

这个问题有点老了，但在Python的statsmodels库中，有很多可以用来进行比较修正的方法。我们可以在这里找到相关信息：

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

Answer 3

如果你想确保从R中得到你想要的结果，你可以指定使用R包中的'stats'这个函数：

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