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

*To follow this exercise exactly be sure you've downloaded the tutorial data set before you begin. If you choose not to download the data set you can follow these instructions on your own data, but you will have to substitute your own specific paths and subject names.

FreeSurfer Tutorial: Group Analysis (25 minutes of exercises)

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

3.1 Using the FreeSurfer Group Descriptor file to create a design matrix BR 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.

3.1.2 Creating a Design Matrix BR

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.

The design matrix is created from the class and variable information in the FSGD file and from the type of design specified when running mri_glmfit (i.e., DODS: different offset different slope; or DOSS: different offset same slope). The design matrix will consist of regressors for intercepts and regressors for slopes. Each class will have an intercept regressor. The intercept regressor is a vector with a 1 for each subject that is a member of the class and 0 otherwise. The slope regressors are handled differently depending upon whether DODS or DOSS is used. For DODS, each class will have a slope regressor for each variable. Like the intercept regressor, the slope regressor for a class will be 0 for subjects not in the class. For subjects in the class, the slope regressor will be the value of the variable. Each variable will have its own set of regressors. For DOSS, each variable will have a single regressor which will be independent of class. This regressor will just be a vector of the variable values listed in the FSGD file. For DODS, the total number of regressors will be Nc*(Nv+1), where Nc is the number of classes and Nv is the number of variables. The first Nc regressors will be the intercepts for each class. The next Nc regressors will be the slopes for each class for the first variable, etc. For DOSS, the total number of regressors will be Nc+Nv. The first Nc regressors will be the intercepts for each class. The next Nv regressors will be the slopes for each variable.

3.1.3 An example FSGDF and design matrices for DODS and DOSS BR This similar FSGD file has two classes (Nc=2) and three variables (Nv=3):

------------------------- cut here ------------------
  GroupDescriptorFile 1
  Class Class1
  CLASS Class2
  Variables             Age  Weight   IQ
  Input subjid1a Class1   10    100   1000
  Input subjid1b Class1   15    150   1500
  Input subjid2a Class2   20    200   2000
  Input subjid2b Class2   25    250   2500
  DefaultVariable Age
------------------------- cut here ------------------

For DODS, the number of regressors is Nc*(Nv+1)=8, and the design matrix would be:

1  0  10  0   100   0   1000    0
1  0  15  0   150   0   1500    0
0  1  0   20   0   200   0    2000
0  1  0   25   0   250   0    2500

For DOSS, the number of regressors is Nc+Nv=5, and the design matrix would be:

1  0  10  100  1000
1  0  15  150  1500
0  1  20  200  2000
0  1  25  250  2500

3.3 Preprocessing steps BR

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, your_average_subject for this example (which is the average you have made for this particular study). In this step you will have to specify your FSGD file, your_fsgd.txt here, your target subject, your_average_subject 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. The following is a sample command line), using mris_preproc, for a thickness study. A similar command will be run in a later exercise, so do not try to execute this sample command line:

mris_preproc --fsgd <fsgd_file.txt> <gd2mtx>\
  --target <average_subject> \
  --hemi ?h \
  --meas thickness \
  --out ?h.fsgd_file.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. Here is a sample command line. Again, a similar command will be run in a later exercise, so do not try to execute this sample command line:

mri_surf2surf --hemi ?h \
  --s <average_subject> \
  --sval ?h.fsgd_file.thickness.mgh \
  --fwhm 10 \
  --tval ?h.fsgd_file.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.

3.4 mri_glmfit BR

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 instance, a contrast vector such as [1 0 0 0 ...] is used to examine the strength of the observed effect from the EV in the first design matrix column. Another contrast vector, [1 -1 0 0 ...], is used to compare the effects between the first two EVs in the design matrix. You can specify your contrast vector as a separate file, which will be read in by mri_glmfit, and used to test your hypotheses.

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. Here is a sample mri_glmfit command. Again, a similar command will be run in a later exercise, so do not try to execute this sample command line:

mri_glmfit --y lh.fsgd_file.thickness.10.mgh \
  --fsgd <fsgd_file.txt> doss\
  --glmdir ?h.fsgd_file.glmdir \
  --surf <average_subject> ?h \
  --C contrast.mat

Where contrast.mat is the contrast vector that you wish to use. 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.

If you had run this command (which you will do in an upcoming exercise) you will have an output directory, ?h.fsgd_file.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  contrast/  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 contrast directory will have the following files:

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

The outputs are as follows:

The following exercise will step you through these 3 commands (mris_preproc, mri_surf2surf, and mri_glmfit) for the tutorial data set:

BR

Note: There are no more exercises for the remainder of this tutorial (just sample command lines).

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