Skip to contents

Thin wrapper over stats::p.adjust(method = "BH").

Usage

benjamini_hochberg(p_values, alpha = 0.05, labels = NULL)

bh(p_values, alpha = 0.05, labels = NULL)

Arguments

p_values

Numeric vector of raw p-values.

alpha

Significance level.

labels

Optional character vector of test labels.

Value

An object of class "morie_multiple_testing_result".

Examples

set.seed(1)
# 100 tests, 20 true effects: FDR control keeps more power than FWER.
p <- c(runif(80), runif(20, 0, 0.005))

res <- benjamini_hochberg(p)
res$n_rejected
#> [1] 20
head(res$adjusted)          # BH-adjusted q-values
#> [1] 0.7222180 0.7799791 0.8813129 0.9509124 0.6428561 0.9509124
sum(res$rejected)           # discoveries at the default alpha = 0.05
#> [1] 20

# BH controls the false discovery rate, so it rejects far more than the
# FWER methods on the same data:
c(BH = benjamini_hochberg(p)$n_rejected,
  holm = holm(p)$n_rejected,
  bonferroni = bonferroni(p)$n_rejected)
#>         BH       holm bonferroni 
#>         20          1          1 

# `bh()` is a shorthand alias; `alpha` sets the FDR level.
identical(bh(p)$adjusted, benjamini_hochberg(p)$adjusted)
#> [1] TRUE
benjamini_hochberg(p, alpha = 0.10)$n_rejected
#> [1] 21