R/qvalue.R

In QidiPeng/eQTLMAPT: eQTL Mediaiton Analysis with accelerated Permutation Testing approaches

#' qvalue is the function to estimate q values given a vector of p values.
#'
#' The code for this function is extracted from the "qvalue" package written by
#' Alan Dabney and John Storey. We include the functions of qvalue here since it
#' is reported than some operating system cannot isntall it, and we are using
#' version 2.8.0.
qvalue <- function(p, fdr.level = NULL, pfdr = FALSE, lfdr.out = TRUE, pi0 = NULL, ...){
	# Argument checks
	p_in <- qvals_out <- lfdr_out <- p
	rm_na <- !is.na(p)
	p <- p[rm_na]
	if(min(p) < 0 || max(p) > 1){
		stop("p-values not in valid range [0, 1].")
	}else if(!is.null(fdr.level) &&(fdr.level <= 0 || fdr.level > 1)){
		stop("'fdr.level' must be in(0, 1].")
	}

	# Calculate pi0 estimate
	if(is.null(pi0)){
		pi0s <- pi0est(p, ...)
	}else{
		if(pi0 > 0 && pi0 <= 1){
			pi0s = list()
			pi0s$pi0 = pi0
		}else{
			stop("pi0 is not(0,1]")
		}
	}

	# Calculate q-value estimates
	m <- length(p)
	i <- m:1L
	o <- order(p, decreasing = TRUE)
	ro <- order(o)
	if(pfdr){
		qvals <- pi0s$pi0 * pmin(1, cummin(p[o] * m / (i *(1 -(1 - p[o]) ^ m))))[ro]
	}else{
		qvals <- pi0s$pi0 * pmin(1, cummin(p[o] * m / i ))[ro]
	}
	qvals_out[rm_na] <- qvals
	# Calculate local FDR estimates
	if(lfdr.out){
		lfdr <- lfdr(p = p, pi0 = pi0s$pi0, ...)
		lfdr_out[rm_na] <- lfdr
	}else{
		lfdr_out <- NULL
	}

	# Return results
	if(!is.null(fdr.level)){
		retval <- list(call = match.call(), pi0 = pi0s$pi0, qvalues = qvals_out,
						pvalues = p_in, lfdr = lfdr_out, fdr.level = fdr.level,
						significant = (qvals <= fdr.level),
						pi0.lambda = pi0s$pi0.lambda, lambda = pi0s$lambda,
						pi0.smooth = pi0s$pi0.smooth)
	}else{
	retval <- list(call = match.call(), pi0 = pi0s$pi0, qvalues = qvals_out,
					pvalues = p_in, lfdr = lfdr_out, pi0.lambda = pi0s$pi0.lambda,
					lambda = pi0s$lambda, pi0.smooth = pi0s$pi0.smooth)
	}
	class(retval) <- "qvalue"
	return(retval)
}


#' lfdr is a function to estimate the local FDR values from p-values.
#'
#' The code for this function is extracted from the "qvalue" package written by
#' Alan Dabney and John Storey.
lfdr <- function(p, pi0 = NULL, trunc = TRUE, monotone = TRUE,
					transf = c("probit", "logit"), adj = 1.5, eps = 10 ^ -8, ...){
	# Check inputs
	lfdr_out <- p
	rm_na <- !is.na(p)
	p <- p[rm_na]
	if(min(p) < 0 || max(p) > 1){
		stop("P-values not in valid range [0,1].")
	}else if(is.null(pi0)){
		pi0 <- pi0est(p, ...)$pi0
	}
	n <- length(p)
	transf <- match.arg(transf)
	# Local FDR method for both probit and logit transformations
	if(transf == "probit"){
		p <- pmax(p, eps)
		p <- pmin(p, 1 - eps)
		x <- qnorm(p)
		myd <- density(x, adjust = adj)
		mys <- smooth.spline(x = myd$x, y = myd$y)
		y <- predict(mys, x)$y
		lfdr <- pi0 * dnorm(x) / y
	}else{
		x <- log((p + eps) / (1 - p + eps))
		myd <- density(x, adjust = adj)
		mys <- smooth.spline(x = myd$x, y = myd$y)
		y <- predict(mys, x)$y
		dx <- exp(x) /(1 + exp(x)) ^ 2
		lfdr <-(pi0 * dx) / y
	}
	if(trunc){
		lfdr[lfdr > 1] <- 1
	}
	if(monotone){
		o <- order(p, decreasing = FALSE)
		ro <- order(o)
		lfdr <- cummax(lfdr[o])[ro]
	}
	lfdr_out[rm_na] <- lfdr
	return(lfdr_out)
}


