- T statistic image,
**output_file_base_tstat.mnc**. The degrees of freedom is. Note that tstat=effect/sdeffect.*df* - effect image,
**output_file_base_effect.mnc** - standard deviation of the effect,
**output_file_base_sdeffect.mnc**. - ratio of random effects standard deviation to fixed effects standard deviation,
**output_file-base_sdratio.mnc**.

For clusters, the method of Cao, 1999, Advances in Applied Probability, 31:577-593 is used.

For our example, suppose
we search the tmap for a local maxima inside a region of 1200cc, which
is on average the volume of a brain acquired with EPI. Also, let us say
that the images were filtered in the preprocessing stage with a 6mm Gaussian
filter and the voxel volume is 2.34375mm*2.34375mm*7mm= 38.4521mm^{3.}
Therefore, in order to calculate the threshold we have to run ** tstat_threshold**
in MATLAB. The inputs to this function are:

** search_volume** is the volume
of the search region in mm

search_volume = 1200000;

** fwhm** is the fwhm of
a smoothing kernel applied to the data. Default is 0.0, i.e. no smoothing,
which is roughly the case for raw fMRI data. For PET data, this would be
the scanner fwhm of about 6mm. In the example above, we blurred the images,
thus:

fwhm = 6;

** voxel_volume** is the volume
of one voxel in mm

voxel_volume = 38.4521;

** df** is the degrees of freedom
of the tstat image (DF), either out as the output of

df = 113;

** significance** is the desired
significance in the range from 0 to 0.2. Default is 0.05.

significance = 0.05;

** cluster_threshold** must be less than the calculated threshold,

The outputs of ** tstat_threshold** are

**N.B.** If ** df** < 1000 then the program uses
random numbers to do certain integrals, so the

At this point we can compute the threshold:

The threshold is 4.9 and the cluster size is 418.57mm^{3} for
this individual case. This means that any cluster of neighbouring voxels above a threshold of 3 with a volume larger than 418.57 mm^{3} is significant at P<0.05.

Now that the individual files
are analyzed we can combine sessions and subjects. This can be done
by using ** multistat** in MATLAB.

**Averaging
Data from Multiple Subjects, or Multiple Sessions**

If you have acquired dynamic
data sets for multiple subjects, you will probably want to analyze the
group average, or for the same subject, if your experiment consisted
in multiple runs, or it spread out on different days, you would like to
combine the sessions. This is currently done by resampling the dynamic
data sets for all subjects in Talairach space (on a standardized sampling
grid), or resampling all your data in the same space and generating
an activation map using ** multistat**.

Next, two examples will be discussed:

**A.
Same subject, multiple sessions.**
**B.
Multiple subjects.**

**A.
Same subject, multiple sessions.**

In the example that was described
before, for the subject **Smith John** there are two functional
scans, both of them done using the same protocol and the same design experiment.
We already analyzed both runs using ** fmridesign** and

To calculate the **tsat** for the combined session, ** multisat**, will
be used.

As it can be observed, as an input, the effect files and the sdeffect file are needed. Those must have the same dimensions, thus they will have to be resampled in the same space. For this example they will be transformed in Talairach space.

**N.B.** If the data of all the runs is already in the same space, you don't need to transform it to Talairach space; i.e if you are averaging the runs of one subject, and all the runs were done in the same day, without moving the subject out of the scanner, you can perform multistat in the native space.

To transform the data from native space to the new space, a transformation must
be found. ** mritotal** can be used:

*mritotal* smith_john_19971024_1_2_mri.mnc.gz
smith_john_19971024_1_2_mri_total.xfm

smith_john_19971024_1_2_mri_total.xfm, is the transformation file that will be used for subject Smith John. The effect and sdeffect files from both dynamic scans will be resampled:

*mincresample* smith_john_19971024_1_105235_mri_MC_hmw_effect.mnc

smith_john_19971024_1_105235_mri_MC_hmw_effect_tal.mnc \

-like /data/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz \

-transformation smith_john_19971024_1_2_mri_total.xfm

*mincresample* smith_john_19971024_1_105235_mri_MC_hmw_sdeffect.mnc

smith_john_19971024_1_105235_mri_MC_sdeffect_tal.mnc \

