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.4521mm3.
Therefore, in order to calculate the threshold we have to run tstat_threshold.m
in Matlab:
THRESH = TSTAT_THRESHOLD ( SEARCH_VOLUME , FWHM , VOXEL_VOLUME , DF , SIGNIFICANCE )
SEARCH_VOLUME is the volume of the search region in mm3. Default is 1000000, i.e. 1000 cc. The method works well for any value of SEARCH_VOLUME, even zero, which simply gives the threshold for the tstat image at a single voxel. For our example :
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 mm3. Default is 5*5*5=125. The voxel size depends on the field of view and the slice thickness that you choose when you scan your subject. If you didn't write down the information when you scanned, you can find this information with mincinfo. In this case:
voxel_volume = 38.4521;
DF is the degrees of freedom
of the tstat image, as follows:
- For a tstat image from
fmristat, DF is fmristat_df, printed out by fmristat at the end.
- For a tstat image from
multistat, it is the average of the individual fmristat_df's, multiplied
by multistat_df printed out by multistat at the end.
- If multistat tstat is
divided by sdfiles, then DF is now just multistat_df.
- If sdfiles is smoothed
by a Gaussian filter using mincblur -fwhm F, then the degrees of freedom
will be the minimum of:
DF = round( multistat_df * (pi / (2 * log(2)))(3 / 2)*F3 / VOXEL_VOLUME ) and
DF = round( multistat_df * ( 2 * (F/FWHM)2 + 1)(3/2)),
where FWHM is the fwhm of the smoothing kernel applied to the data, the second parameter.
Default is Inf., an infinite degrees of freedom, which assumes the data has a Gaussian rather than a t distribution. The degrees of freedom of the individual sdeffect.mnc files is printed out by fmristat, and for our example is 113.
df = 113;
SIGNIFICANCE is the desired significance in the range from 0 to 0.2. Default is 0.05.
significance = 0.05;
At this point we can compute the threshold:
tstat_threshold( 1200000, 6, 38.4521, 113, 0.05 );
The threshold is 4.9 for
this individual case.
Now that the individual files are analyzed we can combine sessions and subjects. This can be done by using multistat.m 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. It is very important that your data has exactly the same dimensions (slices, rows and columns).
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 fmristat.m, from which we got four files for each of the runs:
*mri_MC_effect.mnc , the
effect file
*mri_MC_rho.mnc, the autocorrelation
file
*mri_MC_sdeffect.mnc, the
standard deviation of the effect file
*mri_MC_tstat.mnc, the t
map file.
In Matlab, multisat, will be used:
df =MULTISTAT( INPUT_FILES_EFFECT , INPUT_FILES_SDEFFECT
, INPUT_FILES_DF , INPUT_FILES_FWHM , X , CONTRAST , OUTPUT_FILE_BASE
, NUM_STATS , FWHM_VARATIO )
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. For that the transformation from native space to the new space must be found, using 'mritotal':
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_effect.mnc
smith_john_19971024_1_105235_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_105235_mri_MC_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_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_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
INPUT_FILES_EFFECT is an array that has as elements the effect files, the * _effect.mnc files
input_file_effect = ['smith_john_19971024_1_105235_mri_MC_effect_tal.mnc '; 'smith_john_19971024_1_105957_mri_MC_effect_tal.mnc '];
INPUT_FILES_SDEFFECT is is an array with the sdeffect files, the standard deviations of the dependent variables, the *_sdeffect.mnc files. If INPUT_FILES_SDEFFECT=[], then INPUT_FILES_SDEFFECT is assumed to be 1 for all voxels. This allows multistat to duplicate DOT for analyzing PET data (with voxel sd, rather than pooled sd).
input_files_sdeffect = ['smith_john_19971024_1_105235_mri_MC_sdeffect_tal.mnc '; 'smith_john_19971024_1_105957_mri_MC_sdeffect_tal.mnc'];
N.B. The filenames should have the same length. This is the way Matlab 4 runs, hopefully in the future, it will be changed, but for that Matlab 5 should be used. For now, just make sure your filenames have the same number of characters.
INPUT_FILES_DF is the average degrees of freedom of the input files, as printed out by fmristat.m. It is used to calculate the degrees of freedom, that will be printed out at the end.
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.
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_Y name minus the .mnc extension, or minus the .mnc.gz extension if it is compressed.
output_file_base = 'smith_john_multi';
NUM_STATS: 1: output T statistic
image, OUTPUT_FILE_BASE_tstat.mnc.
2: in addition, output effect image, OUTPUT_FILE_BASE_effect.mnc
3: in addition, output standard deviation of the effect, OUTPUT_FILE_BASE_sdeffect.mnc.
4: in addition, output standard deviation ratio due to random effect files,
OUTPUT_FILE_BASE_sdfiles.mnc.
Besides those files, the number of degrees of freedom is printed at the end of the program, multistat_df. To calculate the threshold, the degrees of freedom that should be used is the degrees of freedom of INPUT_FILES_SDEFFECT multiplied by multistat_df. If INPUT_FILES_SDEFFECT=[], then the degrees of freedom is just multistat_df.
N.B. tstat=effect/sdeffect.
Default value for NUM_STATS is 1.
num_stats = 4; in order to see not only the t maps but also the effect files, the standard deviation of the effect, and the standard deviation of the files.
For more complicated designs, the design matrix can be changed and the contrast array must change also.
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 FWHM_VARATIO. Default is 15.
For our case, we consider only fixed effects, thus the FWHM_VARATIO will be
fwhm_varatio = Inf;
Again, the threshold for the t maps must be calculated. For the multisession example, the only thing that must be taken in account is the degrees of freedom. The degrees of freedom are calculated by multistat.
Assume, 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 fmristat, one at a time. The effect files and the sdeffect files should be resampled to have the same size. For this example, the files will be transformed into Talairach space. To transform the data, a transformation for each subject must be found. This transformation will be used to transform the respective effect files and sdfiles.
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_effect.mnc
doe_joe_19971120_1_165239_mri_MC_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_sdeffect.mnc
doe_joe_19971120_1_165239_mri_MC_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_effect.mnc
citizen_x_19971125_1_122229_mri_MC_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_sdeffect.mnc
citizen_x_19971125_1_122229_mri_MC_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.
df = MULTISTAT( INPUT_FILES_EFFECT
, INPUT_FILES_SDEFFECT , INPUT_FILES_DF, INPUT_FILES_FWHM , X , CONTRAST
, OUTPUT_FILE_BASE , NUM_STATS , FWHM_VARATIO )
df = multistat( INPUT_FILES_EFFECT , INPUT_FILES_SDEFFECT , X , CONTRAST , OUTPUT_FILE_BASE , NUM_STATS )
input_files_effect = ['smith_john_19971024_1_105235_mri_MC_effect_tal.mnc';
'smith_john_19971024_1_105957_mri_MC_effect_tal.mnc';
'citizen_x_19971125_1_122229_mri_MC_effect_tal.mnc';
'doe_joe_19971120_1_165239_mri_MC_effect_tal.mnc'];
input_files_sdeffect = ['smith_john_19971024_1_105235_mri_MC_sdeffect_tal.mnc';
'smith_john_19971024_1_105957_mri_MC_sdeffect_tal.mnc ;
'citizen_x_19971125_1_122229_mri_MC_sdeffect_tal.mnc';
'doe_joe_19971120_1_165239_mri_MC_sdeffect_tal.mnc'];
input_files_df = 113;
input_files_fwhm = 6;
X = [ 1
1
1
1 ];
contrast = [1];
output_file_base = 'hot_warm_multi';
numstats = 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 |
---|