Thresholding the t-map

stat_threshold returns the threshold for local maxima and the cluster size for tstat-maps from fmrilm or multistat. For the local maxima, two methods are used: random field theory (Worsley et al. 1996, "A Unified Statistical Approach for determining Significant Signals in Images of Cerebral Activation" Human Brain Mapping, 4:58-73), and a simple Bonferroni correction based on the number of voxels in the search region. The final threshold is the minimum of the two. For clusters, the method of Cao, 1999, Advances in Applied Probability, 31:577-593 is used. The stat_threshold command is executed as follows:

[peak_threshold, extent_threshold] = stat_threshold ( search_volume, num_voxels, fwhm, df )


The inputs to this function are: search_volume, num_voxels, fwhm, and df. The p-value is left to the default value of 0.05 and the cluster threshold will also be left unspecified so that we use the default value of 3.16. The values of df and fwhm are obtained from the output of fmrilm. We will now discuss how to estimate the search volume and the number of voxels.

To determine both of these values, we use the function called mask_vol. As input, this function requires a mask_file, which is the original, motion-corrected fMRI data. The output of this function is volume and number_voxels:

[volume,  number_voxels] = mask_vol(mask_file)


Let us say for this particular example, we obtain the following output

volume = 680239
number_voxels = 8503

Next we need the fwhm in mm of a smoothing kernel that was applied to the data. For this we could refer to the fwhm image (i.e. study_subj001_20031019_111740_3_mri_MC_fwhm.mnc) which was output from fmrilm. You could simply view the image with register or MINC visualization tools and obtain a rough mean value of fwhm inside the brain, or your region of interest. You may discover that the mean fwhm value inside the brain will be slightly larger than the blurring value used in the pre-processing step of your analysis (there is more information about this phenomenon in an HBM 2002 abstract by Dr. Worsley). In our example data, we find a fwhm of approximately 8 mm. If you did not blur your fmri data, you would enter a value of 0 in this case.

fwhm = 8;

df is the degrees of freedom of the t-stat image, either out as the output of fmrilm, or multistat. The degrees of freedom of the individual effect files are printed out at the matlab command line when fmrilm finishes executing.
For our example the df is 113.

df = 113;

You may also obtain both the df and the fwhm values from the mincheaders of your statistical output files with the following commands:

To obtain the df:
mincinfo  -attvalue  df:df  <input.mnc>

To obtain the fwhm:
mincinfo  -attvalue  fwhm:fwhm  <input.mnc>

The output of stat_threshold is the peak threshold and cluster extent. At this point we can compute the threshold by calling the stat_threshold command:

[peak_threshold, extent_threshold] = stat_threshold ( search_volume, num_voxels, fwhm, df );

The peak threshold is 4.59 and the extent threshold size is 552 mm3 in this example case. This means that any peak with a t-stat magnitude greater than 4.59 can be considered significant and any cluster of neighbouring voxels above a threshold level of 3.16 with a volume larger than 552 mm3 is significant at p < 0.05.

Once the individual files have been analyzed, we can combine sessions and subjects. This can be done by using multistat.

multistat: Averaging Data from Multiple Subjects, or Multiple Sessions

If you have acquired functional images for multiple subjects, you will want to analyze the group average, or if for the same subject, your experiment consisted of multiple runs, or was acquired on different days, you may want to combine the sessions. All of this is done by first resampling the dynamic data statistical output for all subjects into a common space known as the Talairach space (a standardized sampling grid), or more accurately the MNI space and then running multistat on the resampled statistical output. It is very important that your data have exactly the same dimensions (mincinfo output should be the same across all Tal-space images).

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 subj001 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 fmrilm, from which we obtained the *_mag_ef.mnc, the *_mag_sd.mnc and the *_mag_t.mnc files.

To calculate the t-stat for the combined sessions, multistat, will be used. input_files_effect, input_files_sdeffect, input_files_df, input_files_fwhm, X, contrast, output_file_base, which_stats, fwhm_varatio, are the inputs to multistat, and they have to be defined before the function can be used.

The effect files and the sdeffect file are needed, so you must ensure that the which_stats variable of fmrilm is set appropriately. These images will have to be resampled to Talairach space.

If the data for all the runs of a subject originates from the same scanning session, it is already in the same space and does not need to be transformed 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 and you do not need to resample the effect and sdeffect files. Although the example data set we are using here includes functional runs from the same session and do not need to be resampled to Tal space, we shall nevertheless proceed with the resampling step in order to illustrate the procedure.

To transform the data from different sessions in their native space to Tal, a transformation file (*.xfm) must be found. mritotal and the subject's high resolution anatomical is used to generate this transform file:

mritotal   <anatomical file>  <output transform file name>

mritotal study_subj001_20031019_111740_2_mri.mnc.gz  study_subj001_20031019_111740_2_mri_total.xfm

study_subj001_20031019_111740_2_mri_total.xfm, is the transformation file that will be used for subject subj001. The effect and sdeffect files from both dynamic scans will be resampled. You would enter the following commands in their complete form on the command line according to the following format:

