For internal use in function stp. Computes a BCa confidence
upper bound for the FDR following Algorithm 2 in the vignette.
Usage
bcaFDR(Z, design, th = NULL, B = 100,
lambda = 0.5, PER = FALSE, R = 1000,
gamma = 0.95, Q = NULL, ...)
Arguments
Z
a matrix or data.frame representing genes' expression levels. The rows of
Z correspond to the genes in the experiment, and the columns
correspond to the replicates. Treatment replicates are to the left,
control replicates to the right.
design
a vector of length equal to the number of columns in Z with 1's for
the treatment replicates and 2's for the control replicates (1, …, 1, 2, …, 2).
th
Threshold values for estimating the FDR. If NULL, the values from
abs(ac2(Z,design)) are used.
B
Number of bootstrap or permutation replications for estimating the FDR at
each iteration (as passed from stp).
lambda
Parameter for the estimation of pi0 and the FDR as passed
from stp (see Storey, 2002).
R
Number of bootstrap replications for the computation of the FDR's BCa
confidence upper bound (as passed from stp).
gamma
Confidence level for the FDR's BCa upper confidence bound (as passed
from stp).
PER
If FALSE (default), bootstrap replications are used to estimate
the FDR. If TRUE, permutation replications are used instead
(as passed from stp).
Q
Estimated FDR as returned in object $Q from fdr function
(passed from call to stp). For internal use.
...
additional arguments for parallel computation in boot function as
passed from stp (see stp help page for details).
Value
cbound
BCa upper confidence bound for the FDR for each threshold value in th.
warnings
warning messages generated from use of boot.ci function from package
boot.
Acosta, J. P. (2015) Strategy for Multivariate Identification of
Differentially Expressed Genes in Microarray Data. Unpublished MS thesis.
Universidad Nacional de Colombia, Bogot'a.
Storey, J. D. (2002) A direct approach to false discovery rates.
Journal of the Royal Statistical Society: Series B (Statistical Methodology),
64(3): 479–498.
Efron B. and Tibshirani R. J. (1994) An Introduction to the Bootstrap.
Chapman & Hall/CRC, 1993.
See Also
stp.
Examples
## Single time point analysis for 50 genes with 10 treatment
## replicates and 10 control replicates
n <- 50; p <- 20; p1 <- 10
des <- c(rep(1, p1), rep(2, (p-p1)))
mu <- as.matrix(rexp(n, rate=1))
Z <- t(apply(mu, 1, function(mui) rnorm(p, mean=mui, sd=1)))
### 5 up regulated genes
Z[1:5,1:p1] <- Z[1:5,1:p1] + 5
### 10 down regulated genes
Z[6:15,(p1+1):p] <- Z[6:15,(p1+1):p] + 5
resFdr <- fdr(Z, des)
bca <- bcaFDR(Z, des, Q=resFdr$Q, B=50, R=500)
plot(resFdr$th, resFdr$Q, type="l", col="blue")
lines(resFdr$th, bca$cbound, col="green")
legend(x="topright", legend=c("FDR", "BCa upper bound"),
lty=c(1,1), col=c("blue", "green"))
## Note: Discontinuities in the BCa upper bound are due to warnings
## generated during computations with function code{boot.ci}
## from package code{boot}.
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(acde)
Loading required package: boot
> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/acde/bcaFDR.Rd_%03d_medium.png", width=480, height=480)
> ### Name: bcaFDR
> ### Title: BCa Confidence Upper Bound for the FDR
> ### Aliases: bcaFDR
>
> ### ** Examples
>
> ## Single time point analysis for 50 genes with 10 treatment
> ## replicates and 10 control replicates
> n <- 50; p <- 20; p1 <- 10
> des <- c(rep(1, p1), rep(2, (p-p1)))
> mu <- as.matrix(rexp(n, rate=1))
> Z <- t(apply(mu, 1, function(mui) rnorm(p, mean=mui, sd=1)))
> ### 5 up regulated genes
> Z[1:5,1:p1] <- Z[1:5,1:p1] + 5
> ### 10 down regulated genes
> Z[6:15,(p1+1):p] <- Z[6:15,(p1+1):p] + 5
>
> resFdr <- fdr(Z, des)
> bca <- bcaFDR(Z, des, Q=resFdr$Q, B=50, R=500)
> plot(resFdr$th, resFdr$Q, type="l", col="blue")
> lines(resFdr$th, bca$cbound, col="green")
> legend(x="topright", legend=c("FDR", "BCa upper bound"),
+ lty=c(1,1), col=c("blue", "green"))
> ## Note: Discontinuities in the BCa upper bound are due to warnings
> ## generated during computations with function code{boot.ci}
> ## from package code{boot}.
>
>
>
>
>
> dev.off()
null device
1
>