robust.fdr on p-values

rm(list=ls())
source("robust-fdr.R")

setwd("~... ... CNV/pvalues,qvalues/RTFET")
infiles = list.files(, pattern="txt")
  
for( i in 1:length(infiles)){ 
  tb = read.table(infiles[i], sep="\t", header=F)
  
  tmp1 = tb[tb$V4<1.0, ];
  if( length(tmp1[,1]) >2  ) {
    hist(tmp1$V2)
    ret = robust.fdr(tmp1$V4, sides=1, discrete=T )
    tmp1$q1 = ret$q
    tb$q1 = tmp1$q1[match(tb$V1, tmp1$V1)]
  } else {
    tb$q1 = NA;
  }
  summary(tb)
  tb[tb$q1<0.15 & (!is.na(tb$q1)), ]
  
  outfile = paste("output_robustFDR/", infiles[i], ".csv", sep='')
  write.csv(tb, outfile,, row.names=F, quote=F)
}

Ponds Cheng 06, robust fdr

The robust.fdr routine returns a list with the following components:

p: the vector of p-values provided by the user

fdr: the smoothed local FDR estimates

q: the q-value computed from the smoothed local FDR estimates

cdf: the p-value empirical distribution function

loc.fdr: the unsmoothed local FDR estimates

fp: the estimated number of false positives at the p-value threshold in p

fn: the estimated number of false negatives at the p-value threshold in p

te: the estimated total number of errors in the inference at the p-value threshold in p

pi: the estimated proportion of tests with a true null hypothesis

ord: a vector to order results by ascending p-value


#######################################################
# Function: robust.fdr
# Purpose: Implement robust method of FDR estimation (Pounds and Cheng 2006, Bioinformatics)
# Arguments: p - vector of p-values from the analysis
#            sides - indicate whether p-values are 1-sided (set sides=1) or 2-sided (set sides=2), default=1
#            p2 - for one-sided testing, p-values from testing the "other alternative", default=1-p
#            discrete - indicates (T/F) whether p-values are discrete
#            use8 - indicates whether the constant 8 should be used

dclong.spt m0, pi_0, compare with qvalue()

m0.bdn() example run, fdr

#simulate 1000 p-vlaues from a mixture beta and uniform distribution.

p = rbunif(1000, alpha=1, beta=29, gamma=0.7)

#estimate m0 using the method in Story and Tibshirani (2003).

m0.st(p)

#estimate m0 using the method in Nettleton et al. (2006)

m0.nhcw(p)

# ---------- a real example in Bancroft, Du and Nettleton (2013) --------------

# read in C++ source code

tsf_code = read_code("test_statistics.cpp")

# genrate function for sequential permutation test

# for no location difference between two samples

if(!exists("cxxspt.mean")){

    cxxwrapper(tsf_code, "abs_mean_diff0", "abs_mean_diff", c("n1", "", "int")) -> cxxspt.mean

}

# download data

if(!file.exists('spt_data.rda')){

    download.file('http://dclong.github.io/media/spt/spt_data.rda', 'spt_data.rda')

}

load('spt_data.rda')

cxxspt.mean(leukemia, 10, 1000, n1=5) -> cxxspt.out

#estimate m0 using the method in Bancroft, Du and Nettleton (2013).

m0.bdn(cxxspt.out, 20)

Using a example dataset of 1000 points

From qvalue(), $\pi_0$ = 0.662.

uniform distribution, fishers exact test and pFDR

Carlson, Heckerman, Shani, Microsoft

$\pi_0$ is the mixing parameter for H0.
It can be considered as $\pi_0$ is a prior that H0 is true, and $\pi_1$ is the prior that alternative H1 is true. pFDR is a posterior.

q-value calculation based on fisher.test p-values

https://www.biostars.org/p/46187/

Qvalue, how to choose lambda?

How to choose lambda for Qvalue calculation?



Sometimes, qvalue are smaller than uncorrected pvalues. Based on https://www.biostars.org/p/46187/

First, I think you are correct with being concerned about the estimated q-values and the presence of the warning. In my understanding, p-values adjusted for multiple testing should always be greater-equal to the unadjusted p-values. The warning also indicates that something in the estimation process went wrong. It might have to do with your p-values not fully covering the whole interval [0,1]. If the method tries to find an estimate for the percentage of true null-hypotheses in the data, due to the lack of p-values close or equal to 1, it might not be able to do so. At least this is my interpretation of the warning message.

A very nice manual is available here:
http://genomics.princeton.edu/storeylab/qvalue/manual.pdf

One of the expert replied to my email: 

We usually let the qvalue package select the optimal lambda (see the 2002 pnas paper), as the optimal value minimizes the variance of the computed q-values around the estimated mean. 

Other than that, a lower lambda value would mean pi0 estimates would be higher. lambda 0 would result in the conservative BH adjusted p-values. 

FDR, false discovery rate

http://en.wikipedia.org/wiki/False_discovery_rate

http://www.totallab.com/products/samespots/support/faq/pq-values.aspx