[wiki:FsTutorial top] | [wiki:FsTutorial previous] | [wiki:FsTutorial/Visualization next]

FreeSurfer Tutorial: Group Analysis

In this tutorial, you will learn how to perform statistical analysis of group surface-based data, including:BR

Assuming that all surface reconstruction has been completed for all subjects in the study, FreeSurfer's mri_glmfit command can be used to perform inter-subject/group averaging and inference on the cortical surface. Mri_glmfit models the data as a linear combination of effects related to variables of interest, confounds and errors, and permits statistical inferences to be made about effects of interest in relation to error variance. It also allows for certain permutation testing and other means for correcting for mutliple comparisons. For group analysis, this technique fits a general linear model (GLM) at each surface vertex to explain the data from all subjects in the study. In this section, a brief overview of linear modeling is presented and mri_glmfit is described for estimating a linear model and testing hypotheses. The modeling overview can be skipped if this material is already familiar. Other software packages have similar types of programs (e.g., FSL's GFEAT).

1.0 Preparing for Group Analysis

For group analysis, you can create an average subject from all the participants in the study. This average will be used as the target subject upon which the results of your group analysis can be output and viewed. To create this average, use make_average_subject. One has already been created for the later exercise, so there is no need to execute this sample command:

make_average_subject --subjects <subj1> <subj2> ...

The default behavior of this script is to create a subject in the $SUBJECTS_DIR named 'average' using each subjects talairach.xfm transform. This behavior can be modified on the command line. You can specify --out your_named_average to change the name of the average subject and --xform talairach.lta (or talairach.m3z) to specify the use of one of the other transforms.

The average subject is created using the processed volumes and surfaces from the set of subjects you specify following the --subjects flag. The make_average_subject command executes both the make_average_volume and make_average_surface subscripts for you.

The distributed example of an average subject can be found in $FREESURFER_HOME/subjects/buckner_data/tutorial_subjs/group_analysis_tutorial called fsaverage

This average subject was created using the required volumes and surfaces of subjects in:

$FREESURFER_HOME/subjects/buckner_data/group_study

2.0 Linear Modeling overview

Linear modeling describes the observed data as a linear combination of explanatory factors plus noise, and determines how well that description explains the data being analyzed.In order to understand how to perform group analysis in FreeSurfer, you need to understand the general linear model (GLM) and how to construct a GLM in matrix notation. You can click [wiki:FsTutorial/GlmReview here] for a review of this material. The notation we use here is:y=X*beta, where y is the vector observed data (eg, thicknesses for each subject at a vertex), X is the known design matrix (eg, gender, age), and beta is the vector of unknown parameter estimates (PEs). The interpretation of the PEs will depend upon how X is constructed. For example, they could be interpeted as a slope indicating the change of thickness with age. The analysis/estimation is then the process of computing beta given the data y and the design matrix X. A Null Hypothesis (H0) is constructed with a contract matrix C. Inferences are drawn by testing whether the value gamma=Cb is zero.

3.0 Using mri_glmfit for estimating the linear model and hypothesis testing

As stated earlier, mri_glmfit performs inter-subject/group averaging and inference on the surface by fitting a linear model at each vertex. The model consists of subject parameters (e.g., age, gender, etc). The model is the same across all vertices, though the fit will probably be different at each vertex. To specify the model, a design matrix that represents the GLM must be created.

Create an FSGD file

The FreeSurfer Group Descriptor File (FSGDF) provides a way to describe a group of subjects and their accompanying data. This can include class membership and other continuous variables, for example gender or age. When it exists, the FSGDF is used by mri_glmfit, tksurfer and tkmedit. The FSGDF is more than just a way to specify the design matrix. It also keeps track of group membership and covariate definitions. This information is then used to help visualize the results. This is not possible when only a design matrix is used.

The following study variables correspond to all the subjects found in $FREESURFER_HOME/subjects/buckner_data/tutorial_subjs/group_analysis_tutorial/ and are presented in a spreadsheet. You'll need to compose an FSGD file with the appropriate variables in order to specify a design matrix that can be used to examine the relationship between a subject's age and cortical thickness. You can read more about the rules for creating an FSGD file [wiki:FsTutorial/CreateFsgdFile here]. To get started, open the text editor of your choice and add individual tags by following the general example and rules described [wiki:FsTutorial/CreateFsgdFile here]. In general, it's a good idea to name the file something intuitive, such as my_gender_age_fsgd.txt. To create your FSGD file you first need to change to the directory you'd like to create the file in:

cd $FREESURFER_HOME/subjects/buckner_data/tutorial_subjs/group_analysis_tutorial/glm

attachment:demoTable2.jpg

Upon completion, save the FSGD file, which can be viewed in the xterm by typing:

cat my_gender_age_fsgd.txt

in the directory where it has been saved.

