pyprophet/stats.py
def pi0est(p_values, lambda_ = np.arange(0.05,1.0,0.05), pi0_method = "smoother", smooth_df = 3, smooth_log_pi0 = False):
    """ Estimate pi0 according to bioconductor/qvalue """

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

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

    p = np.array(p_values)

    rm_na = np.isfinite(p)
    p = p[rm_na]
    m = len(p)
    ll = 1
    if isinstance(lambda_, np.ndarray ):
        ll = len(lambda_)
        lambda_ = np.sort(lambda_)

    if (min(p) < 0 or max(p) > 1):
        raise click.ClickException("p-values not in valid range [0,1].")
    elif (ll > 1 and ll < 4):
        raise click.ClickException("If lambda_ is not predefined (one value), at least four data points are required.")
    elif (np.min(lambda_) < 0 or np.max(lambda_) >= 1):
        raise click.ClickException("Lambda must be within [0,1)")

    if (ll == 1):
        pi0 = np.mean(p >= lambda_)/(1 - lambda_)
        pi0_lambda = pi0
        pi0 = np.minimum(pi0, 1)
        pi0Smooth = False
    else:
        pi0 = []
        for l in lambda_:
            pi0.append(np.mean(p >= l)/(1 - l))
        pi0_lambda = pi0

        if (pi0_method == "smoother"):
            if smooth_log_pi0:
                pi0 = np.log(pi0)
                spi0 = sp.interpolate.UnivariateSpline(lambda_, pi0, k=smooth_df)
                pi0Smooth = np.exp(spi0(lambda_))
                # spi0 = smoothspline(lambda_, pi0, df = smooth_df) # R reference function
                # pi0Smooth = np.exp(predict(spi0, x = lambda_).rx2('y')) # R reference function
            else:
                spi0 = sp.interpolate.UnivariateSpline(lambda_, pi0, k=smooth_df)
                pi0Smooth = spi0(lambda_)
                # spi0 = smoothspline(lambda_, pi0, df = smooth_df) # R reference function
                # pi0Smooth = predict(spi0, x = lambda_).rx2('y')  # R reference function
            pi0 = np.minimum(pi0Smooth[ll-1],1)
        elif (pi0_method == "bootstrap"):
            minpi0 = np.percentile(pi0,0.1)
            W = []
            for l in lambda_:
                W.append(np.sum(p >= l))
            mse = (np.array(W) / (np.power(m,2) * np.power((1 - lambda_),2))) * (1 - np.array(W) / m) + np.power((pi0 - minpi0),2)
            pi0 = np.minimum(pi0[np.argmin(mse)],1)
            pi0Smooth = False
        else:
            raise click.ClickException("pi0_method must be one of 'smoother' or 'bootstrap'.")
    if (pi0<=0):
        raise click.ClickException("The estimated pi0 <= 0. Check that you have valid p-values or use a different range of lambda.")

    return {'pi0': pi0, 'pi0_lambda': pi0_lambda, 'lambda_': lambda_, 'pi0_smooth': pi0Smooth}
def pi0est(p_values, lambda_ = np.arange(0.05,1.0,0.05), pi0_method = "smoother", smooth_df = 3, smooth_log_pi0 = False):
    """ Estimate pi0 according to bioconductor/qvalue """

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

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

    p = np.array(p_values)

    rm_na = np.isfinite(p)
    p = p[rm_na]
    m = len(p)
    ll = 1
    if isinstance(lambda_, np.ndarray ):
        ll = len(lambda_)
        lambda_ = np.sort(lambda_)

    if (min(p) < 0 or max(p) > 1):
        raise click.ClickException("p-values not in valid range [0,1].")
    elif (ll > 1 and ll < 4):
        raise click.ClickException("If lambda_ is not predefined (one value), at least four data points are required.")
    elif (np.min(lambda_) < 0 or np.max(lambda_) >= 1):
        raise click.ClickException("Lambda must be within [0,1)")

    if (ll == 1):
        pi0 = np.mean(p >= lambda_)/(1 - lambda_)
        pi0_lambda = pi0
        pi0 = np.minimum(pi0, 1)
        pi0Smooth = False
    else:
        pi0 = []
        for l in lambda_:
            pi0.append(np.mean(p >= l)/(1 - l))
        pi0_lambda = pi0

        if (pi0_method == "smoother"):
            if smooth_log_pi0:
                pi0 = np.log(pi0)
                spi0 = sp.interpolate.UnivariateSpline(lambda_, pi0, k=smooth_df)
                pi0Smooth = np.exp(spi0(lambda_))
                # spi0 = smoothspline(lambda_, pi0, df = smooth_df) # R reference function
                # pi0Smooth = np.exp(predict(spi0, x = lambda_).rx2('y')) # R reference function
            else:
                spi0 = sp.interpolate.UnivariateSpline(lambda_, pi0, k=smooth_df)
                pi0Smooth = spi0(lambda_)
                # spi0 = smoothspline(lambda_, pi0, df = smooth_df) # R reference function
                # pi0Smooth = predict(spi0, x = lambda_).rx2('y')  # R reference function
            pi0 = np.minimum(pi0Smooth[ll-1],1)
        elif (pi0_method == "bootstrap"):
            minpi0 = np.percentile(pi0,0.1)
            W = []
            for l in lambda_:
                W.append(np.sum(p >= l))
            mse = (np.array(W) / (np.power(m,2) * np.power((1 - lambda_),2))) * (1 - np.array(W) / m) + np.power((pi0 - minpi0),2)
            pi0 = np.minimum(pi0[np.argmin(mse)],1)
            pi0Smooth = False
        else:
            raise click.ClickException("pi0_method must be one of 'smoother' or 'bootstrap'.")
    if (pi0<=0):
        raise click.ClickException("The estimated pi0 <= 0. Check that you have valid p-values or use a different range of lambda.")

    return {'pi0': pi0, 'pi0_lambda': pi0_lambda, 'lambda_': lambda_, 'pi0_smooth': pi0Smooth}
