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

FreeSurfer Tutorial: Group Analysis (25 minutes of exercises)

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:

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

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.

For group morphometric analysis, the observed data is comprised of a set of surface measures (such as cortical thickness) at each vertex in a surface model, for each subject in the group. This data can be organized as a set of vectors, each associated with a different vertex in the surface model, and containing a surface measurement for every subject in the group at the corresponding vertex.

First, a linear model must be designed to include all explanatory variables (EVs) that may account for each vector's values. A simple linear model is given by y=a*x+b+e, where the observed data y is a one-dimensional vector of surface measures -- one measurement per subject at a vertex; x is a one-dimensional vector containing a variable, such as age, describing each subject; a is the parameter estimate (PE) for x, for instance the value that a subject's age must be multiplied by to fit the data in y; b is a constant, and in this example, would correspond to the baseline measurement present in the data; and e is the error in the model fitting. If an additional explanatory variable is added to explain the observed data, the model would be given as y=a1*x1+a2*x2+b+e, containing two different model waveforms, a1*x1 and a2*x2, corresponding to two variables, such as age and gender, describing all subjects in the study.

2.1 Estimation overview BR Once the model is specified, an estimation step follows, in which the model is fit to each vertex's vector separately; no interactions between vertices are taken into account in the examples presented here. This step generates the estimate of the "goodness of fit" of each of the explanatory variables to each vector of surface measurements. Thus if a particular vertex responds strongly to the explanatory variable x1, a large value for a1 will be produced by model-fitting; if the data appears unrelated to x2 then a2 will have a very small value.

This kind of linear modeling is commonly expressed in matrix notation, where the the matrix X contains all the explanatory variables (designed effects and confounds) in the model, and the matrix A contains all the PEs. The matrix X is also commonly called the design matrix and it can be user-specified in FreeSurfer in the form of an FSGD (FreeSurfer Group Descriptor) file, as the exercises below illustrate. Each column of X corresponds to a different explanatory variable (also called a regressor or a covariate). As typically formulated and solved, the estimation step produces a set of estimates of the PEs, which in turn are used in hypothesis testing.

2.2 Inference overview BR Estimates of the PEs can be converted into statistical parametric maps, which are commonly visualized as a color-coded surface overlay. The overlay assigns each vertex a value based on the likelihood that the null hypothesis is false at that vertex. A linear combination of estimates of PEs is used to encode the particular hypothesis of interest. This encoding is accomplished with a user-specified contrast vector, which assigns a contrast weight to each column of the design matrix. A simplest example of a contrast vector that tests the null hypothesis for the explanatory variable associated with the first design matrix column would be[ 1 0 0 0...]. To compute this particular contrast at each vertex, the PE value associated with the first design matrix column at that vertex is divided by the error in its estimate, yielding a t-value. The t-value provides a good measure confidence in the estimate of the PE value, and can be converted into a probability (P) or Z statistic at that vertex via a standard statistical transformation. T, P and Z values all convey the same information about how significantly the observed data is related to a given explanatory variable.