A correct FSGD file is presented [wiki:FsTutorial/CorrectFsgdFile here] for comparison. It is needed for a later exercise, so create it now if you have not already done so.

Creating a Design Matrix

The FSGDF is specified in the command-line for mri_glmfit with the option --fsgd fname <gd2mtx> where gd2mtx is the method by which the group description is converted into a design matrix. Legal values for gd2mtx are:

If you do not specify one of the above methods, dods will be used by default.

Note: It is not necessary to run mri_glmfit now to create the design matrix, as mri_glmfit will create it for you later in this exercise.

Please click [wiki:FsTutorial/DesignMatrix here] for further explanation of the design matrix and how it is used with the FSGD file.

Preprocessing steps

Once an FSGD file is set up and the average subject is made the preprocessing steps can be followed. The first step will use mris_preproc to assemble your data into a single file in the common surface space, average for this example (which is the average that has been made for this particular study). In this step you will have to specify your FSGD file, gender_age.txt here, your target subject, average here, the hemisphere and measure you are using. You will also name the output file - it's a good idea to use a naming convention that will make it obvious what comparison you are working with. Before running the command you want to be sure your SUBJECTS_DIR is set appropriately and you are in the glm directory you wish to run the analyses in:

setenv SUBJECTS_DIR $FREESURFER_HOME/subjects/buckner_data/tutorial_subjs/group_analysis_tutorial
cd $FREESURFER_HOME/subjects/buckner_data/tutorial_subjs/group_analysis_tutorial/glm

now you can run the mris_preproc command:

mris_preproc --fsgd gender_age.txt doss\
  --target average \
  --hemi lh \
  --meas thickness \
  --out lh.gender_age.thickness.mgh

The next step is to do surface smoothing. Smooth input with a Gaussian kernel with the given full-width/half-maximum (fwhm) specified in mm. For all examples to follow we will use a fwhm = 10mm. To do this mri_surf2surf will be used along with the output from mris_preproc and your average subject.

mri_surf2surf --hemi lh \
  --s average \
  --sval lh.gender_age.thickness.mgh \
  --fwhm 10 \
  --tval lh.gender_age.thickness.10.mgh

You can do the surface smoothing as part of the first step with mris_preproc, but if you do it afterwards as a separate step you can smooth to many different levels without having to rebuild the data each time.

mri_glmfit

mri_glmfit performs the general linear model (GLM) analysis on the volume or the surface. Options include simulation for correction for multiple comparisons, weighted LMS, variance smoothing, PCA/SVD analysis of residuals, per-voxel design matrices, and 'self' regressors. This program performs both the estimation and inference. The framework for testing specific hypotheses is specified in the form of a contrast vector. For each comparison you want to run you will need to design a separate contrast vector. For information on how to set up a contrast vector please click [wiki:FsTutorial/CreateContrastVector here].

To use this in your mri_glmfit command create a file that contains just the one line of your contrast vector using the text editor of your choice. Again, using an appropriate naming convention is a good idea. Make sure this contrast vector has been created for you and is called age.mat:

cd $FREESURFER_HOME/subjects/buckner_data/tutorial_subjs/group_analysis_tutorial/glm
cat age.mat

Create age.mat if it does not exist. It should contain this line:

0 0 1

As an additional exercise, try constructing the contrast vector for testing the same data for difference between males and females independent of age. [wiki:FsTutorial/MaleFemaleContrastSolution Click here for the answer.]

mri_glmfit will take the output from your smoothing step above, your fsgd file, your average subject and the contrast vector as inputs. You will also have to specify a glm directory name, and in this directory all of the outputs will be saved. It is a good idea to use a descriptive name, so you can easily recognize which outputs are in which directory.

mri_glmfit --y lh.gender_age.thickness.10.mgh \
  --fsgd gender_age.txt doss\
  --glmdir lh.gender_age.glmdir \
  --surf average lh \
  --C age.mat

The flag --surf is used to specify that the input has a surface geometry from the hemisphere of the given FreeSurfer subject. If --surf is not specified, then mri_glmfit will assume that the data are volume-based and use the geometry as specified in the header to make spatial calculations.

When this command is finished you will have an lh.gender_age.glmdir. There will be a number of output files in this directory, as well as two other directories. If you did an ls in your glmdir there will be:

ar1.mgh   eres.mgh    mri_glmfit.log  rstd.mgh  age/  y.fsgd
beta.mgh  fsgd.X.mat  rvar.mgh  Xg.dat

the outputs are as follows:

There will be a subdirectory for each contrast that you specify. The name of the directory will be that of the contrast matrix file (without the .mat extension). The age directory will have the following files:

C.dat  F.mgh  gamma.mgh  sig.mgh

The outputs are as follows:

4.0 Using mri_glmfit to correct for multiple comparisons

