Hi all,
I appear to be getting erroneous output from 3dFDR, but wonder if the problem may be with my syntax or some other issue like roundoff error.
As I had been previously been doing some offline analysis with 1D files (including BH FDR correction), I wanted to verify that my values corroborated with 3dFDR output. So I generated random pvals on a uniform distribution u~[0,0.5] and tried it out. I believe 3dFDR uses the sum(1/n) constant in the FDR equation (for arbitrary pval distributions), but I manually calculated FDR both with and without this constant.
In short, the 3dFDR output did not match my calculated FDR values, and oddly, the 3dFDR output q-values were mostly hovering slightly below 0.5 (the max value on my uniform distribution of p-values).
Any thoughts? Perhaps I am making a simple mistake. First few p-values are below (column 1: "p-values" input, column 2: my FDR vals with no scaling constant, column 3: my FDR w/ constant, column 4: FDR output), see bottom for syntax and my code.
n=2560
pval, fdr1 (no scale), fdr2 (scale), 3dFDR
0.00017073,0.00017073,0.0014384,0.22705
0.00026119,0.00026129,0.0022014,0.22705
0.00032933,0.00032959,0.0027769,0.22705
0.00057553,0.0005762,0.0048546,0.22705
0.0005992,0.00060014,0.0050563,0.22705
0.00065028,0.00065156,0.0054895,0.22705
0.00067891,0.00068051,0.0057334,0.22705
0.00070953,0.00071147,0.0059943,0.22705
0.0015453,0.0015502,0.013061,0.43446
0.0016971,0.001703,0.014348,0.43446
0.0022898,0.0022988,0.019368,0.49432
0.0023171,0.0023271,0.019606,0.49432
0.0028265,0.0028398,0.023926,0.49462
0.0029172,0.0029321,0.024703,0.49462
Command-line:
~/Desktop$ /usr/local/afni/3dFDR -input1D pvals.txt -float -qval -prefix FDR_afni.1D -overwrite -old
++ 3dFDR: AFNI version=AFNI_2008_07_18_1710 (Aug 28 2009) [64-bit]
++ Authored by: B. Douglas Ward
Reading input data: pvals.txt
*+ WARNING: over-writing file FDR_afni.1D
N.B. I tried both 'new' and 'old' methods -- pretty much identical output.
Thanks much,
John Sheppard
Here's my FDR Matlab algorithm for reference, but regardless of my own algorithm, the 3dFDR output definitely seems wrong.
for jj=1:length(pvals)
c(jj)=1/jj;
end
const=sum(c);
pvals = reshape(pvals,length(pvals),1);
pvals = [ [1:length(pvals)]' pvals];
pvals = sortrows(pvals,2);
for jj = 1:length(pvals)
FDRvals(jj,1) = min(1,const*pvals(jj,2)*length(pvals)/(length(pvals)-jj+1)); %B-H, 1995
end
FDRvals = [pvals(:,1) FDRvals]; % concatenate with the indices
FDRvals = sort(FDRvals,1); % put back in the right order
FDRvals = FDRvals(:,2); % get only FDRvals, not indices
end
