SPM5 Gem 6: Corrected cluster size threshold

This is SPM5 version of SPM2 Gem 13.

This is a script that will tell you the corrected cluster size threshold for given cluster-defining threshold: CorrClusTh.m

The usage is pretty self explanatory:

 Find the corrected cluster size threshold for a given alpha
 function [k,Pc] =CorrClusTh(SPM,u,alpha,guess)
 SPM   - SPM data structure
 u     - Cluster defining threshold
         If less than zero, u is taken to be uncorrected P-value
 alpha - FWE-corrected level (defaults to 0.05)
 guess - Set to NaN to use a Newton-Rhapson search (default)
         Or provide a explicit list (e.g. 1:1000) of cluster sizes to
         search over.
         If guess is a (non-NaN) scalar nothing happens, except the the
         corrected P-value of guess is printed. 

 Finds the corrected cluster size (spatial extent) threshold for a given
 cluster defining threshold u and FWE-corrected level alpha. 
  

To find the 0.05 (default alpha) corrected cluster size threshold for a 0.01 cluster-defining threshold:

>> load SPM
>> CorrClusTh(SPM,0.01)
  For a cluster-defining threshold of 2.4671 the level 0.05 corrected
  cluster size threshold is 7860 and has size (corrected P-value) 0.0499926
  

Notice that, due to the discreteness of cluster sizes you cannot get an exact 0.05 threshold.

The function uses an automated search which may sometimes fail. If you specify a 4th argument you can manually specify the cluster sizes to search over:

>> CorrClusTh(SPM,0.01,0.05,6000:7000)

  WARNING: Within the range of cluster sizes searched (6000...7000)
  a corrected P-value <= alpha was not found (smallest P: 0.0819589)

  Try increasing the range or an automatic search.

Ooops… bad range.

Lastly, you can use it as a look up for a specific cluster size threshold. For example, how much over the 0.05 level would a cluster size of 7859 be?

>> CorrClusTh(SPM,0.01,0.05,7859)
  For a cluster-defining threshold of 2.4671 a cluster size threshold of
  7859 has corrected P-value 0.050021

Just a pinch!

SPM5 Gem 1: Introduction to SPM5 scripting

Date: Tue, 13 Dec 2005 17:03:13 +0100
From: John Ashburner <john@FIL.ION.UCL.AC.UK>
To: SPM@JISCMAIL.AC.UK
Subject: Re: [SPM] Where can we find some development materials for SPM ?

>    As you know,we usually need to modify the code of SPM to fit our
> problem.but we can not find some relevant development  tutorials. Would you
> please tell me how to learn the framework of SPM step by step ?
>   In addition, I want to know where I can find the details of the SPM
> structure.

It may be easiest to learn by example.  If you want to develop a new 
user-interface for SPM5, then you would create a file called spm_config_*.m, 
similar to the other spm_config.m files (if you strip out the documentation 
parts, you will see that these are actually quite small).  Your spm_config* 
file can then be added to the toolbox subdirectory and accessed through the 
TOOLS pulldown.

The help button allows you to navigate through the documentation of each of 
the Matlab functions, which you may find useful.

For reading and writing images, you would use these functions.
    spm_vol
    spm_slice_vol
    spm_sample_vol

    spm_create_vol
    spm_write_plane
    spm_write_vol

    spm_get_space

There is Matlab help on all these functions.  Alternatively, you could use the 
NIFTI routines directly.  There is no documentation on this, but here are a 
few examples of how you can use them:
======================================

% Example of creating a simulated .nii file.
dat            = file_array;
dat.fname = 'junk.nii';
dat.dim     = [64 64 32];
dat.dtype  = 'FLOAT64-BE';
dat.offset  = ceil(348/8)*8;

% alternatively:
% dat = file_array( 'junk.nii',dim,dtype,off,scale,inter)

disp(dat)

% Create an empty NIFTI structure
N = nifti;

fieldnames(N) % Dump fieldnames

% Creating all the NIFTI header stuff
N.dat = dat;
N.mat = [2 0 0 -110 ; 0 2 0 -110; 0 0 -2 92; 0 0 0 1];
N.mat_intent = 'xxx'
N.mat_intent = 'Scanner';
N.mat0 = N.mat;
N.mat0_intent = 'Aligned';

N.diminfo.slice = 3;
N.diminfo.phase = 2;
N.diminfo.frequency = 2;
N.diminfo.slice_time.code='xxx';
N.diminfo.slice_time.code = 'sequential_increasing';
N.diminfo.slice_time.start = 1;
N.diminfo.slice_time.end = 32;
N.diminfo.slice_time.duration = 3/32;

N.intent.code='xxx' ; % dump possibilities
N.intent.code='FTEST'; % or N.intent.code=4;
N.intent.param = [4 8];

N.timing.toffset = 28800;
N.timing.tspace=3;
N.descrip = 'This is a NIFTI-1 file';
N.aux_file='aux-file-name.txt';
N.cal = [0 1];

create(N); % Writes hdr info

dat(:,:,:)=0;

[i,j,k] = ndgrid(1:64,1:64,1:32);
dat(find((i-32).^2+(j-32).^2+(k*2-32).^2 < 30^2))=1;
dat(find((i-32).^2+(j-32).^2+(k*2-32).^2 < 15^2))=2;


% displaying a slice
imagesc(dat(:,:,12));colorbar

% get a handle to 'junk.nii';
M=nifti('junk.nii');

imagesc(M.dat(:,:,12));
======================================

Best regards,
-John

SPM2 Gem 13: Corrected cluster size threshold

This is a script that will tell you the corrected cluster size threshold for given cluster-defining threshold: corrclusth.m

The usage is pretty self explanatory:

 Find the corrected cluster size threshold for a given alpha
 function [k,Pc] =CorrClusTh(SPM,u,alpha,guess)
 SPM   - SPM data structure
 u     - Cluster defining threshold
         If less than zero, u is taken to be uncorrected P-value
 alpha - FWE-corrected level (defaults to 0.05)
 guess - Set to NaN to use a Newton-Rhapson search (default)
         Or provide a explicit list (e.g. 1:1000) of cluster sizes to
         search over.
         If guess is a (non-NaN) scalar nothing happens, except the the
         corrected P-value of guess is printed. 

 Finds the corrected cluster size (spatial extent) threshold for a given
 cluster defining threshold u and FWE-corrected level alpha. 
  

To find the 0.05 (default alpha) corrected cluster size threshold for a 0.01 cluster-defining threshold:

>> load SPM
>> CorrClusTh(SPM,0.01)
  For a cluster-defining threshold of 2.4671 the level 0.05 corrected
  cluster size threshold is 7860 and has size (corrected P-value) 0.0499926
  

Notice that, due to the discreteness of cluster sizes you cannot get an exact 0.05 threshold.

The function uses an automated search which may sometimes fail. If you specify a 4th argument you can manually specify the cluster sizes to search over:

>> CorrClusTh(SPM,0.01,0.05,6000:7000)

  WARNING: Within the range of cluster sizes searched (6000...7000)
  a corrected P-value <= alpha was not found (smallest P: 0.0819589)

  Try increasing the range or an automatic search.

Ooops… bad range.

Lastly, you can use it as a look up for a specific cluster size threshold. For example, how much over the 0.05 level would a cluster size of 7859 be?

>> CorrClusTh(SPM,0.01,0.05,7859)
  For a cluster-defining threshold of 2.4671 a cluster size threshold of
  7859 has corrected P-value 0.050021

Just a pinch!