qqplotPval: Q-Q plot for p values
In timflutre/rutilstimflutre: Timothee Flutre's personal R code

qqplotPval

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

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.