# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# Adapted from [Thresholding with false discovery rate](https://matthew-brett.github.io/teaching/fdr.html) :
# https://matthew-brett.github.io/teaching/fdr.html

# Load libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats as sst


# +
# Generate samples
def p_samples(N,loc=0, scale=1., N_signal = 0, seed = 12345 ):
    
    np.random.seed(seed) # make it replicable
    
    # Data from a random normal distribution
    normal_distribution = sst.norm(loc=loc,scale=scale) #loc is the mean, scale is the variance.
    
    N_noise = N - N_signal
    # simulate low z scores
    noise_z_values = np.random.normal (size=N_noise)
    signal_z_values = np.random.normal (loc=-2.5, size=N_signal)
    mixed_z_values = np.sort (np.concatenate ((noise_z_values, signal_z_values)))
    mixed_p_values = normal_distribution.cdf (mixed_z_values)
    return mixed_p_values

p_values = p_samples(100, N_signal=20, seed=42)

# +
# initialize parameters
N = 100
alpha = 0.05
N_signal = 20
q = alpha

# Sort p-values
p_values = np.sort(p_values)
n = len(p_values)
i = np.arange(1, n+1)

# Rank to handle equal p-values 
helper_df = pd.DataFrame(p_values)
rank = round(helper_df.rank(axis=0, method = 'min')[0])

# -

# Test the significant values
max_threshold_pval = p_values < (q * rank / n)

# +
# Adjusted p-values that must be below the significance level
adj_pvals = p_values * n / rank

# Adjusted significance level
critical_val = (q * rank / n)
# -

# Dataframe output
df = pd.DataFrame()
df['index'] = i 
df['rank'] = rank
df['p_value'] = p_values
df['adjusted_pval'] = adj_pvals
df

# Verify the greatest threshold value's index
max_threshold_ind = p_values[p_values < critical_val].argmax(axis=0)
max_threshold_ind

# Showing the data just for fun, unnecessary
p_correction = zip(i, rank, np.round_(critical_val,3), np.round_(p_values,3), np.round_(adj_pvals,3))
for tup in (p_correction):
    print(tup)

# +
data = df
pv_index = 0
method = 'bh'
if isinstance(data, pd.DataFrame):
    data.rename({pv_index: "p_value"})
    ## error for non-numeric data frame column
    if not (np.issubdtype(data['p_value'].dtypes, np.number)):
        raise TypeError("Please ensure you have specified the column index of numeric p-values.")
else:
    data = pd.DataFrame({"p_value": data})



# +
def p_adjust(data, pv_index, method, alpha=0.05):
    """
    A summary dataframe with columns for the p-values, adjusted p-values for both Bonferroni and
    Benjamini-Hochberg (BH), adjusted significancelevel for Bonferroni and the critical value for BH

   Args:
        - data (dataframe): dataframe containing at least a column of p-values to be adjusted
        - pv_index (int): original p-value column index from existing input dataframe
        - alpha (int): significance level as a value between 0 and 1

   Returns:
        Dataframe: appends to input dataframe both adjusted p-values and significance levels (Bonferroni and BH)
        in ascending raw p-value order.Includes following columns:
            - bonf_val (int): Bonferroni adjusted significance level (same for all)
            - Bonf_significant (bool): True if significant p-value or False if not
            - bh_val (int): Benjamini-Hochberg (BH) critical value
            - BH_significant (bool): True if significant p-value or False if not
    """

#     ####if it's a pd.dataframe, rename to col header
#     if isinstance(data, pd.DataFrame):
#         data.rename({pv_index: "p_value"})
#         if np.issubdtype(data['p_value'].dtypes, np.number):
            
#     ###or make a vector a pd.dataframe
#     else:
#         data = pd.DataFrame({"p_value": data})
    if isinstance(data, pd.DataFrame):
        data.rename({pv_index: "p_value"})
        ## error for non-numeric data frame column
        if not (np.issubdtype(data['p_value'].dtypes, np.number)):
            raise TypeError("Please ensure you have specified the column index of numeric p-values.")
    else:
        data = pd.DataFrame({"p_value": data})
        # set the size of the data
    
    m = data.shape[0]

    # sort p-values
    df =data.sort_values(by=['p_value'])
    df["rank"]=round(df.rank(axis=0, method = 'min')["p_value"])
    df["bh_value"] = alpha*df["rank"]/m


    ### generate final data frame
    df["bonf_pvalue"] = np.where(df['p_value']*m<1,df['p_value']*m,1)
    df["bh_pvalue"] = df['p_value']*rank/m
    df = df.drop(['rank'], axis=1)
    
    
    if method == 'bh' or method == 'fdr':
        return(df[['p_value','bh_pvalue']])
    if method =='bonf' or method == 'bonferroni':
        return(df[['p_value','bonf_pvalue']])
    else:
        raise ValueError("Method should be set as 'bonf' or 'bh' corrections.")
    

