pr_fdr

 False Discovery Rate: adaptive algorithm.
 FORMAT H = pr_fdr(p,alfa,No,sortf)
 
 Input:
 p             - pvalues
 alfa          - FDR
 No            - number of Null Hypotheses
 sortf         - is 1 if p is already sorted ascending
 
 Output:
 H             - vector of rejections 
                 1 hypothesis rejected
                 0 otherwise

 References:
 
 Benjamini Y, Hochberg Y (2000) "On the adaptive control of the false
 discovery rate in multiple testing with independent statistics", Journal of
 Educational and Behavioral Statistics, 25:(1),60-83 (2000).
 
 Benjamini Y, Hochberg Y (1995) "Controlling the False Discovery Rate:
 a Practical and Powerful Approach to Multiple Testing", Journal
 of the Royal Statistical Society B, 57 289-300.


 Federico E. Turkheimer 16/2/2001

 $Id: pr_fdr.m,v 1.1 2005/06/01 09:24:23 matthewbrett Exp $

0001 function H = pr_fdr(p,alfa,No,sortf)
0002 % False Discovery Rate: adaptive algorithm.
0003 % FORMAT H = pr_fdr(p,alfa,No,sortf)
0004 %
0005 % Input:
0006 % p             - pvalues
0007 % alfa          - FDR
0008 % No            - number of Null Hypotheses
0009 % sortf         - is 1 if p is already sorted ascending
0010 %
0011 % Output:
0012 % H             - vector of rejections
0013 %                 1 hypothesis rejected
0014 %                 0 otherwise
0015 %
0016 % References:
0017 %
0018 % Benjamini Y, Hochberg Y (2000) "On the adaptive control of the false
0019 % discovery rate in multiple testing with independent statistics", Journal of
0020 % Educational and Behavioral Statistics, 25:(1),60-83 (2000).
0021 %
0022 % Benjamini Y, Hochberg Y (1995) "Controlling the False Discovery Rate:
0023 % a Practical and Powerful Approach to Multiple Testing", Journal
0024 % of the Royal Statistical Society B, 57 289-300.
0025 %
0026 %
0027 % Federico E. Turkheimer 16/2/2001
0028 %
0029 % $Id: pr_fdr.m,v 1.1 2005/06/01 09:24:23 matthewbrett Exp $
0030 
0031 if nargin < 1
0032   error('Need p values');
0033 end
0034 if nargin < 2
0035   alfa = 0.05;
0036 end
0037 n    = length(p);
0038 if nargin < 3
0039   No = n;
0040 elseif No > n
0041   No = n;
0042 end
0043 if nargin < 4
0044   sortf = 0;
0045 end
0046 if ~sortf
0047   q = sort(p(:));
0048 else
0049   q = p(:);
0050 end
0051 
0052 H = zeros(size(p));
0053 
0054 idx = 1:n;
0055 brk = max(find(q<=idx'*alfa/No));
0056 if isempty(brk), return, end
0057 idx = 1:(brk-1);
0058 if any(q(idx)<=idx'*alfa/n)
0059   th = q(brk);
0060   H(p<=th) = 1;
0061 end
0062 
0063 return