mincresample  <infile>  <outfile>  -transform <transform file>  -like <model>

mincresample study_subj001_20031019_111740_3_mri_MC_hmw_mag_ef.mnc
        study_subj001_20031019_111740_3_mri_MC_hmw_mag_ef_tal.mnc
        -like /data/avgbrain/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz
        -transformation study_subj001_20031019_111740_2_mri_total.xfm

mincresample study_subj001_20031019_111740_3_mri_MC_hmw_mag_sd.mnc
        study_subj001_19971024_1_105235_mri_MC_mag_hmw_sd_tal.mnc
        -like /data/avgbrain/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz
        -transformation study_subj001_20031019_111740_2_mri_total.xfm

mincresample study_subj001_20031019_111740_4_mri_MC_hmw_mag_ef.mnc
        study_subj001_20031019_111740_4_mri_MC_mag_hmw_ef_tal.mnc
        -like /data/avgbrain/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz
        -transformation study_subj001_20031019_111740_2_mri_total.xfm

mincresample study_subj001_20031019_111740_4_mri_MC_hmw_mag_sd.mnc
        study_subj001_20031019_111740_4_mri_MC_mag_hmw_sd_tal.mnc
        -like /data/avgbrain/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz
        -transformation study_subj001_20031019_111740_2_mri_total.xfm

Instead of using mincresample multiple times as we did above, for the same, subject, you could use the resample_tal command to apply the same transformation to resample your images into the new space. Similar to the mincresample command, you will have to provide the transformation and the '-like' file (the model to be used) for the resample_tal command to work. If you type resample_tal -h, on the command-line, you will see all the options for the command.

Note: resample_tal is useful only for the case where you use the same transformation for all the files that you would like to transform at a time.

~mferre/bin/resample_tal  [options]  <file.mnc>  [<file2.mnc> ...]

~mferre/bin/resample_tal  -like /data/avgbrain/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz \
    -transformation   study_subj001_20031019_111740_2_mri_total.xfm \
    study_subj001_20031019_111740_3_mri_MC_hmw_mag_sd.mnc \
    study_subj001_20031019_111740_3_mri_MC_hmw_mag_ef.mnc \
    study_subj001_20031019_111740_4_mri_MC_hmw_mag_ef.mnc \
    study_subj001_20031019_111740_4_mri_MC_hmw_mag_sd.mnc

Now that your files are in the same space, we can define the inputs for multistat. The inputs for multistat are the following:

input_files_effect is an array that has, as elements, the effect files, the *_mag_ef.mnc files, that were created after running fmrilm. The effect filenames should be padded with extra studys if they do not have the same number of characters. For our example:

input_file_effect = ['study_subj001_20031019_111740_3_mri_MC_hmw_mag_ef_tal.mnc';
'study_subj001_20031019_111740_4_mri_MC_hmw_mag_ef_tal.mnc' ] ;

input_files_sdeffect is an array of the sdeffect filenames, *_mag_sd.mnc files:

input_files_sdeffect = ['study_subj001_20031019_111740_3_mri_MC_hmw_mag_sd_tal.mnc';
'study_subj001_20031019_111740_4_mri_MC_hmw_mag_sd_tal.mnc' ] ;

Note: Due to a Matlab syntax restriction, the above variables for the input_files_effect and input_files_sdeffect must have the same number of character spaces between the single quotes. If this is not the case for your list of input files, you should add empty spaces after the last character of the shorter filenames, up to the closing quote mark. Basically, the closing quote character should line up in the same column across the multiple rows of input files listed in these variables.


input_files_df is a row vector of the degrees of freedom of each of the input files, as printed out by fmrilm. It is used to calculate the degrees of freedom which will be printed out at the end of multistat. In the example discussed here, the degrees of freedom printed was 113, so:

input_files_df = [  113  113  ];

Note: you may omit the writing of df values, since they are now included in the mincheaders of your statistical output files. These mincheader entries will be read by multistat automatically, so you may write this variable as follows:

input_files_df = [ ];

input_files_fwhm is the fwhm in mm of the original fMRI data. Here we shall use 8mm. But as is the case for the df values, the fwhm is also read from the mincheaders, so this variable can be left study as well (i.e. input_files_fwhm = [ ]; )

input_files_fwhm = 8;

X is the design matrix, whose rows are the files, and columns are the response variables. 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, which is done very quickly. In this example:

X = [ 1   1 ]' ;

contrast is a row vector of contrasts for the T statistic image. The default is [  1  ]; which simply averages the files referred to by X.

contrast = [  1  ];

Please keep in mind the following WARNING which may apply if you are attempting a subtraction in multistat:

multistat can be very slow if the number of columns of X is more than 1 and less than the number of input files, and input_files_sd is not empty, since it loops over voxels, rather than doing calculations in parallel.

output_file_base: file name base for the output statistics. The default is the first INPUT_FILES_mag_ef name minus the .mnc extension, or minus the .mnc.gz  extension if it is compressed.

