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(p, 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?