Back to the main page.
Bug 3008 - FDR thresholding wrong
| Status | CLOSED FIXED |
| Reported | 2015-11-19 18:37:00 +0100 |
| Modified | 2016-06-14 16:14:49 +0200 |
| Product: | FieldTrip |
| Component: | core |
| Version: | unspecified |
| Hardware: | All |
| Operating System: | All |
| Importance: | P5 critical |
| Assigned to: | Robert Oostenveld |
| URL: | |
| Tags: | |
| Depends on: | |
| Blocks: | |
| See also: |
- 2015-11-19 18:37:03 +0100
Hi, I realized that the fdr.m file gives a wrong output in the (maybe less frequent) case that the smallest p values are above the FDR threshold. The current script only declares the p values below or equal to the threshold as significant, whereas in fact all p values that are smaller or equal to the largest p value below or equal to this threshold should be declared significant. This is also stated in the original paper that was implemented in the fdr.m file: "find the largest P value that is below the line [threshold]. All voxels with P values less than or equal to this are declared active." An example: p = [0.008 0.008 0.012 0.018 0.03]; q = 0.05; % h=fdr(p,q) gives me the following, when setting a breakpoint at line 39: pi = [0.0044 0.0088 0.0131 0.0175 0.0219]; % slope of the FDR threshold ps = p; % because already sorted from smallest to largest % plot plot(p,'ro'); hold on; plot(pi); % illustrating that the first 3 p values should be significant because 0.012 is the largest p value below the threshold % line 36 gives the following output h = (ps<=pi); h = [0 1 1 0 0]; % this is incorrect, the first sample should also be significant, i.e. 1 % this can be achieved by updating line 36 to the following: h = ps<=(max(ps(ps<=pi))); Cheers, Kathrin
Robert Oostenveld - 2015-11-23 13:18:16 +0100
I think we have been discussing this with various people (including experts) in the past and that there were two slightly different models that could be used in the FDR correction. Let me look up the old discussion thread...
- 2015-11-24 00:05:20 +0100
There are two types of FDR correction (Benjamini-Hochberg & Benjamini-Yekutieli), of which the second is currently implemented. The reported bug would affect both methods. My proposed fix needs to be updated for the case that none of the tests is significant under the FDR correction. The following if-loop replacing the h = ... statement should do (line 36): if ~isempty(ps(ps<=pi)) % no significant tests h = ps<=(max(ps(ps<=pi))); else h = zeros(size(ps)); end
Robert Oostenveld - 2015-11-24 12:29:42 +0100
(In reply to kathrin.muesch from comment #2) thanks, you beat me to finding those references. So it is not related to the difference, which I had suspected. In me reading your suggestion if ~isempty(ps(ps<=pi)) % no significant tests h = ps<=(max(ps(ps<=pi))); else h = zeros(size(ps)); end I would say that it is the same as if any(ps<=pi) h = ps<=(max(ps(ps<=pi))); else h = zeros(size(ps)); end right? but the "if" applies in case there is something significant, and the "else" in case there is not. So I don't understand the placement of your comment (although it is probably not relevant, so don't bother). The "else" clause has to do with no significant tests. Looking up the Genovese paper, figure 1, I see the potential difference between the present FT implementation and the paper. I threshold each sorted p value ("ps") relative to the corresponding "qi" (which gets smaller towards the origin, whereas the paper states that the largest "ps" under the line is to be taken as a fixed threshold. In principle it could have been with the existing implementation that large p values would have reached significance, whereas smaller ones that (after sorting) decrease less fast than the "qi" line (e.g. think of a plateau) would not fall below the "qi" line and hence not have been marked significant. I have fixed it. mac011> svn commit private/fdr.m Sending private/fdr.m Transmitting file data . Committed revision 10935.
- 2015-11-24 16:40:14 +0100
Hi Robert, Your fix is totally appropriate and more elegant than mine. My comment in the code was supposed to go with the else statement. Thanks for fixing it! Kathrin