One method for correcting for multiple comparisons is to perform simulations under the null hypothesis and see how often the value of a statistic from the 'true' analysis is exceeded. This frequency is then interpreted as a p-value which has been corrected for multiple comparisons. This is especially useful with surface-based data as traditional random field theory is harder to implement. This simulator is roughly based on FSLs permuation simulator (randomise) and AFNIs null-z simulator (AlphaSim). Note that FreeSurfer also offers False Discovery Rate (FDR) correction in tkmedit and tksurfer.

The estimation, simulation, and correction are done in three distinct phases:

  1. Estimation: run the analysis on your data without simulation. At this point you can view your results with FDR thresholding in tksurfer. FDR is often conservative relative to cluster-based thresholding.
  2. Simulation: run the simulator with the same parameters as the estimation to get the Cluster Simulation Data (CSD).
  3. Clustering: run mri_surfcluster, passing it the CSD from the simulator and the output of the estimation. These programs will print out clusters along with their p-values.

The estimation step has been described above, using mri_glmfit with no simulation. The simulation step is run using mri_glmfit again, adding in a simulation flag and parameters. If a design is non-orthogonal the permutation simulation can not be run, instead a simple monte carlo simulation can be run. The clustering step is run with mri_surfcluster.

4.1 Simulations BR

The simulation is invoked by calling mri_glmfit and specifying a simulation type and it's associated parameters, with the flag --sim which is to be followed by 4 parameters:

--sim nulltype nsim thresh csdbasename

The first parameter the nulltype, which is the method of generating the null data to be tested. Valid options are:

The next parameter is nsim which corresponds to the number of simulations to run. You can run multiple simulations in parallel, if you have multiple processors, to cut down on processing time.

The next parameter is thresh which corresponds to your threshold and is specified as a -log10(pvalue). Eg, for a p-value threshold of .01, use thresh=2.

The last parameter is csdbasename which corresponds to the base name of the file which will store the Cluster Simulation Data (CSD). Each contrast will get it's own file. When running multiple simulations in parallel be sure to use a unique csdbasename for each run.

Here's a sample command to run the mc-full simulation with 10000 iterations and a p-value threshold of .001:

mri_glmfit --y ?h.fsgd_file.thickness.10.mgh \
  --glmdir ?h.fsgd_file.glmdir \
  --fsgd fsgd_file.txt doss \
  --surf <average_subject> ?h \
  --fwhm 10 --C contrast.mat \
  --sim mc-full 10000 3 ?h.fsgd_file.glmdir/csd1

this will create ?h.fsgd_file.glmdir/csd1-contrast.csd

If you want to split this into multiple runs you could use the following two commands:

mri_glmfit --y ?h.fsgd_file.thickness.10.mgh \
  --fsgd fsgd_file.txt doss\
  --fwhm 10 \
  --glmdir ?h.fsgd_file.glmdir \
  --surf <average_subject> ?h \
  --C contrast.mat \
  --sim mc-full 5000 3 ?h.fsgd_file.glmdir/csd1
mri_glmfit --y ?h.fsgd_file.thickness.10.mgh \
  --fsgd fsgd_file.txt doss\
  --fwhm 10 \
  --glmdir ?h.fsgd_file.glmdir \
  --surf <average_subject> ?h \
  --C contrast.mat \
  --sim mc-full 5000 3 ?h.fsgd_file.glmdir/csd2

which will generate csd1-contrast.csd and csd2-contrast.csd

4.2 Clustering BR

Using the outputs from the estimation step and the simulations, mri_surfcluster (or mri_volcluster) will create two outputs: the summary file with a table of the clusters it found, and an output surface map of the clusters wth the cluster-wise p-value. The sample mri_surfcluster command is:

mri_surfcluster --src ?h.fsgd_file.glmdir/contrast/sig.mgh \
  --csd csd1-contrast.csd \
  --csd csd2-contrast.csd \
  --sum ?h.fsgd_file.glmdir/contrast/sig.cluster.sum \
  --ocp ?h.fsgd_file.glmdir/contrast/sig.cluster.mgh

you can pass all the CSD files that were created through this command by adding as many --csd as you need.

The surfcluster summary file, sig.cluster.sum, will look like this:

#ClusterNo     Max     VtxMax     Size(mm^2)     TalX     TalY     TalZ     CWP     CWPLow     CWPHi
  1            3.561   105964     241.68         31.3     -42.8    26.4     0.06140  0.05830   0.06450
  2           -3.048   86718      9.78           30.4     -66.5    21.6     0.29340  0.28760   0.29920

CWP stands for cluster-wise probability, this is the probability after correction for multiple comparisons. The CWP column is the nomial p-value. CWPLow and CWPHi are the 90% confidence intervals on the p-value. Each cluster gets its own p-value, which depends upon its size. The output surface map, sig.cluster.mgh, will be a map of these clusters with thei CWP. This can be viewed with:

tksurfer <average_subject> ?h inflated \
  -overlay ?h.fsgd_file.glmdir/contrast/sig.cluster.mgh