utility for computing plug-in FDR

Description

utility for computing plug-in FDR

Usage

pifdr( obs, perms, legacy=FALSE, trimToUnit = TRUE, ... ) 

Arguments

Details

Revised 12/30/13 to employ hist() to rapidly bin the permuted values.

Use legacy=TRUE to obtain the approximate implementation, for which the following remarks held: “As currently implemented the algorithm is quadratic in length(obs). While it is possible to get a unique FDR value for every element of obs, an approximate approach yields practically identical precision and by default this will be used for obs with length 2000 or more. In this case, approx is used with rule=2 to interpolate from the grid-based FDR estimates back to the data values.” Additional parameters npts and applier may be supplied if legacy is set to TRUE.

npts defined the number of points spanning the range of obs to be used for a lossy grid-based computation only used if length(obs)>npts.

applier is to be an sapply-like function.

Value

vector of plug-in FDR estimates congruent to obs

References

Hastie Tibshirani and Friedman Elements of Statistical Learning ch 18.7

Examples

set.seed(1234)
op = par(no.readonly=TRUE)
par(mfrow=c(2,2))
X = c(rchisq(30000,1),rchisq(300,10))
Y = rchisq(30300*3,1)
qqplot(Y, X, xlab="null", ylab="observed")
hist(pp <- pifdr(X,Y), xlab="plug-in FDR", main=" ")
library(multtest)
rawp = 1-pchisq(X, 1)
MT <- mt.rawp2adjp(rawp) 
MT2 = MT[[1]][order(MT[[2]]),]
plot(MT2[,"BH"], pp, xlab="BH FDR", ylab="plug-in FDR")
par(op)

Results

R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(GGtools)
Loading required package: GGBase
Loading required package: snpStats
Loading required package: survival
Loading required package: Matrix
Loading required package: data.table
Loading required package: parallel
Loading required package: Homo.sapiens
Loading required package: AnnotationDbi
Loading required package: stats4
Loading required package: BiocGenerics

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, cbind, colnames, do.call, duplicated, eval, evalq,
    get, grep, grepl, intersect, is.unsorted, lapply, lengths, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, rank,
    rbind, rownames, sapply, setdiff, sort, table, tapply, union,
    unique, unsplit

Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: IRanges
Loading required package: S4Vectors

Attaching package: 'S4Vectors'

The following objects are masked from 'package:Matrix':

    colMeans, colSums, expand, rowMeans, rowSums

The following objects are masked from 'package:base':

    colMeans, colSums, expand.grid, rowMeans, rowSums


Attaching package: 'IRanges'

The following object is masked from 'package:data.table':

    shift

Loading required package: OrganismDbi
Loading required package: GenomicFeatures
Loading required package: GenomeInfoDb
Loading required package: GenomicRanges
Loading required package: GO.db

Loading required package: org.Hs.eg.db

Loading required package: TxDb.Hsapiens.UCSC.hg19.knownGene

Attaching package: 'GGtools'

The following object is masked from 'package:stats':

    getCall

> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/GGtools/pifdr.Rd_%03d_medium.png", width=480, height=480)
> ### Name: pifdr
> ### Title: utility for computing plug-in FDR
> ### Aliases: pifdr
> ### Keywords: models
> 
> ### ** Examples
> 
> set.seed(1234)
> op = par(no.readonly=TRUE)
> par(mfrow=c(2,2))
> X = c(rchisq(30000,1),rchisq(300,10))
> Y = rchisq(30300*3,1)
> qqplot(Y, X, xlab="null", ylab="observed")
> hist(pp <- pifdr(X,Y), xlab="plug-in FDR", main=" ")
> library(multtest)
> rawp = 1-pchisq(X, 1)
> MT <- mt.rawp2adjp(rawp) 
> MT2 = MT[[1]][order(MT[[2]]),]
> plot(MT2[,"BH"], pp, xlab="BH FDR", ylab="plug-in FDR")
> par(op)
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>