-like /data/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz \

-transformation smith_john_19971024_1_2_mri_total.xfm

*mincresample* smith_john_19971024_1_105957_mri_MC_hmw_effect.mnc

smith_john_19971024_1_105957_mri_MC_effect_tal.mnc \

-like /data/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz \

-transformation smith_john_19971024_1_2_mri_total.xfm

*mincresample* smith_john_19971024_1_105957_mri_MC_hmw_sdeffect.mnc

smith_john_19971024_1_105957_mri_MC_sdeffect_tal.mnc \

-like /data/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz \

-transformation smith_john_19971024_1_2_mri_total.xfm

Instead of using ** mincresample** multiple times, for the same, subject, you could use

Now that your files are in the same space, let's define the inputs for ** multistat**:

** input_files_effect** is an
array that has as elements the effect files, the * _effect.mnc files that were created after running

For our example:

input_file_effect = [ | 'smith_john_19971024_1_105235_mri_MC_hmw_effect_tal.mnc'; |

'smith_john_19971024_1_105957_mri_MC_hmw_effect_tal.mnc' ] ; |

For our example:

input_files_sdeffect = [ | 'smith_john_19971024_1_105235_mri_MC_hmw_sdeffect_tal.mnc'; |

'smith_john_19971024_1_105957_mri_MC_hmw_sdeffect_tal.mnc' ] ; |

**N.B.** The filenames should
have the same length. Just
make sure your filenames have the same number of characters, for that you can just add the correct number of blanks.

** input_files_df** is the average degrees of
freedom of the input files, as printed out by

In the example disscussed before, the degrees of freedom printed was 113, so:

input_files_df = 113;

** input_files_fwhm** is the fwhm in mm of the
original fMRI data. It is only used to calculate the degrees of freedom,
printed out at the beginning.

Default is 6.

input_files_fwhm = 6;

** X** is the design matrix, whose
rows are the files, and columns are the explanatory (independent) variables
of interest.

Default is X=[1; 1; 1; ..1] which just averages the files. If the rank of X equals the number of files, e.g. if X is square, then only a fixed effects analysis is possible. This is done very quickly.

In this example,

X= [ 1

1];

** contrast** is a row vector
of contrasts for the T statistic image. Default is [1 0 ... 0], i.e.
it picks out the first column of X. For our particular situation,

contrast = [1];

** output_file_base**: Base for
output statistics. Default is the first INPUT_FILES_effect name minus the .mnc
extension, or minus the .mnc.gz extension if it is compressed.

output_file_base = 'smith_john_multi';

** which_stats**: logical row vector indicating output statistics by 1.

Columns correspond to:

If the number of columns is less than 4 it is padded with zeros. Default is 1.

For our case:

which_stats = [1 1 1 1];

** fwhm_varatio** is the fwhm in mm of the Gaussian
filter used to smooth the ratio of the random effects variance divided
by the fixed effects variance.

- 0 will do no smoothing, and give a purely
random effects analysis;

- Inf will do complete smoothing to a
global ratio of one, giving a purely fixed effects analysis.

The higher the ** fwhm_varatio**, the higher
the ultimate degrees of freedom of the tstat image (printed out before
a pause), and the more sensitive the test. However too much smoothing
will bias the results. The program pauses to allow inspection of the df;
if it is too low, press control c to cancel (no files are created), and
try increasing

For our case, we consider only fixed effects,
thus the ** fwhm_varatio** will be

fwhm_varatio = Inf;

We can now calculate the statistics:

Again, the threshold for the t maps must be calculated. The way to calculate it is similar to the one run, one subject case , that we disscussed before

For the multisession example, the only thing that must be taken in account is the degrees of freedom. The degrees of freedom are calculated byAssume, that in addition to John Smith, we had two other subjects perform the same task that he carried out in his first dynamic scan. The names of the two other subjects are `Joe Doe' and `Citizen X' and we have two corresponding 3D T1-weighted anatomic volumes:

**citizen_x_19971125_1_2_mri.mnc.gz**

and two dynamic functional
image series (following motion correction):

**citizen_x_19971125_1_122229_mri_MC.mnc**

