qqplotPval: Q-Q plot for p values

View source: R/stats.R

qqplotPvalR Documentation

Q-Q plot for p values

Description

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.

Usage

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
)

Arguments

pvalues

vector of raw p values (missing values will be omitted)

use.density

if TRUE, uses smoothScatter

nrpoints

if use.density=TRUE, number of points to be superimposed on the density image

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 qvalue()

lim.custom

custom limit on p values (optional)

lty.custom

line style for the custom limit (compulsory if lim.custom is supplied)

lgd.custom

legend for the custom limit (compulsory if lim.custom is supplied)

plot.signif

show line(s) corresponding to the significance threshold (compulsory if lim.custom is supplied)

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)

Value

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

Author(s)

Timothee Flutre (inspired by an anonymous comment at http://gettinggeneticsdone.blogspot.fr/2009/11/qq-plots-of-p-values-in-r-using-ggplot2.html)

See Also

plotHistPval, significantTests

Examples

## 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)

timflutre/rutilstimflutre documentation built on Aug. 18, 2024, 7:43 p.m.