FDR "corrected" vs "adjusted" p-values

Demonstration of the difference, following a query about the different results from SPM's spm_P_FDR and FSL's fdr.

See http://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=ind0704&L=fsl&P=22195

Contents

p-values, from random t-values

rand('state', 0); randn('state', 0); % (both used by trnd)
df = 20; % arbitrary (smallish) degrees of freedom
% t non-centrality parameter, 50 null voxels, 50 activated:
delta = [zeros(50,1); ones(50,1)];
t = nctrnd(df, delta);
hist(t)
p = tcdf(-t, df); % right-sided p-values, as SPM

"Corrected" p-values

[Ps Is] = sort(p(:));
S = length(Ps);
cV = 1;
Qs = Ps*S./(1:S)'*cV;
[I Iu] = sort(Is); % indexing to from sorted to original unsorted
clear I;
if ~all(p == Ps(Iu)), error('inverse sorting not successful'), end
Q = Qs(Iu);
plot(Ps, Qs, 'bx-')

"Adjusted" p-values

Code adapted slightly from spm_P_FDR.

Ensures monotonicity of relationship between pUnc and pFDR, see http://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=ind0704&L=fsl&P=22195

P     = zeros(size(p));
for i = 1:length(p)
	% Find PZ(i) in Ps, or smallest Ps(j) s.t. Ps(j) >= PZ(i)
	%---------------------------------------------------------------
	I = min(find(Ps>=p(i)));
	% -"Adjusted" p-values
	%---------------------------------------------------------------
	if isempty(I)
		P(i) = 1;
	else
		P(i) = min(Qs(I:S));
	end
end
plot(Ps, P(Is), 'bx-')

Critical p-value

Pc = spm_uc_FDR(0.5, [1 df], 'P', 1, Ps)

Pc =

    0.2643

Write results

dim = [5 5 4];
V = struct('fname', 'fdr_query_p.nii',...
           'dim',   dim,...
           'mat',   eye(4),...
           'pinfo', [1;0;0],...
           'dt',    [spm_type('float32') spm_platform('BigEnd')],...
           'n',     [1 1]...
           );
spm_write_vol(V, reshape(p, dim));
V.fname = 'fdr_query_P.nii';
spm_write_vol(V, reshape(P, dim));
V.fname = 'fdr_query_Q.nii';
spm_write_vol(V, reshape(Q, dim));
V.fname = 'fdr_query_mask.nii';
spm_write_vol(V, ones(dim));

Use FSL's fdr

[status stdout] = system('fsl'); % (ignore me...)
if status ~= 0, error('FSL not working'), end
cmd = 'setenv FSLOUTPUTTYPE NIFTI; ';
cmd = [cmd 'fdr -i fdr_query_p.nii -m fdr_query_mask.nii -o fdr_query_fslQ -q 0.5'];
[status stdout] = system(cmd);
fslPc = str2double(regexp(stdout, '0\.\d+', 'match'))
fslQ = spm_read_vols(spm_vol('fdr_query_fslQ.nii'));
fslQ = fslQ(:);

fslPc =

    0.2643

Compare

errPc = abs(Pc - fslPc)
maxQerr = max(abs(Q-fslQ)) % FSL's fdr p-values (q-rate) correspond to Q
plot(1:100,p, 1:100,P, 1:100,fslQ);
legend('Uncorrected', 'SPM P', 'FSL Q')

errPc =

   3.4523e-07


maxQerr =

   5.7571e-08

Note Q is conservative in the following sense:

all(Q >= P)

ans =

     1

Ged Ridgway (drc.spm@gmail.com), with thanks to Tom Nichols.

Published with MATLAB® 7.1