To analyze the multisubject data
we'll have to discuss :

The way to analyze multiple
subjects data is very similar to single subject, multiple sessions.
The motion corrected files, from each subject, should be analyzed with
** fmridesign** and

The data for Smith John is already resampled, the data from the other two subjects should be resampled also.

mritotal doe_joe_19971120_1_2_mri.mnc.gz joe_doe_19971120_1_2_mri_total.xfm

The doe_joe_19971120_1_2_mri_total.xfm file will be used to transform the effect files and the sdeffect files, of subject Joe Doe from native into Talairach space.

mincresample doe_joe_19971120_1_165239_mri_MC__hmw_effect.mnc

doe_joe_19971120_1_165239_mri_MC_hmw_effect_tal.mnc
\

-like /data/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz \

-transformation doe_joe_19971120_1_2_mri_total.xfm

mincresample doe_joe_19971120_1_165239_mri_MC_hmw_sdeffect.mnc

doe_joe_19971120_1_165239_mri_MC_hmw_sdeffect_tal.mnc
\

-like /data/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz \

-transformation doe_joe_19971120_1_2_mri_total.xfm

Same thing should be done
for **citizen_x_19971125_1_122229_mri_MC.mnc**

mritotal citizen_x_19971125_1_2_mri.mnc.gz
citizen_x_19971125_1_2_mri_total.xfm

mincresample citizen_x_19971125_1_122229_mri_MC__hmw_effect.mnc

citizen_x_19971125_1_122229_mri_MC_hmw_effect_tal.mnc
\

-like /data/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz \

-transformation citizen_x_19971125_1_2_mri_total.xfm

mincresample citizen_x_19971125_1_122229_mri_MC__hmw_sdeffect.mnc

citizen_x_19971125_1_122229_mri_MC_hmw_sdeffect_tal.mnc
\

-like /data/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz \

-transformation citizen_x_19971125_1_2_mri_total.xfm

After resampling all the files, multistat can be used.

input_files_effect = [ | 'smith_john_19971024_1_105235_mri_MC_hmw_effect_tal.mnc'; | |

'smith_john_19971024_1_105957_mri_MC_hmw_effect_tal.mnc'; | ||

'citizen_x_19971125_1_122229_mri_MC_hmw_effect_tal.mnc '; | ||

'doe_joe_19971120_1_165239_mri_MC_hmw_effect_tal.mnc ' | ] ; |

input_files_sdeffect = [ | 'smith_john_19971024_1_105235_mri_MC_hmw_sdeffect_tal.mnc' ; | |

'smith_john_19971024_1_105957_mri_MC_hmw_sdeffect_tal.mnc' ; | ||

'citizen_x_19971125_1_122229_mri_MC_hmw_sdeffect_tal.mnc ' ; | ||

'doe_joe_19971120_1_165239_mri_MC_hmw_sdeffect_tal.mnc ' | ] ; |

As you can see,the file names are padded with blanks, to make the names the equal number of characters. The blanks will be taken out automatically when you use

input_files_df = 113;

input_files_fwhm = 6;

X = [ | ||

1 | ||

1 | ||

1 | ||

1 | ||

] ; |

contrast = [ 1 ];

output_file_base = 'hot_warm_multi';

which_stats = ones ( 1, 4 ) ;

fwhm_varatio = 15;

Because we have multiple subjects, the effects may vary from session to session and from subject to subject. To take this extra source of variability into account, we should do a random effects analysis.

The extra variability is measured by sdfiles,
outputted from ** multistat**. If there is no random effect, then sdfiles should
be roughly 1. If there is a random effect, sdfiles will be greater than
1. To allow for the random effects, we multiply the sdeffect by sdfiles,
and divide the tstat by sdfiles. The only problem is that sdfiles is itself
extremely noisy due to its low degrees of freedom (multistat_df), which
results in a loss of sensitivity. One way around this is to assume that
the underlying effect is fairly smooth, so it can be better estimated by
smoothing sdfiles. This will increase its degrees of freedom, thus increasing
sensitivity. Note however that too much smoothing could introduce bias.
As a guess, try 15mm smoothing, which is acttually the default.

Previous |
Next |
Back to table of contexts |
---|