Print

---

Hi all,

As pointed out to me by Marko Wilke and Cornelius Werner and other
Germans, that handy script that I created for finding the critical
cluster size threshold was wrong.  It was finding the critical cluster
size in *RESELS*.

Since we don't usually think in terms of cluster size in RESELS, here's
the revised script that produces the critical cluster size threshold
(for a given cluster-defining threshold, and a SPM.mat) in *voxels*.

Sorry for the trouble.

-Tom


%%%%%%%%%%%%%%%%% CorrClusTh.m %%%%%%%%%%%%%%%%%%%%%

function [k,Pc] =CorrClusTh(SPM,u,alpha)
% function [k,Pc] =CorrClusTh(SPM,u,alpha)
% u     - Cluster defining threshold
% SPM   - SPM data structure
% alpha - FWE-corrected level (defaults to 0.05)
%
% Finds the corrected cluster size (spatial extent) threshold for a given
% cluster defining threshold u and FWE-corrected level alpha.
%
%_________________________________________________________________________
% $Id: CorrClusTh.m,v 1.4 2005/04/14 13:42:00 nichols Exp $ Thomas Nichols, Marko Wilke

if nargin<1 | isempty(SPM)
  load(spm_get(1,'SPM.mat', 'Select SPM.mat to check'));
end
df   = [1 SPM.xX.erdf];
STAT = 'T';
R    = SPM.xVol.R;
S    = SPM.xVol.S;
FWHM = SPM.xVol.FWHM;
v2r  = 1/prod(FWHM(~isinf(FWHM)));                      %-voxels to resels

if nargin<2 | isempty(u)
  u = spm_input(['Cluster defining th. {',STAT,' or p value}'],'+0','r',0.001,1);
  if u <= 1; u = spm_u(u,df,STAT); end
end
if nargin<3 | isempty(alpha)
  alpha = 0.05;
end

% Repeatedly run these lines with different cluster size thresholds
% until you get the corrected P you want

Pc = 1;

for k = 1:1000
        Pc = spm_P(1,k*v2r,u,df,'T',R,1,S);
        fprintf('k=%d Pc=%g\n',k,Pc)
        if Pc <= alpha, break; end
end;

fprintf(['  For a cluster-defining threshold of %0.4f the %0.2f-corrected\n'...
        '   cluster size threshold is %d and has size (corrected P-value)' ...
                                     ' %g\n'],u,alpha,k,Pc);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%