Homec4science
Diffusion Slepians (master) Browse

Slepi/Code/Utilities/spm12/spm_uc_clusterFDR.m
8a186cd82567master

Last Change
Raw File

spm_uc_clusterFDR.m
View Options

function [u, Ps, ue] = spm_uc_clusterFDR(q,df,STAT,R,n,Z,XYZ,V2R,ui,G)
% Cluster False Discovery critical extent threshold
% FORMAT [u, Ps, ue] = spm_uc_clusterFDR(q,df,STAT,R,n,Z,XYZ,ui[,G])
%
% q - prespecified upper bound on False Discovery Rate
% df - [df{interest} df{residuals}]
% STAT - statistical field
% 'Z' - Gaussian field
% 'T' - T - field
% 'X' - Chi squared field
% 'F' - F - field
% R - RESEL Count {defining search volume}
% n - conjunction number
% Z - height {minimum over n values}
% or mapped statistic image(s)
% XYZ - locations [x y x]' {in voxels}
% or vector of indices of elements within mask
% or mapped mask image
% V2R - voxel to resel
% ui - feature-inducing threshold
% G - patch structure (for surface-based inference)
%
% u - critical extent threshold
% Ps - sorted p-values
% ue - critical extent threshold for FWE
%__________________________________________________________________________
%
% References
%
% J.R. Chumbley and K.J. Friston, "False discovery rate revisited: FDR and
% topological inference using Gaussian random fields". NeuroImage,
% 44(1):62-70, 2009.
%
% J.R. Chumbley, K.J. Worsley, G. Flandin and K.J. Friston, "Topological
% FDR for NeuroImaging". NeuroImage, 49(4):3057-3064, 2010.
%__________________________________________________________________________
% Copyright (C) 2009-2013 Wellcome Trust Centre for Neuroimaging
% Justin Chumbley & Guillaume Flandin
% $Id: spm_uc_clusterFDR.m 5809 2013-12-20 14:30:22Z guillaume $
% Read statistical value from disk if needed
%--------------------------------------------------------------------------
if isstruct(Z)
Vs = Z;
Vm = XYZ;
[Z, XYZmm] = spm_read_vols(Vs(1),true);
for i=2:numel(Vs)
Z = min(Z, spm_read_vols(Vs(i),true));
end
Z = Z(:)';
XYZ = round(Vs(1).mat \ [XYZmm; ones(1, size(XYZmm, 2))]);
try
M = spm_vol(spm_file(Vs(1).fname,'basename','mask'));
I = spm_get_data(M,XYZ,false) > 0;
catch
I = ~isnan(Z) & Z~=0;
end
XYZ = XYZ(1:3,I);
Z = Z(I);
if isstruct(Vm)
Vm = spm_get_data(Vm,XYZ,false) > 0;
end
XYZ = XYZ(:,Vm);
Z = Z(:,Vm);
end
% Threshold the statistical field
%--------------------------------------------------------------------------
XYZ = XYZ(:,Z >= ui);
% Extract size of excursion sets
%--------------------------------------------------------------------------
if nargin == 9
N = spm_clusters(XYZ);
else
tmp = false(size(G.vertices,1),1);
tmp(XYZ(1,:)) = true;
N = spm_mesh_clusters(G,tmp);
N = N(XYZ(1,:));
end
if ~isempty(N)
N = histc(N,(0:max(N))+0.5);
end
N = N(1:end-1);
N = N .* V2R;
% Compute uncorrected p-values based on N using Random Field Theory
%--------------------------------------------------------------------------
Ps = zeros(1,numel(N));
Pk = zeros(1,numel(N));
ws = warning('off','SPM:outOfRangePoisson');
for i = 1:length(N)
[Pk(i), Ps(i)] = spm_P_RF(1,N(i),ui,df,STAT,R,n);
end
warning(ws);
[Ps, J] = sort(Ps, 'ascend');
S = length(Ps);
% Calculate FDR inequality RHS
%--------------------------------------------------------------------------
cV = 1; % Benjamini & Yeuketeli cV for independence/PosRegDep case
Fi = (1:S)/S*q/cV;
% Find threshold
%--------------------------------------------------------------------------
I = find(Ps <= Fi, 1, 'last');
if isempty(I)
u = Inf;
else
u = N(J(I)) / V2R;
end
% As we're there, also determine the FWE critical extent threshold
%--------------------------------------------------------------------------
if nargout == 3
[Pk, J] = sort(Pk, 'ascend');
I = find(Pk <= q);
if isempty(I)
ue = Inf;
else
ue = min(N(J(I))) / V2R;
end
end
c4science · Help