plotHistPval: P values

In timflutre/rutilstimflutre: Timothee Flutre's personal R code

P values

Description

Plot the histogram of p values with freq=FALSE to imitate figure 1 of Storey & Tibshirani (2003). The x-axis goes from 0 to 1 with equidistant bins of width 0.05, and the y-axis is in the density scale. For a given bin, the density of p values in this bin corresponds to the relative frequency in this bin divided by the bin width, where the relative frequency is the absolute frequency (i.e. raw counts in this bin) divided by the sum of all counts (i.e. over all bins). In the density scale, the area of a given bin is equal to its density times the bin width, thus the area of the whole histogram is equal to 1. Assuming all null hypotheses are true, the "density" histogram represents the density of the discrete version of the Uniform distribution, where the height of each bin is 1, so that the whole histogram area remains equal to 1.

Usage

plotHistPval(
  pvalues,
  breaks = seq(0, 1, 0.05),
  freq = FALSE,
  main = NULL,
  col = "grey",
  border = "white",
  pi0 = NULL,
  ylim = NULL,
  legend.x = "topright",
  verbose = 1
)

Arguments

pvalues	vector of raw p values (missing values will be omitted)
breaks	vector giving the breakpoints between histogram cells (see hist)
freq	as in hist; if TRUE, the histogram graphic is a representation of frequencies, the "counts" component of the result; if FALSE, probability densities, component "density", are plotted (so that the histogram has a total area of one)
main	an overall title for the plot (default: "Density of <nb of non-NA p values> p values")
col	a colour to be used to fill the bars (see hist)
border	the color of the border around the bars (see hist)
pi0	estimate of the proportion of null hypotheses
ylim	optional limits for the y-axis (only if freq=TRUE)
legend.x	x coordinate to be used to position the legend; not shown if set to NULL
verbose	verbosity level (0/1)

Value

invisible output of hist

Author(s)

Timothee Flutre

See Also

qqplotPval

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)
pval.bonf <- p.adjust(pvalues, method="bonferroni")
sum(pval.bonf <= thresh, na.rm=TRUE) # 9
pval.bh <- p.adjust(pvalues, method="BH")
sum(pval.bh <= thresh, na.rm=TRUE) # 104
if(require(qvalue)){
  qv <- qvalue(p=pvalues, fdr.level=thresh, pfdr=TRUE)
  summary(qv)
}

## End(Not run)

timflutre/rutilstimflutre documentation built on Feb. 12, 2025, 11:35 p.m.