R/FDRplot.R

In fdrci: Permutation-Based FDR Point and Confidence Interval Estimation

#' Plot results of FDR table generated by fdrTbl()
#'
#' This function plots FDR point and CI estimates over a sequence of 
#' possible significance thresholds. Results from fdrTbl() can be plotted directly as input to FDRplot.
#'
#' @param plotdat a table that is returned from fdrTbl(), or results formated in the same way.
#' @param lowerbound -log10(p-value) lower bound for the x-axis of the plot.
#' @param upperbound -log10(p-value) upper bound for the x-axis of the plot.
#' @param ymax upper limit for range of the y-axis.
#' @param annot annotation text to be added to plot area.
#' @param xpos x-axis position for annot
#' @param ypos y-axis position for annot
#' @author Joshua Millstein, \email{joshua.millstein@@usc.edu}
#' @references Millstein J, Volfson D. 2013. Computationally efficient permutation-based confidence interval estimation for tail-area FDR. Frontiers in Genetics | Statistical Genetics and Methodology 4(179):1-11.
#' @importFrom ggplot2 ggplot geom_ribbon aes scale_x_continuous scale_y_continuous theme_bw geom_line xlab ylab labs unit theme geom_point geom_text annotate
#' @importFrom grDevices pdf dev.off
#' @return ggplot2 object
#'
#' @author Joshua Millstein
#' @references Millstein J, Volfson D. 2013. Computationally efficient
#' permutation-based confidence interval estimation for tail-area FDR.
#' Frontiers in Genetics | Statistical Genetics and Methodology 4(179):1-11.
#' @examples
#' ss = 100
#' nvar = 100
#' X = as.data.frame(matrix(rnorm(ss*nvar),nrow=ss,ncol=nvar))
#' e = as.data.frame(matrix(rnorm(ss*nvar),nrow=ss,ncol=nvar))
#' Y = .1*X + e
#' nperm = 10
#' 
#' myanalysis = function(X,Y){
#' 	ntests = ncol(X)
#' 	rslts = as.data.frame(matrix(NA,nrow=ntests,ncol=2))
#' 	names(rslts) = c("ID","pvalue")
#' 	rslts[,"ID"] = 1:ntests
#' 	for(i in 1:ntests){
#' 		fit = cor.test(X[,i],Y[,i],na.action="na.exclude",
#' 			alternative="two.sided",method="pearson")
#' 		rslts[i,"pvalue"] = fit$p.value
#' 	}
#' 	return(rslts)
#' } # End myanalysis
#' 
#' # Generate observed results
#' obs = myanalysis(X,Y)
#' 
#' # Generate permuted results
#' perml = vector('list',nperm)
#' for(perm in 1:nperm){
#' 	X1 = X[order(runif(nvar)),]
#' 	perml[[perm]] = myanalysis(X1,Y)
#' }
#' 
#' # FDR results table
#' myfdrtbl = fdrTbl(obs$pvalue,perml,"pvalue",nvar,0,3)
#' # Plot results
#' FDRplot(myfdrtbl,0,3,annot="A. An Example")
#' 
#' @export
FDRplot <-
function(plotdat,lowerbound,upperbound,ymax=1, annot="",xpos=.8, ypos=.8){
	# plotdat is a table that is returned from fdrTbl()
	# pname contains a string that is the name of the permutation p-value column
	# lowerbound and upperbound are -log10(p-value) bounds for the x-axis of the plot
	# mn = main title

	##### myfdrplot.pdf
	threshold <- ll <- ul <- fdr <- S <- NULL
	p1 = ggplot2::ggplot(plotdat, aes(x=threshold)) +
		ggplot2::geom_ribbon(aes(ymin=ll, ymax=ul), alpha=.4, fill="darkseagreen2") + 
		ggplot2::scale_x_continuous(limits = c(lowerbound, upperbound)) + 
		ggplot2::scale_y_continuous(limits = c(0, ymax)) +
		ggplot2::theme_bw(base_size = 16) + 
		ggplot2::geom_line(aes(y = fdr), color="aquamarine4", size=1.5) +
		ggplot2::xlab(expression(paste('detection threshold (-',log[10],'(p-value))'))) +
		ggplot2::ylab("FDR") + 
		ggplot2::labs(title="") +
		ggplot2::theme(plot.margin=unit(c(0,20,10,10),"pt")) +
		ggplot2::geom_point(aes( y=fdr)) + 
		ggplot2::geom_text(aes(y=fdr, label=S), hjust=-.2, vjust=-.2) +
		ggplot2::annotate("text", x=xpos, y=ypos, label=annot, size=6)
	p1	
	return(p1)
} # End FDRplot function