FDR correction matlab script

March 14, 2016 40 sec read

The following script is to determine is a channel (or voxel) survives the FDR corrected threshold. Assume you have 48 channels and you already calculated the (uncorrected) p-value of each channel. You put it into an array called p. Now you want to know which channel will survive if you do a FDR correction at q=0.05. You may simply use

fdr0(p, 0.05)

The output is an array of 48 numbers, with 1s in the channel where it survives the FDR threshold.

function y = fdr0(p, q)
% y = fdr0(p, q)
%
% to calculate whether a pvalue survive FDR corrected q
%
% p: an array of p values. (e.g. p values for each channel)
% q: desired FDR threshold (typically 0.05)
% y: an array of the same size with p with only two possible values. 0
% means this position (channel) does not survive the threshold, 1 mean it
% survives
%
% Ref:
% Genovese et al. (2002). Thresholding statistical maps in functional
% neuroimaging using the false discovery rate. Neuroimage, 15:722-786.
%
% Example:
%   y = fdr0(rand(10,1),0.5);
%
% Xu Cui
% 2016/3/14
%

pvalue = p;
y = 0 * p;

[sortedpvalue, sortedposition] = sort(pvalue);
v = length(sortedposition);
for ii=1:v
    if q*ii/v >= sortedpvalue(ii)
        y(sortedposition(ii)) = 1;
    end
end

return;