#' pi0est is a function to estimates the proportion of true null p-values.
#'
#' The code for this function is extracted from the "qvalue" package written by
#' Alan Dabney and John Storey.
pi0est <- function(p, lambda = seq(0.05,0.95,0.05), pi0.method = c("smoother", "bootstrap"),
									 smooth.df = 3, smooth.log.pi0 = FALSE, ...){
	# Check input arguments
	rm_na <- !is.na(p)
	p <- p[rm_na]
	pi0.method = match.arg(pi0.method)
	m <- length(p)
	lambda <- sort(lambda) # guard against user input

	ll <- length(lambda)
	if(min(p) < 0 || max(p) > 1){
		stop("ERROR: p-values not in valid range [0, 1].")
	}else if(ll > 1 && ll < 4){
		stop(sprintf(paste("ERROR:", paste("length(lambda)=", ll, ".", sep=""),
						 					"If length of lambda greater than 1,",
						 					"you need at least 4 values.")))
	}else if(min(lambda) < 0 || max(lambda) >= 1){
		stop("ERROR: Lambda must be within [0, 1).")
	}
	# Determines pi0
	if(ll == 1){
		pi0 <- mean(p >= lambda)/(1 - lambda)
		pi0.lambda <- pi0
		pi0 <- min(pi0, 1)
		pi0Smooth <- NULL
	}else{
			ind <- length(lambda):1
			pi0 <- cumsum(tabulate(findInterval(p, vec=lambda))[ind]) /(length(p) *(1-lambda[ind]))
			pi0 <- pi0[ind]
##	pi0 <- sapply(lambda, function(l) mean(p >= l) /(1 - l))
		pi0.lambda <- pi0
		# Smoother method approximation
		if(pi0.method == "smoother"){
			if(smooth.log.pi0){
				pi0 <- log(pi0)
				spi0 <- smooth.spline(lambda, pi0, df = smooth.df)
				pi0Smooth <- exp(predict(spi0, x = lambda)$y)
				pi0 <- min(pi0Smooth[ll], 1)
			}else{
				spi0 <- smooth.spline(lambda, pi0, df = smooth.df)
				pi0Smooth <- predict(spi0, x = lambda)$y
				pi0 <- min(pi0Smooth[ll], 1)
			}
		}else if(pi0.method == "bootstrap"){
			# Bootstrap method closed form solution by David Robinson
			minpi0 <- quantile(pi0, prob = 0.1)
			W <- sapply(lambda, function(l) sum(p >= l))
			mse <-(W /(m ^ 2 *(1 - lambda) ^ 2)) *(1 - W / m) +(pi0 - minpi0) ^ 2
			pi0 <- min(pi0[mse == min(mse)], 1)
			pi0Smooth <- NULL
		}else{
			stop('ERROR: pi0.method must be one of "smoother" or "bootstrap".')
		}
	}
	if(pi0 <= 0){
		stop("ERROR: The estimated pi0 <= 0. Check that you have valid p-values or use a different range of lambda.")
	}
	return(list(pi0 = pi0, pi0.lambda = pi0.lambda, lambda = lambda, pi0.smooth = pi0Smooth))
}