# -

test = p_adjust(data, pv_index, method='error?', alpha=0.05)
# test[['bh_significant', 'bonf_significant']]

test
# ERROR METHOD: input invalid method


# +

def p_bonferroni_helper():
    """
    Applies Bonferroni correction to the original p-values

    Args:
    - pvals (int): original p-value column index from existing input dataframe
    - alpha (int): significance level as a value between 0 and 1

    Returns:
    - vector: returns the Bonferroni adjusted p-value
    """
    # vector output
    adj_bonf = alpha / pvals

    # dataframe output
    bonf_df = pd.DataFrame()
    bonf_df['p_value'] = p_values
    bonf_df['adjusted_pval'] = adj_bonf

    return bonf_df

def p_adjust():
    """
    This function executes an specific p-value adjustment method.

    Args:
        - data: a Pandas dataframe or a 1-d Numpy array.
        - col: if the input is a dataframe, col refers to the column name of the dataframe that has the p-values.
        - alpha: significance level for both methods.
        - method: method used for the adjustment ("bh" or "bonf")

    Returns:
        - data frame: data frame with the following columns:
            raw_p_value: original p-values.
            adjusted_p_value: p-values after the adjustment.
            signficant: boolean values determining if each p-value is significant.
            critical_value: it's the calculated critical value to compare with the cut-off.
    """
    # dataframe check
    if isinstance(data, pd.DataFrame):
        # numeric check
        if is_numeric_array(data):
            data.rename({pv_index: "p_value"})

    else:
        ## error for non-numeric data frame column
        raise ValueError("Please ensure you have specified the column index of numeric p-values.")

    else if (is.vector(data) | is.matrix(data)):
        if is.numeric(data):
            data = pd.DataFrame({"p_value": data})

    ###set the size of the data
    m = data.shape[0]

    return(m)

    pass




if self.data is None or type(self.data) is not np.ndarray:
            raise TypeError("Data is missing or of wrong object type. Please specify data as a Numpy array.")

        if self.K is None or self.K%1 != 0:
            raise TypeError("K is missing or is not an integer. Please specify K as an integer.")

        if self.data.shape[0] < self.K:
            raise ValueError("Cannot choose more initialize values than data observations.")
	   ####if it's a pd.dataframe, rename to col header
    if isinstance(data, pd.DataFrame):
        data.rename({pv_index: "p_value"})
    ###or make a vector a pd.dataframe
    else:
        data = pd.DataFrame({"p_value": data})



# +
def p_adjust(data, pv_index=0, method='bonf', alpha=0.05):
    """
    A summary dataframe with columns for the p-values, adjusted p-values for both Bonferroni and
    Benjamini-Hochberg (BH), adjusted significancelevel for Bonferroni and the critical value for BH

   Args:
        - data (dataframe): dataframe containing at least a column of p-values to be adjusted
        - pv_index (int): original p-value column index from existing input dataframe
        - alpha (int): significance level as a value between 0 and 1

   Returns:
        Dataframe: appends to input dataframe both adjusted p-values and significance levels (Bonferroni and BH)
        in ascending raw p-value order.Includes following columns:
            - bonf_val (int): Bonferroni adjusted significance level (same for all)
            - Bonf_significant (bool): True if significant p-value or False if not
            - bh_val (int): Benjamini-Hochberg (BH) critical value
            - BH_significant (bool): True if significant p-value or False if not
    """

    #     ####if it's a pd.dataframe, rename to col header
    #     if isinstance(data, pd.DataFrame):
    #         data.rename({pv_index: "p_value"})
    #         if np.issubdtype(data['p_value'].dtypes, np.number):

    #     ###or make a vector a pd.dataframe
    #     else:
    #         data = pd.DataFrame({"p_value": data})
    if isinstance(data, pd.DataFrame):
        data.rename({pv_index: "p_value"})
        ## error for non-numeric data frame column
        if not (np.issubdtype(data['p_value'].dtypes, np.number)):
            raise TypeError("Please ensure you have specified the column index of numeric p-values.")
    else:
        data = pd.DataFrame({"p_value": data})
        # set the size of the data

    m = data.shape[0]

    # sort p-values
    df = data.sort_values(by=['p_value'])
    df["rank"] = round(df.rank(axis=0, method='min')["p_value"])
    df["bh_value"] = alpha * df["rank"] / m

    ### generate final data frame
    df["bonf_pvalue"] = np.where(df['p_value'] * m < 1, df['p_value'] * m, 1)
    df["bh_pvalue"] = df['p_value'] * rank / m
    df = df.drop(['rank'], axis=1)

    if method == 'bh' or method == 'fdr':
        return (df[['p_value', 'bh_pvalue']])
    if method == 'bonf' or method == 'bonferroni':
        return (df[['p_value', 'bonf_pvalue']])
    else:
        raise ValueError("Method should be set as 'bonf' or 'bh' corrections")


