I'm developing a Genome-Wide Selection Screening using a statistic called D’2IS or Dvalue. I have an experimental distribution of 70,411,837 Dvalues and a null empirical distribution of the same size. To account for the multiple testing problem I calculated the false discovery rate (FDR) for every experimental Dvalue (Dexp) as

FDR(Dexp) = Snull/Sexp

where Snull is the number of values ≥ Dexp in the null distribution and Sexp is the number of values ≥ Dexp in the experimental distribution (Noble, 2009).

I did not use methods to estimate FDR from p-values such as the Benjamini-Hochberg procedure since I found at bibliography that are only suitable if an analytical null model is available. With only an empirical null model, it will be necessary the estimation of p-values in an intermediate step. This could provide too pessimistic or too optimistic p-values, and in turn, lead to biased FDR values (Strimmer, 2008). Instead, obtaining the FDR directly from both experimental and null distributions (as specified above) seems a faster and more reliable procedure on this situation (Noble, 2009).

For my study three FDR threshold levels were set as FDR <= 0.15, 0.1 and 0.05. The problem comes for the most strict level, since the non-monotone behaviour of the FDR arises around this FDR level:

enter image description here

I know there are methods to achieve monotone FDR values. However, the processing of p-values and/or other parameters subtracted from the data is mandatory in most of them. Therefore:

How I interpret the FDR?

Should I take the minimum Dvalue to achieve FDR <= 0.05 as a threshold? I know that Storey defined the q-value as the minimum FDR attained at or above a given statistic value, but the algorithm proposed by Storey to obtain q-values needs intermediate steps to obtain p-values (Storey, 2003).

Does anyone know a method to achieve monotone FDR values without a theoretical null distribution?


1 Answer 1


I don't think it's an issue that your FDR is non-monotonically decreasing over $D$ values , why should it be? You have an empirical null, and only two parameters you're playing with. It's bound to jump around a little, and the trend is reasonable, I would say to just do as you suggest and take the minimum $D$ s.t. $\text{FDR} \leq 0.05$.

I do have some issue with this approach though, can you provide any other citation for this simple estimate of FDR? I've never seen it before, and it seems a bit too naive to me.

The FDR is defined as the expected proportion of false discoveries over all discoveries. Using some notation from Sun et. al 2006, there are $V$ false positive events (true $H_0$ and declared significant) and $S$ true positive events (true $H_1$ and declared significant).

$$ FDR \equiv E \left[ \frac{V}{V+S} \right] $$

In your method from Noble 2009, he defines the FDR as $s_{null}$ the number of events from the empirical null which rise above the significance threshold divided by the number of events from the real data which rise above the threshold. My interpretation of this is that he is using the empirical null to create an estimator of the number of false positives, therefore saying that $s_{null}$ estimates $V$ (in and of itself I don't necessarily have a problem with this, but I do not think it is a trivial result, and I have not seen it reproduced elsewhere, so I would like to see a formal proof that it is true). Assuming that this is okay (and I'm not totally sure that it is), he goes forward and estimates $V+S$ with the number of positives from the actual data, because this implicitly contains both true and false positives. Dividing these two, he obtains an estimate of the FDR.

There's a problem with this, because FDR is not defined as $\frac{E[V]}{E[V+S]}$, but rather as $E \left[ \frac{V}{V+S} \right]$, which is not a trivial difference. In fact, Storey found that these are only the same when all tests are independent, which is why Sun et. al 2006 only proved stratification for independent tests. So by separately estimating $V$ and $V+S$ (and assuming both are correct!), his FDR estimate is only conservative to the case of independent tests, and will overestimate it. Isn't this the problem you were looking to avoid in the first place?

Personally I think that working with the $P$ values and using an established approach like Storey's $q$ value (as you mentioned) would serve you much better. They do rely on an analytic null distribution, yes, but Strimmer's critism was only of mis-specified null distributions. Do you have evidence that your null $P$ values do not follow $U(0,1)$? If you don't have strong evidence to support the fact that your null is not uniform, I would go with Storey or BH FDR estimates, which will be much more well known and much easier to defend to reviewers (and easier for you too!).

  • $\begingroup$ Thank you very much for your answer, i'm still newbie at statistics so i'm a bit confused. Indeed, this aproximation is naīve and I only found it at Noble (2009), but it seemed easy and suitable for my situation. I don't have any proof that my null p-values do not follow U(0,1), there is a way to find it? As you recomend me, I will implement the Storey aproximation and control the FDR with q-values. $\endgroup$ Aug 25, 2017 at 20:23
  • $\begingroup$ @MiguelRodriguez Not unless, say, there has been previous work on the distribution of null P values for D statistics or something. I would definitely say you're safe with the $q$-value approach, but perhaps you may want to try both approaches and see if the estimates are relatively similar. On a side note, is there anything else that I can add to the answer to better address your question? If it helped you, it's protocol on this site to accept an answer so that others in the future know that it helped you (it also helps with site statistics). $\endgroup$
    – Chris C
    Aug 25, 2017 at 23:29
  • $\begingroup$ Indeed, I will try both approaches to compare the FDR output. As a matter of fact, I have just noticed that Storey's qvalue R package could use a collection of p-values or a collection of test-statistics with corresponding empirical null statistics as input. I am going to read the documentation to check if it is suitable for my situation. $\endgroup$ Aug 27, 2017 at 15:48
  • $\begingroup$ @MiguelRodriguez That's really cool, if Storey ends up supporting your approach I will gladly eat my words. It does seem as though the set up would offer some advantages, perhaps he has thought through it more thoroughly than we have! $\endgroup$
    – Chris C
    Aug 28, 2017 at 14:24
  • $\begingroup$ Extracted from the documentation of qvalue R package: Additionally, the qvalue object can be calculated given a set of empirical null statistics: library(qvalue); data(hedenfalk); obs_stats <- hedenfalk$stat; null_stats <- hedenfalk$stat0; pvalues <- empPvals(stat = obs_stats, stat0 = null_stats); qobj <- qvalue(p = pvalues) It seems that even with empirical null statistics it is mandatory the estimation of p-values in an intermediate step. Anyway, i'll try it. $\endgroup$ Aug 28, 2017 at 19:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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