pyprophet/stats.py
def lfdr(p_values, pi0, trunc = True, monotone = True, transf = "probit", adj = 1.5, eps = np.power(10.0,-8)):
    """ Estimate local FDR / posterior error probability from p-values according to bioconductor/qvalue """
    p = np.array(p_values)

    # Compare to bioconductor/qvalue reference implementation
    # import rpy2
    # import rpy2.robjects as robjects
    # from rpy2.robjects import pandas2ri
    # pandas2ri.activate()

    # density=robjects.r('density')
    # smoothspline=robjects.r('smooth.spline')
    # predict=robjects.r('predict')

    # Check inputs
    lfdr_out = p
    rm_na = np.isfinite(p)
    p = p[rm_na]

    if (min(p) < 0 or max(p) > 1):
        raise click.ClickException("p-values not in valid range [0,1].")
    elif (pi0 < 0 or pi0 > 1):
        raise click.ClickException("pi0 not in valid range [0,1].")

    # Local FDR method for both probit and logit transformations
    if (transf == "probit"):
        p = np.maximum(p, eps)
        p = np.minimum(p, 1-eps)
        x = scipy.stats.norm.ppf(p, loc=0, scale=1)

        # R-like implementation
        bw = bw_nrd0(x)
        myd = KDEUnivariate(x)
        myd.fit(bw=adj*bw, gridsize = 512)
        splinefit = sp.interpolate.splrep(myd.support, myd.density)
        y = sp.interpolate.splev(x, splinefit)
        # myd = density(x, adjust = 1.5) # R reference function
        # mys = smoothspline(x = myd.rx2('x'), y = myd.rx2('y')) # R reference function
        # y = predict(mys, x).rx2('y') # R reference function

        lfdr = pi0 * scipy.stats.norm.pdf(x) / y
    elif (transf == "logit"):
        x = np.log((p + eps) / (1 - p + eps))

        # R-like implementation
        bw = bw_nrd0(x)
        myd = KDEUnivariate(x)
        myd.fit(bw=adj*bw, gridsize = 512)

        splinefit = sp.interpolate.splrep(myd.support, myd.density)
        y = sp.interpolate.splev(x, splinefit)
        # myd = density(x, adjust = 1.5) # R reference function
        # mys = smoothspline(x = myd.rx2('x'), y = myd.rx2('y')) # R reference function
        # y = predict(mys, x).rx2('y') # R reference function

        dx = np.exp(x) / np.power((1 + np.exp(x)),2)
        lfdr = (pi0 * dx) / y
    else:
        raise click.ClickException("Invalid local FDR method.")

    if (trunc):
        lfdr[lfdr > 1] = 1
    if (monotone):
        lfdr = lfdr[p.ravel().argsort()]
        for i in range(1,len(x)):
            if (lfdr[i] < lfdr[i - 1]):
                lfdr[i] = lfdr[i - 1]
        lfdr = lfdr[scipy.stats.rankdata(p,"min")-1]

    lfdr_out[rm_na] = lfdr
    return lfdr_out
def lfdr(p_values, pi0, trunc = True, monotone = True, transf = "probit", adj = 1.5, eps = np.power(10.0,-8)):
    """ Estimate local FDR / posterior error probability from p-values according to bioconductor/qvalue """
    p = np.array(p_values)

    # Compare to bioconductor/qvalue reference implementation
    # import rpy2
    # import rpy2.robjects as robjects
    # from rpy2.robjects import pandas2ri
    # pandas2ri.activate()

    # density=robjects.r('density')
    # smoothspline=robjects.r('smooth.spline')
    # predict=robjects.r('predict')

    # Check inputs
    lfdr_out = p
    rm_na = np.isfinite(p)
    p = p[rm_na]

    if (min(p) < 0 or max(p) > 1):
        raise click.ClickException("p-values not in valid range [0,1].")
    elif (pi0 < 0 or pi0 > 1):
        raise click.ClickException("pi0 not in valid range [0,1].")

    # Local FDR method for both probit and logit transformations
    if (transf == "probit"):
        p = np.maximum(p, eps)
        p = np.minimum(p, 1-eps)
        x = scipy.stats.norm.ppf(p, loc=0, scale=1)

        # R-like implementation
        bw = bw_nrd0(x)
        myd = KDEUnivariate(x)
        myd.fit(bw=adj*bw, gridsize = 512)
        splinefit = sp.interpolate.splrep(myd.support, myd.density)
        y = sp.interpolate.splev(x, splinefit)
        # myd = density(x, adjust = 1.5) # R reference function
        # mys = smoothspline(x = myd.rx2('x'), y = myd.rx2('y')) # R reference function
        # y = predict(mys, x).rx2('y') # R reference function

        lfdr = pi0 * scipy.stats.norm.pdf(x) / y
    elif (transf == "logit"):
        x = np.log((p + eps) / (1 - p + eps))

        # R-like implementation
        bw = bw_nrd0(x)
        myd = KDEUnivariate(x)
        myd.fit(bw=adj*bw, gridsize = 512)

        splinefit = sp.interpolate.splrep(myd.support, myd.density)
        y = sp.interpolate.splev(x, splinefit)
        # myd = density(x, adjust = 1.5) # R reference function
        # mys = smoothspline(x = myd.rx2('x'), y = myd.rx2('y')) # R reference function
        # y = predict(mys, x).rx2('y') # R reference function

        dx = np.exp(x) / np.power((1 + np.exp(x)),2)
        lfdr = (pi0 * dx) / y
    else:
        raise click.ClickException("Invalid local FDR method.")

    if (trunc):
        lfdr[lfdr > 1] = 1
    if (monotone):
        lfdr = lfdr[p.ravel().argsort()]
        for i in range(1,len(x)):
            if (lfdr[i] < lfdr[i - 1]):
                lfdr[i] = lfdr[i - 1]
        lfdr = lfdr[scipy.stats.rankdata(p,"min")-1]

    lfdr_out[rm_na] = lfdr
    return lfdr_out
