Estimates the proportion of true null hypotheses (pi0) and
tightens the BH thresholds by that factor. Native implementation of
Storey and Tibshirani (2003) that reproduces qvalue::qvalue
with its defaults: pi0 from a cubic smoothing spline (3 df, via
stats::smooth.spline) through \(\hat\pi_0(\lambda)\) over the
grid seq(0.05, 0.95, 0.05), read at the largest lambda, and
q-values \(q_{(i)} = \hat\pi_0 \min_{j \ge i} \min(1, m p_{(j)} / j)\).
Usage
storey_q(
p_values,
alpha = 0.05,
lambda_param = seq(0.05, 0.95, 0.05),
labels = NULL
)Arguments
- p_values
Numeric vector of raw p-values.
- alpha
Significance level.
- lambda_param
Tuning value(s) in
[0, 1)for the pi0 estimator, as qvalue'slambda: a single value gives the cutoff estimator \(N(\lambda)/(m(1-\lambda))\) (\(N(\lambda)\) = number of p-values \(\ge \lambda\)); a grid of at least four values (the default) gives the smoother estimator.- labels
Optional character vector of test labels.
Value
An object of class "morie_multiple_testing_result"
with the q-values in adjusted, plus pi0 and
lambda_param.
References
Storey, J. D. and Tibshirani, R. (2003). Statistical significance for genomewide studies. PNAS 100, 9440-9445.
Examples
set.seed(1)
p <- c(runif(80), runif(20, 0, 0.005))
res <- storey_q(p)
res$pi0 # estimated proportion of true nulls
#> [1] 0.6053452
head(res$adjusted) # q-values
#> [1] 0.4371912 0.4721566 0.5334985 0.5756302 0.3891499 0.5756302
res$n_rejected
#> [1] 21
# A single `lambda_param` gives the cutoff estimator at that lambda.
storey_q(p, lambda_param = 0.5)$pi0
#> [1] 0.78
storey_q(p, lambda_param = 0.8)$pi0
#> [1] 0.75
# Adaptive FDR rejects at least as many as BH (pi0 <= 1 tightens BH):
c(storey = storey_q(p)$n_rejected, BH = benjamini_hochberg(p)$n_rejected)
#> storey BH
#> 21 20