benjamini_hochberg

benjamini_hochberg.R

benjamini_hochberg <- function(df, fdr){
  #input:
    #df = dataframe where the rows are the tests and one of the column named "pval"" contains the p-values
    #fdr = false discovery rate threshold
  # output
    # "threshold_BH" is the benjamini_hochberg threshold
    # "reject" indicates whether we reject the null or not in the output
  n = nrow(df)
  bh = df %>% mutate(rank = dense_rank(pval))  %>% mutate(threshold_BH = (rank/n)*fdr, reject = pval < threshold_BH) %>% arrange(pval)
  return(bh)
}

example = data.frame(pval = c(0.05,0.03,0.09), test_name = c("test1", "test2", "test3"))
View(benjamini_hochberg(example,0.05))