def lfdr ( p_values , pi0 , trunc = True , monotone = True , transf = "probit" , adj = 1.5 , eps = np . power ( 10.0 , - 8 ) ) : p = np . array ( p_values ) # Compare to bioconductor/qvalue reference implementation # import rpy2 # import rpy2.robjects as robjects # from rpy2.robjects import pandas2ri # pandas2ri.activate() # density=robjects.r('density') # smoothspline=robjects.r('smooth.spline') # predict=robjects.r('predict') # Check inputs lfdr_out = p rm_na = np . isfinite ( p ) p = p [ rm_na ] if ( min ( p ) < 0 or max ( p ) > 1 ) : raise click . ClickException ( "p-values not in valid range [0,1]." ) elif ( pi0 < 0 or pi0 > 1 ) : raise click . ClickException ( "pi0 not in valid range [0,1]." ) # Local FDR method for both probit and logit transformations if ( transf == "probit" ) : p = np . maximum ( p , eps ) p = np . minimum ( p , 1 - eps ) x = scipy . stats . norm . ppf ( p , loc = 0 , scale = 1 ) # R-like implementation bw = bw_nrd0 ( x ) myd = KDEUnivariate ( x ) myd . fit ( bw = adj * bw , gridsize = 512 ) splinefit = sp . interpolate . splrep ( myd . support , myd . density ) y = sp . interpolate . splev ( x , splinefit ) # myd = density(x, adjust = 1.5) # R reference function # mys = smoothspline(x = myd.rx2('x'), y = myd.rx2('y')) # R reference function # y = predict(mys, x).rx2('y') # R reference function lfdr = pi0 * scipy . stats . norm . pdf ( x ) / y elif ( transf == "logit" ) : x = np . log ( ( p + eps ) / ( 1 - p + eps ) ) # R-like implementation bw = bw_nrd0 ( x ) myd = KDEUnivariate ( x ) myd . fit ( bw = adj * bw , gridsize = 512 ) splinefit = sp . interpolate . splrep ( myd . support , myd . density ) y = sp . interpolate . splev ( x , splinefit ) # myd = density(x, adjust = 1.5) # R reference function # mys = smoothspline(x = myd.rx2('x'), y = myd.rx2('y')) # R reference function # y = predict(mys, x).rx2('y') # R reference function dx = np . exp ( x ) / np . power ( ( 1 + np . exp ( x ) ) , 2 ) lfdr = ( pi0 * dx ) / y else : raise click . ClickException ( "Invalid local FDR method." ) if ( trunc ) : lfdr [ lfdr > 1 ] = 1 if ( monotone ) : lfdr = lfdr [ p . ravel ( ) . argsort ( ) ] for i in range ( 1 , len ( x ) ) : if ( lfdr [ i ] < lfdr [ i - 1 ] ) : lfdr [ i ] = lfdr [ i - 1 ] lfdr = lfdr [ scipy . stats . rankdata ( p , "min" ) - 1 ] lfdr_out [ rm_na ] = lfdr return lfdr_out
Estimate local FDR / posterior error probability from p-values according to bioconductor/qvalue
Estimate local FDR / posterior error probability from p - values according to bioconductor / qvalue