20

p 値のリストがあり、 FDRの複数の比較のために調整 p 値を計算したいと思います。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 の助けを借りて Python で実行した私の試みです。

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'エラーをスローします。誰でも私の問題を指摘できますか? 助けてくれてありがとう！

4

7 に答える 7

16

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')
于 2011-09-17T07:46:37.807 に答える
16

この質問は少し古いですが、Python の statsmodels で利用できる複数の比較修正があります。我々は持っています

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

于 2012-12-06T13:58:08.707 に答える
11

これが私が使用する社内関数です。

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
于 2014-02-12T21:05:53.247 に答える
9

R をまったく呼び出さずに Python の numpy ライブラリを使用すると、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]

(BondedDust掲載のRコードより)

于 2015-11-04T21:37:35.403 に答える
2

（これが答えではないことはわかっています...ただ役に立ちたいだけです。）Rのp.adjustのBHコードは次のとおりです。

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
      }
于 2011-09-16T22:53:44.407 に答える
1

古い質問ですが、Python での R FDR コードの翻訳は次のとおりです (これはおそらくかなり非効率的です)。

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
于 2014-01-27T14:55:39.720 に答える
0

さて、あなたのコードを機能させるために、私はこのようなものが機能すると思います：

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

p.adjustが十分に単純な場合は、Pythonで記述できるため、Rを呼び出す必要はありません。また、p.adjustを頻繁に使用する場合は、単純なPythonラッパーを作成できます。

def adjust_pvalues(pvalues, method='BH'):
    return R['p.adjust'](R.FloatVector(pvalues), method=method)
于 2011-09-16T22:47:49.393 に答える