<jupyter_start><jupyter_text>Benjamini-Hochberg False Discovery Rate---  Publication: <NAME>., & <NAME>. (1995). Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. Journal of the Royal Statistical Society. Series B (Methodological), 57(1), 289-300. Retrieved from http://www.jstor.org/stable/2346101  In this notebook, I will demonstrate a geometric representation of BH FDR correction.  As describe in Wikipedia (https://en.wikipedia.org/wiki/False_discovery_rateBenjamini%E2%80%93Hochberg_procedure): --- Formula for BH FDR For a given $\alpha$, find the largest $k$ such that $P_{(k)}\leq {\frac {k}{m}}\alpha$ where $P_{(k)}$ is our p-value, $k$ is the rank of that p-value, $m$ is the number of p-values, and $\alpha$ is our significance cutoff. Reject the null hypothesis (i.e., declare discoveries) for all $H_{(i)}$ for $i = 1 , \ldots, k$. --- What we'll see: Geometrically, this FDR procedure corresponds to plotting $P_{(k)}$ vs. $k$ (on the y and x axes respectively), drawing the line through the origin with slope ${\frac {\alpha }{m}}$, and rejecting all null hypotheses that are to the left of where the points cross the line.<jupyter_code># First, we'll import all of the necessary libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.sandbox.stats.multicomp import multipletests
# Next, we'll open some example p-values -- these are unsorted and we won't worry about what they correspond to
# We will read them, sort them, and reset the index
df = pd.read_csv("p-values.txt").sort_values(by="p_values")
df.reset_index(inplace=True, drop=True)
print(len(df.p_values))
df.head()
# Here is a function I designed to perform the correction
def fdr_bh(df, alpha=0.05):
    """Given a dataframe with p_values, return dataframe with additional column of True/False if significant/not"""
    df["significant"]=[True if i <= ((df.index[x]+1)/len(df.index))*alpha else False for x,i in enumerate(df.p_values)]
    return df
df = fdr_bh(df)
df.head()
# y = slope*x + intercept
# slope = alpha/m
# intercept = 0
# Formula now:
# y = alpha/m*x + b

# x values go from 0-3361 (number of p-values)
x = np.linspace(0,3361,3362)
y = [(0.05/3362)*i for i in x]

plt.figure(figsize=(15,10))
plt.scatter(df.index, df.p_values, marker='o',c=df.significant, cmap=plt.cm.coolwarm, edgecolors='none')
plt.plot(x,y, color='g', linewidth=1)
plt.xlim(0,3390)
plt.ylim(0,1.1)
plt.ylabel("p-value")
plt.xlabel("rank")
plt.show()
# So, the points cross the line between dataframe index (or rank or x value) 1194 and 1195
# See below
df[1193:1197]

# As mentioned, everything to the left of this crossing -- all points in red -- 
# are signficantly different from the null
# Next, we'll re-open the data and use an existing package to perform FDR correction
df = pd.read_csv("p-values.txt").sort_values(by="p_values")
df.reset_index(inplace=True, drop=True)
print(len(df.p_values))
df.head()
# We'll correct our p-values using multipletests from statsmodels
# Returns array of True/False where True for hypothesis that can be rejected for given alpha
#         array of adjusted p-values
#         corrected alpha for Sidak method -- we'll ignore
#         corrected alpha for Bonferroni method -- we'll ignore
adj_p_values = multipletests(df.p_values, alpha=0.05, method="fdr_bh")
df["adj_p_values"]=adj_p_values[1]
df["significant"]=adj_p_values[0]
df.head()
# y = slope*x + intercept
# slope = alpha/m
# intercept = 0
# Formula now:
# y = alpha/m*x + b

# x values go from 0-3361 (number of p-values)
x = np.linspace(0,3361,3362)
y = [(0.05/3362)*i for i in x]

plt.figure(figsize=(15,10))
plt.scatter(df.index, df.p_values, marker='o',c=df.significant, cmap=plt.cm.coolwarm, edgecolors='none')
plt.plot(x,y, color='g', linewidth=1)
plt.xlim(0,3390)
plt.ylim(0,1.1)
plt.ylabel("p-value")
plt.xlabel("rank")
plt.show()
# Just as we saw above, the points cross the line between dataframe index (or rank or x value) 1194 and 1195
# See below
df[1193:1197]

# Everything to the left of this crossing -- all points in red -- 
# are signficantly different from the null<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12884753042233357, 'fraction_numerical': 0.032689095681221664, 'mean_word_length': 4.657219973009447, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 7, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '235561', 'license': 'MIT', 'n_tokens_mistral': 1525, 'n_tokens_neox': 1383, 'n_words': 493, 'path': 'Benjamini-Hochberg_FDR.ipynb', 'repo_name': 'jlboat/ipython-notebooks', 'size': '4202'}