A t-value map can be produced for each explanatory variable of interest. Each map indicates how strongly vertices on the surface are related to one explanatory variable. Parameter estimates can also be compared to see if one explanatory variable is more strongly related to the data than another. To encode this kind of hypothesis, one PE is subtracted from another using a "contrast" vector such as [1 -1 0 0 ...], a combined standard error is computed, and a new t-map is generated. In a similar fashion, to test for a more complicated collection of effects, a matrix of contrast weights can be specified. A more rigorous description of single and multiple linear regression and GLM, types of analyses, estimation and hypothesis testing is available at [http://www.statsoft.com/textbook/stglm.html http://www.statsoft.com/textbook/stglm.html].

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.

3.1 Using the FreeSurfer Group Descriptor file to create a design matrix BR The FreeSurfer Group Descriptor File (FSGDF) Format (Version 1) 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.

3.1.1 Formatting the FSGDF BR The FSGDF format uses tags to identify information, as shown below:

  Example of a legal file:
  ------------------------- cut here ------------------
  GroupDescriptorFile 1
  Title MyTitle
  Class Class1 plus blue
  CLASS Class2 circle green
  SomeTag
  Variables             Age  Weight   IQ
  Input subjid1 Class1   10    100   1000
  Input subjid2 Class2   20    200   2000
  #Input subjid3 Class2   20    200   2000
  DefaultVariable Age
  ------------------------- cut here ------------------

Notes:

• The first line of the file must be "GroupDescriptorFile 1".

• Title is not necessary. This will be used for display.

• Class only needs the class name, the next two items, if present, will be used in the display.

• The third input will be treated as a comment, due to the # at the beginning of the line.

• DefaultVariable is the default variable for display.

• SomeTag is not a valid tag, so it will be ignored.

General rules:

• Tags are NOT case sensitive (for instance, 'Class' and 'CLASS' are the same tag).
• Labels are case sensitive.
• When multiple items appear on a line, they can be separated by any white space (i.e., blank or tab).
• Any line where # appears as the first non-white space character is treated as a comment (ignored).

• The Variables line should appear before the first Input line.

• All Class lines should appear before the first Input line.
• Variable label replications are not allowed.
• Class label replications are not allowed.
• Subject Id replications are not allowed.
• If a class label is not used, a warning is printed out.

• DefaultVariable must be a member of the Variable list.

• No error is generated if a tag does not match.
• Empty lines are OK.
• A class label can optionally be followed by a class marker.

• A class marker can optionally be followed by a class color.BR

First exercise:

• [wiki:FsTutorial/CreateFsgdFile Exercise A. Create an FSGD file]

3.1.2 Creating a Design Matrix BR

The FSGDF is specified in the command-line for mris_glm 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:

• doss (different offset, same slope): this will create a design matrix in which each class has its own offset but forces all classes to have the same slope.

• dods (different offset, different slope): this value models each class with its own offset and slope.

• none: this value is used if neither of the previous models work for your particular analysis. Using this value requires that you specify the design matrix manually.

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

Note: It is not necessary to run mris_glm 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 mris_glm (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> \
  --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.

• [wiki:FsTutorial/CreateContrastVectors Exercise B. Specify contrast vectors to test 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> \
  --glmdir ?h.fsgd_file.glmdir \
  --pca \
  --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. The flag --pca is used to perform PCA/SVD analysis on the residual.

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  pca-eres/       rvar.mgh  Xg.dat

the outputs are as follows:

• ar1.mgh - BR beta.mgh - all regression coefficientsBR eres.mgh - residual errorBR fsgd.X.mat - BR mri_glmfit.log - execution parametersBR rstd.mgh - residual error stddev (just sqrt of rvar)BR rvar.mgh - residual error varianceBR Xg.dat - BR y.fsgd - fsgd fileBR

The two subdirectories created, contrast/ (will be named after your contrast file) and pca-eres/, contain some additional outputs. In the contrast/ directory there will be:

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

The outputs are as follows:

• C.dat - copy of contrast matrixBR F.mgh - F-ratioBR gamma.mgh - contrastBR sig.mgh - significance from F-test (actually -log10(p))BR

In the pca-eres/ directory there will be:

sdiag.mat  stats.dat  u.mat  v.mgh

The outputs are as follows:

• sdiag.mat - singular valuesBR v.mgh - spatial eigenvectorsBR u.mat - frame eigenvectorsBR

In addition, there is stats.dat with 5 columns:BR

• (1) component numberBR (2) variance spanned by that componentBR (3) cumulative variance spanned up to that componentBR (4) percent variance spanned by that componentBR (5) cumulative percent variance spanned up to that componentBR

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

• [wiki:FsTutorial/ComputeContrast Exercise C. Construct command lines for group analysis]

BR

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 (see if FDR is sufficient:).
2. Simulation: run the simulator with the same parameters
   • as the estimation to get the Cluster Simulation Data (CSD).
3. Clustering: run mri_surfcluster or mri_volcluster with 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 (or mri_volcluster).

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. Useable options are:

• (1) perm - perumation, randomly permute rows of X (cf FSL randomise) BR (2) mc-full - replace input with white gaussian noiseBR (3) mc-z - do not actually do analysis, just assume the output is z-distributed (cf ANFI AlphaSim)BR

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).

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.

A sample command, to run the permutation simulation with mri_glmfit is:

mri_glmfit --y ?h.fsgd_file.thickness.10.mgh --fsgd fsgd_file.txt --surf <average_subject> ?h --C contrast.mat --sim perm 10000 3 csd1

this will create 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 --surf <average_subject> ?h --C contrast.mat --sim perm 5000 3 csd1
mri_glmfit --y ?h.fsgd_file.thickness.10.mgh --fsgd fsgd_file.txt --surf <average_subject> ?h --C contrast.mat --sim perm 5000 3 csd2

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

If your design matrix is only one column of all ones then add the flag --perm-1 to the command line.

If running the monte carlo simulation instead, you will need to specify your smoothing level again with the flag --fwhm <your_fwhm>. The command will look like this:

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

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