def pi0est ( p_values , lambda_ = np . arange ( 0.05 , 1.0 , 0.05 ) , pi0_method = "smoother" , smooth_df = 3 , smooth_log_pi0 = False ) : # Compare to bioconductor/qvalue reference implementation # import rpy2 # import rpy2.robjects as robjects # from rpy2.robjects import pandas2ri # pandas2ri.activate() # smoothspline=robjects.r('smooth.spline') # predict=robjects.r('predict') p = np . array ( p_values ) rm_na = np . isfinite ( p ) p = p [ rm_na ] m = len ( p ) ll = 1 if isinstance ( lambda_ , np . ndarray ) : ll = len ( lambda_ ) lambda_ = np . sort ( lambda_ ) if ( min ( p ) < 0 or max ( p ) > 1 ) : raise click . ClickException ( "p-values not in valid range [0,1]." ) elif ( ll > 1 and ll < 4 ) : raise click . ClickException ( "If lambda_ is not predefined (one value), at least four data points are required." ) elif ( np . min ( lambda_ ) < 0 or np . max ( lambda_ ) >= 1 ) : raise click . ClickException ( "Lambda must be within [0,1)" ) if ( ll == 1 ) : pi0 = np . mean ( p >= lambda_ ) / ( 1 - lambda_ ) pi0_lambda = pi0 pi0 = np . minimum ( pi0 , 1 ) pi0Smooth = False else : pi0 = [ ] for l in lambda_ : pi0 . append ( np . mean ( p >= l ) / ( 1 - l ) ) pi0_lambda = pi0 if ( pi0_method == "smoother" ) : if smooth_log_pi0 : pi0 = np . log ( pi0 ) spi0 = sp . interpolate . UnivariateSpline ( lambda_ , pi0 , k = smooth_df ) pi0Smooth = np . exp ( spi0 ( lambda_ ) ) # spi0 = smoothspline(lambda_, pi0, df = smooth_df) # R reference function # pi0Smooth = np.exp(predict(spi0, x = lambda_).rx2('y')) # R reference function else : spi0 = sp . interpolate . UnivariateSpline ( lambda_ , pi0 , k = smooth_df ) pi0Smooth = spi0 ( lambda_ ) # spi0 = smoothspline(lambda_, pi0, df = smooth_df) # R reference function # pi0Smooth = predict(spi0, x = lambda_).rx2('y')  # R reference function pi0 = np . minimum ( pi0Smooth [ ll - 1 ] , 1 ) elif ( pi0_method == "bootstrap" ) : minpi0 = np . percentile ( pi0 , 0.1 ) W = [ ] for l in lambda_ : W . append ( np . sum ( p >= l ) ) mse = ( np . array ( W ) / ( np . power ( m , 2 ) * np . power ( ( 1 - lambda_ ) , 2 ) ) ) * ( 1 - np . array ( W ) / m ) + np . power ( ( pi0 - minpi0 ) , 2 ) pi0 = np . minimum ( pi0 [ np . argmin ( mse ) ] , 1 ) pi0Smooth = False else : raise click . ClickException ( "pi0_method must be one of 'smoother' or 'bootstrap'." ) if ( pi0 <= 0 ) : raise click . ClickException ( "The estimated pi0 <= 0. Check that you have valid p-values or use a different range of lambda." ) return { 'pi0' : pi0 , 'pi0_lambda' : pi0_lambda , 'lambda_' : lambda_ , 'pi0_smooth' : pi0Smooth }
Estimate pi0 according to bioconductor/qvalue
Estimate pi0 according to bioconductor / qvalue