qqplotPval | R Documentation |
Produce a quantile-quantile plot for p values and display its confidence interval.
A quantile is an order statistic, and the j-th order statistic from a Uniform(0,1) sample has a Beta(j,N-j+1) distribution (Casella & Berger, 2001, 2nd edition, p230).
Let us assume we have N independent p values, \{p_1,\ldots,p_N\}
.
Under the null, they are independent and identically uniformly distributed: \forall i \; p_i \sim \mathcal{U}_{[0,1]}
.
Therefore, the 95% confidence interval for the j-th quantile of the set of p values can be calculated with: qbeta(0.95, j, N-j+1).
See also the qqman package.
qqplotPval(
pvalues,
use.density = FALSE,
nrpoints = 1000,
pch = 1,
plot.conf.int = TRUE,
xlab = expression(Expected ~ ~-log[10](italic(p) ~ values)),
ylab = expression(Observed ~ ~-log[10](italic(p) ~ values)),
thresh = 0.05,
ctl.fwer.bonf = FALSE,
ctl.fwer.sidak = FALSE,
ctl.fdr.bh = FALSE,
ctl.fdr.storey = FALSE,
pfdr = TRUE,
lim.custom,
lty.custom,
lgd.custom,
plot.signif = FALSE,
main = NULL,
col = NULL,
verbose = 1
)
pvalues |
vector of raw p values (missing values will be omitted) |
use.density |
if TRUE, uses |
nrpoints |
if |
pch |
vector of point symbol(s) (default is 1 for all points) |
plot.conf.int |
show the confidence interval |
xlab |
a title for the x axis (see default) |
ylab |
a title for the x axis (see default) |
thresh |
significance threshold at which to control the FWER and the FDR |
ctl.fwer.bonf |
control the family-wise error rate with the Bonferroni procedure |
ctl.fwer.sidak |
control the family-wise error rate with the Sidak procedure |
ctl.fdr.bh |
control the false discovery rate with the Benjamini-Hochberg procedure |
ctl.fdr.storey |
control the false discovery rate with Storey's procedure in the qvalue package |
pfdr |
logical passed to |
lim.custom |
custom limit on p values (optional) |
lty.custom |
line style for the custom limit (compulsory if |
lgd.custom |
legend for the custom limit (compulsory if |
plot.signif |
show line(s) corresponding to the significance threshold (compulsory if |
main |
an overall title for the plot (default: "Q-Q plot (<length(pvalues)> p values)") |
col |
vector of plotting color(s) for the points (default is all points in black) |
verbose |
verbosity level (0/1) |
invisible data.frame with a column of p values (NA omitted) and the adjusted p values if any of ctl.fwer.bonf
and ctl.fdr.bh
is set; look also at the attributes
Timothee Flutre (inspired by an anonymous comment at http://gettinggeneticsdone.blogspot.fr/2009/11/qq-plots-of-p-values-in-r-using-ggplot2.html)
plotHistPval
, significantTests
## Not run: set.seed(1859)
P <- 1000; P1 <- 100; thresh <- 0.05
pvalues.0 <- setNames(runif(n=P-P1, min=0, max=1),
paste0("null",1:(P-P1)))
pvalues.1 <- setNames(rbeta(n=P1, shape1=1, shape2=10^3),
paste0("alt", 1:P1))
pvalues <- c(pvalues.0, pvalues.1)
pvalues[c(1,10)] <- NA
plotHistPval(pvalues, pi0=(P-P1)/P, freq=TRUE)
qqplotPval(pvalues)
out <- qqplotPval(pvalues, ctl.fwer.bonf=TRUE, ctl.fdr.bh=TRUE, plot.signif=TRUE)
pvalues <- out$pvalues # NA omitted
sum(out$pv.bonf <= thresh)
sum(out$pv.bh <= thresh)
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.