6
$\begingroup$

I am trying to understand how to correctly apply False Discovery Rate when comparing multiple hypothesis tests. Although I use here R code, my doubts are about the procedure, rather than the programming.

I built a toy model in R made of 10000 hypothesis (e.g. genes expression) made over two 5 samples populations:

set.seed(620)
x = matrix(rnorm(10000*5),nrow=10000)
y = matrix(rnorm(10000*5),nrow=10000)

With these datasets x and y I know that all the null hypothesis are true.

I now evaluate the p-values:

p = sapply(1:10000, function(i) t.test(x[i,],y[i,])$p.val)

As expected the number of p-values below 0.05 (or any other number) is 453, i.e. about 5% false positives as expected.

Next I adjust the p-values using False Discovery Rate adjustment and estimate the q-values:

q = p.adjust(x, method = "fdr")

Now, if I understood correctly, selecting the hypothesis that have a q value of 0.05 one should get 5% of false discoveries (number of false positives divided by the number of discoveries).

The number of hypothesis with q < 0.05 is 0. I think that it might be because, since all the null hypothesis are true, no matter how I choose q the false discoveries will always be 100% among the discoveries (this is also how I explain to myself that most of the q-values are close to 1).

Next I substituted the last one hundred rows of y with numbers sampled by a normal distribution with mean 3, and estimate the p and the q-values:

y[9901:10000,] = rnorm(500, mean=3) 
p = sapply(1:10000, function(i) t.test(x[i,],y[i,])$p.val)
q = p.adjust(p, method = "fdr")

After these modifications, the number of p-values < 0.05 increases to 544 and 98 out of the 100 hypothesis that should be rejected are detected.

However the number of hypothesis having q-values < 0.05 is surprisingly low: only 9. They are all hypothesis that should be rejected, so it seems to me that the false discovery rate has been kept to 0 rather than to 0.05.

If I accept, for example, the hypothesis with q-value = 0.5 I end up with accepting 95 hypothesis. Of these 95, 67 are true discoveries and 28 are false discoveries. Therefore the FDR is of 28/95 = 0.3 and not 0.5 as I would expect.

Is there something I have not understood properly? Why are the result I get so different from the one I would expect theoretically?

$\endgroup$
  • $\begingroup$ You should probably clarify exactly what method you are using to calculate your q-values $\endgroup$ – M. Berk Feb 7 '14 at 13:28
  • 1
    $\begingroup$ thanks! p.adjust with "fdr" method uses the Benjamini–Hochberg procedure. $\endgroup$ – lucacerone Feb 7 '14 at 14:14
  • $\begingroup$ While you are at it, you should also probably clarify whether you are conducting one-sided or two-sided t tests. $\endgroup$ – Alexis Apr 24 '14 at 17:33
3
$\begingroup$

Note that Benjamini-Hochberg controls the false discovery rate to be less than or equal to the specified level, which has been satisfied by all of your scenarios.

In terms of the procedure you have followed, there are no faults there. However, it is worth pointing out that for real gene expression data you may get better results using a specific method such as Storey, John D., and Robert Tibshirani. "Statistical significance for genomewide studies." Proceedings of the National Academy of Sciences 100, no. 16 (2003): 9440-9445. This particular method exploits the fact that in such studies it is virtually impossible that the number of significant genes is zero.

$\endgroup$
  • $\begingroup$ Thanks for the paper! While I find it an interesting read (I work in a biology related field, but not with gene expression data), what bothers me of the example I posted is that in that case everything follows a normal distribution, the tests are independent and yet the fdr only allow to find very few results.. Of course it satisfies the constraints, but still it doesn't have the power I would expect... $\endgroup$ – lucacerone Feb 7 '14 at 15:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.