# -

test = p_adjust(data, pv_index, method = 'bh', alpha=0.05)

test


# +
def p_adjust(data, pv_index=0, method='bonf', alpha=0.05):                                                                     
    """                                                                                                                        
    A summary dataframe with columns for the p-values, adjusted p-values for both Bonferroni and                               
    Benjamini-Hochberg (BH), adjusted significancelevel for Bonferroni and the critical value for BH                           
                                                                                                                               
   Args:                                                                                                                       
        - data (dataframe): dataframe containing at least a column of p-values to be adjusted                                  
        - pv_index (int): original p-value column index from existing input dataframe                                          
        - alpha (int): significance level as a value between 0 and 1                                                           
                                                                                                                               
   Returns:                                                                                                                    
        Dataframe: appends to input dataframe both adjusted p-values and significance levels (Bonferroni and BH)               
        in ascending raw p-value order.Includes following columns:                                                             
            - bonf_val (int): Bonferroni adjusted significance level (same for all)                                            
            - Bonf_significant (bool): True if significant p-value or False if not                                             
            - bh_val (int): Benjamini-Hochberg (BH) critical value                                                             
            - BH_significant (bool): True if significant p-value or False if not                                               
    """                                                                                                                        
                                                                                                                               
    #     ####if it's a pd.dataframe, rename to col header                                                                     
    #     if isinstance(data, pd.DataFrame):                                                                                   
    #         data.rename({pv_index: "p_value"})                                                                               
    #         if np.issubdtype(data['p_value'].dtypes, np.number):                                                             
                                                                                                                               
    #     ###or make a vector a pd.dataframe                                                                                   
    #     else:                                                                                                                
    #         data = pd.DataFrame({"p_value": data})                                                                           
    if isinstance(data, pd.DataFrame):                                                                                         
        data.rename({pv_index: "p_value"})                                                                                     
        ## error for non-numeric data frame column                                                                             
        if not (np.issubdtype(data['p_value'].dtypes, np.number)):                                                             
            raise TypeError("Please ensure you have specified the column index of numeric p-values.")                          
    else:                                                                                                                      
        data = pd.DataFrame({"p_value": data})                                                                                 
        # set the size of the data                                                                                             
                                                                                                                               
    m = data.shape[0]                                                                                                          
                                                                                                                               
    # sort p-values                                                                                                            
    df = data.sort_values(by=['p_value'])                                                                                      
    df["rank"] = round(df.rank(axis=0, method='min')["p_value"])                                                               
    df["bh_value"] = alpha * df["rank"] / m                                                                                    
                                                                                                                               
    ### generate final data frame                                                                                              
    df["bonf_pvalue"] = np.where(df['p_value'] * m < 1, df['p_value'] * m, 1)                                                  
    df["bh_pvalue"] = df['p_value'] * df['rank'] / m                                                                           
                                                                                                                               
                                                                                                                               
#     if method == 'bh' or method == 'fdr':                                                                                      
#         df["adjusted"] = np.where(df['p_value'] * m < 1, df['p_value'] * m, 1)                                                 
#         return (df[['p_value', 'adjusted']])                                                                                   
#     if method == 'bonf' or method == 'bonferroni':                                                                             
#         df["adjusted"] = df['p_value'] * df['rank'] / m                                                                        
#         return (df[['p_value', 'adjusted']])                                                                                   
#     else:                                                                                                                      
#         raise ValueError("Method should be set as 'bonf' or 'bh' corrections")                                                 
test1 = p_adjust(data, pv_index, method = 'bh', alpha=0.05)
test1
# -

d = { "p_value": [0.07],"adjusted": [0.07]}
df = pd.DataFrame(data=d)
assert df.equals(p_adjust(data=[0.07], method="bonf")) == False
test1


# +
# assert p_adjust(data=[0.07], method="bh") == df, "p_adjust single value vector for bh"
# -

p_adjust(data=[0.07], method="bonf")

df

# +
d = {"p_value": [0.07], "adjusted": [0.07]}

df = pd.DataFrame(data=d)

df = df[["p_value", "adjusted"]]

assert df.equals(p_adjust(data=[0.07], method="bonf"))
# -

df[["p_value", "adjusted"]]

p_adjust(data=[0.07], method="bonf")

d = {"p_value": [0.07, 0.2], "adjusted": [0.14, 0.4]}
df = pd.DataFrame(data=d)
df = df[["p_value", "adjusted"]]
assert df.equals(p_adjust(data=[0.07, 0.2], method="bonf"))

d = {"p_value": [0.07, 0.2], "adjusted": [0.14, 0.2]}
df = pd.DataFrame(data=d)
df = df[["p_value", "adjusted"]]
df

p_adjust(data=[0.07, 0.2], method="bh")

0.20 * 2 /2

StarcoderdataJupyterScriptsDedupFiltered