estimate.m0=function(p, B = 20) { # #This function estimates the number of true null hypotheses given a vector of p-values #using the method of Nettleton et al. (2006) JABES 11, 337-356. #The estimate obtained is identical to the estimate obtained by the iterative #procedure described by Mosig et al. Genetics 157:1683-1698. #The number of p-values falling into B equally sized bins are counted. #The count of each bin is compared to the average of all the bin counts associated #with the current bins and all bins to its right. Working from left to right, #the first bin that has a count less than or equal to the average is identified. #That average is multiplied by the total number of bins to obtain an estimate of m0, #the number of tests for which the null hypothesis is true. # m <- length(p) m0 <- m bin <- c(-0.1, (1:B)/B) bin.counts=rep(0,B) for(i in 1:B){ bin.counts[i]=sum((p>bin[i])&(p<=bin[i+1])) } tail.means <- rev(cumsum(rev(bin.counts))/(1:B)) temp <- bin.counts - tail.means index <- min((1:B)[temp <= 0]) m0 <- B * tail.means[index] return(m0) } jabes.q=function(p,B=20) { # #This function computes q-values using the approach of Nettleton et al. #(2006) JABES 11, 337-356. # #Author: Dan Nettleton # m = length(p) m0=estimate.m0(p,B) k = 1:m ord = order(p) p[ord] = (p[ord] * m0)/(1:m) qval = p qval[ord]=rev(cummin(rev(qval[ord]))) return(qval) } ub.mix= function(x,start=c(.6,1,4)){ # #The method of Allison et al. (2002) CSDA is used to model the #p-value distribution as a mixture of a uniform and another beta #distribution. The mixing proportion on the uniform (pi0) is an #estimate of the proportion of true null hypotheses. # #This code is adapted from code provided by Gary Gadbury. # #The starting values were provided by Gary's group, but they #may need to be changed dependending on the problem. #Always check the fit of the estimated distribution using #the function plot.ub.mix. # #The function returns a vector that contains estimates of #pi0, alpha, and beta in that order. # mix.obj1<-function(p,x){ # This object is the objective function for a mixture of a uniform # and one beta distribution # e=p[1] + (1-p[1])*(gamma(p[2]+p[3]) /(gamma(p[2])*gamma(p[3])))*x^(p[2]-1)*(1-x)^(p[3]-1) sle=-sum(log(e)) return(sle) } return(optim(start,mix.obj1,lower=c(.001,.01,.01), upper=c(1,170,170),method="L-BFGS-B",x = x)\$par) } plot.ub.mix= function (p,parms=NULL,nclass=20) { # #This function plots the estimated uniform-beta mixture density #over a histogram of p-values. # #The first argument p is a required vector of p-values. #The second argument is a vector which contains estimates #of pi0, alpha, and beta. If no vector is supplied, one #will be estimated using the function ub.mix. The third #argument specifies the number of bins for the histogram. # if(is.null(parms)){ parms=ub.mix(p) } x=0:1000/1000 hist(p,probability=T,col=4,nclass=nclass) box() lines(x,parms[1]+(1-parms[1])*dbeta(x,parms[2],parms[3]),col=2,lwd=2) lines(c(0,1),rep(parms[1],2),lty=2,col=2,lwd=2) } ppde= function (p,parms=NULL) { # #p is a vector of p-values #parms is a vector containing estimates of pi0, alpha, and beta #from a uniform-beta mixture model fit the the p-values. #If parms is not specified, it will be estimated using ub.mix. # if(is.null(parms)){ parms=ub.mix(p) } pi0=parms[1] pi1=1-pi0 a=parms[2] b=parms[3] return(1/(pi0/(pi1*dbeta(p,a,b))+1)) } q.value= function(p, lambda=seq(0,0.95,0.05)) { # #This is Storey and Tibshirani's (PNAS, 2003) default method #for determining the q-values. # #Code was originally obtained from John Storey's Web site. #It has been edited slightly. # m<-length(p) pi0 <- rep(0,length(lambda)) for(i in 1:length(lambda)) { pi0[i] <- mean(p >= lambda[i])/(1-lambda[i]) } spi0 <- smooth.spline(lambda,pi0,df=3) pi0 <- max(predict(spi0,x=1)\$y,0) pi0 <- min(pi0,1) u <- order(p) v <- rank(p) qvalue <- pi0*m*p/v qvalue[u[m]] <- min(qvalue[u[m]],1) for(i in (m-1):1) { qvalue[u[i]] <- min(qvalue[u[i]],qvalue[u[i+1]],1) } return(qvalue) } bh.fdr= function(p) { # #This function computes q-values using Benjamini and Hochberg's (1995) #approach for controlling FDR. # m = length(p) k = 1:m ord = order(p) p[ord] = (p[ord] * m)/(1:m) qval = p qval[ord]=rev(cummin(rev(qval[ord]))) return(qval) }