output_file_base = 'study_subj001_multi';

which_stats: this is a row vector indicating desired output statistics in a similar method as the fmrilm case that we discussed earlier. You may use character strings to specify the statistics that you want. However, the which_stats strings used in the multistat case differ slightly from the fmrilm case:

In this case, we select the t-stat, effect and sd-effect output:

which_stats = ' _t  _ef  _sd ';

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.

The higher the fwhm_varatio, the higher the ultimate degrees of freedom of the t-stat image, and the more sensitive the test. However too much smoothing will bias the results.

If fwhm_varatio is negative, the absolute value will be taken to be the desired df, and the smoothing value is internally chosen by multistat to get as close to this as possible (if fwhm > 50, fwhm = Inf).

In this example, since we are within a subject we consider only fixed effects, and the fwhm_varatio will be Inf.

fwhm_varatio = Inf;

We can now execute multistat to obtain the output statistic images, and the df of the t-stat test:

df = multistat ( input_files_effect, input_files_sdeffect, input_files_df, input_files_fwhm, X, contrast, output_file_base, which_stats, fwhm_varatio )

The threshold for the t-stat map should be determined. The way to calculate it is similar to the one run, one subject case which was discussed above. For the multistat example, the only change that must be taken in account is the degrees of freedom.

B. Multiple subjects.

Assume, that in addition to subj001, we had two other subjects, subj002 and subj003. We have two corresponding 3D T1-weighted anatomical volumes:

and two dynamic functional image series (following motion correction) in which they performed the same task:

The method to analyze multiple subject data is very similar to single subject, multiple sessions. The motion corrected files, from each subject, should be analyzed with fmridesign and fmrilm, one at a time. The effect files and the sd effect files should be resampled to MNI space. In order to do this, we will need to obtain the transformation file with mritotal:

mritotal  study_subj002_20031120_103700_2_mri.mnc.gz study_subj002_20031120_103700_2_mri_total.xfm

The study_subj002_20031120_103700_2_mri_total.xfm file will be used to transform the effect files and the sdeffect files, of subject subj002 with mincresample.

mincresample study_subj002_20031120_103700_3_mri_MC__hmw_mag_ef.mnc
    study_subj002_20031120_103700_3_mri_MC_hmw_mag_ef_tal.mnc
        -like /data/avgbrain/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz
        -transformation study_subj002_20031120_103700_2_mri_total.xfm

mincresample study_subj002_20031120_103700_3_mri_MC_hmw_mag_sd.mnc
    study_subj002_20031120_103700_3_mri_MC_hmw_mag_sd_tal.mnc
        -like /data/avgbrain/avgbrain1/brain/images/norm_avg_305_mri_2mm_unfilt.mnc.gz
        -transformation study_subj002_20031120_103700_2_mri_total.xfm

These same steps would also be repeated for subject subj003. Once the resampling of all the effect and sdeffect files are performed, multistat can be used.

input_files_effect = ['study_subj001_20031019_111740_3_mri_MC_hmw_mag_ef_tal.mnc';
'study_subj001_20031019_111740_4_mri_MC_hmw_mag_ef_tal.mnc';
'study_subj003_20031122_094500_3_mri_MC_hmw_mag_ef_tal.mnc';
'study_subj002_20031120_103700_3_mri_MC_hmw_mag_ef_tal.mnc;];

input_files_sdeffect = ['study_subj001_20031019_111740_3_mri_MC_hmw_mag_sd_tal.mnc';
'study_subj001_20031019_111740_4_mri_MC_hmw_mag_sd_tal.mnc';
'study_subj003_20031122_094500_3_mri_MC_hmw_mag_sd_tal.mnc';
'study_subj002_20031120_103700_3_mri_MC_hmw_mag_sd_tal.mnc'];

You may have noticed the padding of empty spaces at the end of the shorter filenames above. This is necessary in order for Matlab to process this matrix of input filenames correctly. When entering a matrix of input effect and sd-effect filenames, ensure that the quote characters line up with each other across the matrix rows.

The following definitions of input file degrees of freedom and fwhm could be simply left as an empty matrix, i.e. [ ], since multistat will read the df and fwhm from the mincheader.

input_files_df = [  113  113  113  113  ];

input_files_fwhm = 8;

X = [  1  1  1  1  ]' ;

contrast = [  1  ];

output_file_base = 'grouped_hot_warm_multi';

which_stats = ' _t  _ef  _sd ';

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. Setting the fwhm_varatio to 0 for a completely random effects analysis may introduce bias in the results. Therefore, we perform what is considered a mixed effects analysis and allow multistat to determine the proper amount of smoothing necessary to achieve a df = 100 by setting the fwhm_varatio to -100. This is the default value of fwhm_varatio and can be omitted from the multistat call.

fwhm_varatio = -100;

At this point we repeat the multistat call in the same way that was done previously in part (A) above. Once the multistat results have been obtained we may use visualization tools to display them.





Previous Back to table of contents Next