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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|