Differences between revisions 95 and 96
Deletions are marked like this. Additions are marked like this.
Line 1: Line 1:
#acl AdminGroup:read,write,delete,revert AdrianKCLee:read,write All:

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

= TUTORIAL LOCATION HAS MOVED =
Good news for those of you who are confused by FSGD files, and mri_glmfit commands: There is a new group analysis tool in Freesurfer, Qdec. To accompany this tool there is also a new group analysis tutorial.

'''Please see the new group analysis tutorial, using Qdec. Click''' [[FsTutorial/QdecGroupAnalysis|here]]

= FreeSurfer Tutorial: Group Analysis =
 *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.

In this tutorial, you will learn how to perform statistical analysis of group surface-based data, including:<<BR>>
 *Making an average subject from your set of subjects
 *Constructing a FreeSurfer Group Descriptor File (FSGD)<<BR>>
 *Preprocessing the group data<<BR>>
 *Constructing the design matrix<<BR>>
 *Constructing contrast matrices to test hypotheses<<BR>>
 *Correcting for multiple comparisons<<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 volumes and surfaces of subjects in:

$FREESURFER_HOME/subjects/buckner_data/group_study

(Note: the data in the group_study directory is not required to complete the tutorial. However, [[FsTutorial/Data|it is available for download]] if you wish to pursue your own group analysis.)

== 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 [[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 (e.g., thicknesses for each subject at a vertex), '''X''' is the known design matrix (e.g., 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 [[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 [[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/stats
}}}

{{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 [[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:

 * ''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 (default).

 * ''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 mri_glmfit now to create the design matrix, as mri_glmfit will create it for you later in this exercise.'''

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

=== Preprocessing steps ===

'''These are sample commands that you should run for the group analysis. In interest of time these steps have already been run for you, and the output is found in {{{$FREESURFER_HOME/subjects/buckner_data/tutorial_subjs/group_analysis_tutorial/stats}}}'''

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/stats
}}}

now you can run the mris_preproc command:
{{{
mris_preproc --fsgd gender_age.txt \
  --target fsaverage \
  --hemi lh \ 		[[FsTutorial|Back to Tutorial Top]]

<<TableOfContents>>

= Introduction =

This tutorial is designed to introduce you to the "command-line" group
analysis stream in FreeSurfer (as opposed to QDEC which is
GUI-driven), including correction for multiple comparisons. While this
tutorial shows you how to perform a surface-based thickness study, it
is important to realize that most of the concepts learned here apply
to any group analysis in FreeSurfer, surface or volume, thickness or
fMRI.

Here are some useful Group Analysis Links:<<BR>>
[[FsTutorial/QdecGroupAnalysis|QDEC Tutorial]]<<BR>>
[[FsgdFormat | FSGD Format]]<<BR>>
[[FsgdExamples|FSGD Examples]]<<BR>>
[[DodsDoss| DODS vs DOSS]]<<BR>>
[[FsTutorial/GroupAnalysisJP|Old Group Analysis Tutorial]]<<BR>>

= This Data Set =

The data used for this tutorial is 40 subjects from Randy Buckner's
lab. It consists of males and females ages 18 to 93. You can see the
demographics [[FsTutorial/GroupAnalysisDemographics|here]]. You will
perform an analysis looking for the effect of age on cortical
thickness accounting for the effects of gender in the analysis. This
is the same data set used in the
[[FsTutorial/QdecGroupAnalysis|QDEC Tutorial]],
and you will get the same result.

= General Linear Model (GLM) DODS Setup =

== Design Matrix/FSGD File ==

For this example, we will model the thickness as a straight line. A
line has two parameters: intercept (or offset) and a slope.
  1. The slope is the change of thickness with age
  1. The intercept/offset is interpreted as the thickness at age=0.
Parameter estimates are also called "regression coefficients".

To account for effects of gender, we will model each sex with its
own line, meaning that there will be four linear parameters (also
called "betas"):
  1. Intercept for Females
  1. Intercept for Males
  1. Slope for Females
  1. Slope for Males
In FreeSurfer, this type of design is called [[DodsDoss| DODS]] (for
"Different-Offset, Different-Slope").

In FreeSurfer, you can either create your own design matrices, or, if
you can specify your design in terms of a FreeSurfer Group Descriptor
File [[FsgdFormat| (FSGD)]], FreeSurfer will create them for you. The
FSGD file is a simple text file. See [[FsgdFormat|this page]] for the
format. The [[FsTutorial/GroupAnalysisDemographics|demographics page]]
also has an example FSGD file for this data.

Things to do:
  1. Create an FSGD file for the above design. One (gender_age.fsgd)
  already exists so that you can continue with the exercises.

== Contrasts ==

A contrast is a vector that embodies the hypothesis we want to
test. In this case, we wish to test the change in thickness
with age after removing the effects of gender. To do this, create a
simple text file with the following numbers:

{{{
0 0 0.5 0.5
}}}

Notes:
  1. There is one value for each parameter (so 4 values)
  1. The intercept/offset values are 0 (nuisance)
  1. The slope values are 0.5 so as to average the F and M slopes
  1. A file called lh-Avg-thickness-age-Cor.mtx already exists

= Analysis =

== Assemble the Data (mris_preproc) ==

Assembling the data simply means:
  1. Resampling each subject's data into a common space, and
  1. Concatenating all the subjects' into a single file
  1. Spatial smoothing (can be done between 1 and 2)
This can be done in two equivalent ways:

=== Previously Cached (qcached) Data ===

When you have run recon-all with the -qcache option, recon-all will
resample data onto the average subject (fsaverage) and smooth it at
various FWHM (full-width/half-max), usually 0, 5, 10, 10, 20, and
25mm. This can speed later processing. The data for this tutorial have
been cached, so run:

{{{
mris_preproc --fsgd gender_age.fsgd \
  --cache-in thickness.fwhm10.fsaverage \
  --target fsaverage --hemi lh \
  --out lh.gender_age.thickness.10.mgh
}}}

Notes:
  1. Only takes a few seconds because the data have been cached
  1. The FSGD file lists all the subjects, helps keep order.
  1. The independent variable is the thickness smoothed to 10mm FWHM.
  1. The data are for the left hemisphere
  1. The output is lh.gender_age.thickness.10.mgh (which already exists)

Things to do:
  1. Run mri_info on lh.gender_age.thickness.10.mgh to find its
  dimensions.

=== Uncached Data ===

In the case that you have no cached the data, you can perform the same
operations manually using the two commands below. There is no need to
do this for this tutorial.

{{{
OPTIONAL: THIS WILL TAKE ABOUT 10 MINUTES
mris_preproc --fsgd gender_age.fsgd \
  --target fsaverage --hemi lh \
Line 101: Line 128:
  --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.
{{{ 		  --out lh.gender_age.thickness.00.mgh
}}}
Notes:
  1. This resamples each subjects data to fsaverage
  1. Output is lh.gender_age.thickness.00.mgh, which is unsmoothed

{{{
OPTIONAL: THIS WILL TAKE ABOUT 5 MINUTES
Line 111: Line 141:
  --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 ===

'''As stated earlier, these are sample commands that you should run for the group analysis. In interest of time these steps have already been run for you, and the output is found in {{{$FREESURFER_HOME/subjects/buckner_data/tutorial_subjs/group_analysis_tutorial/stats}}}'''

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 [[FsTutorial/CreateContrastVectors|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. [[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\ 		  --tval lh.gender_age.thickness.10B.mgh
}}}
  1. This smooths by 10mm FWHM
  1. "--cortex" means only smooth areas in cortex (exclude medial
  wall). This is automatically done with qcache. You can also specify
  other labels.
  1. Output is lh.gender_age.thickness.10B.mgh.

== GLM Analysis (mri_glmfit) ==

{{{
mri_glmfit \
  --y lh.gender_age.thickness.10.mgh \
  --fsgd gender_age.fsgd dods\
  --C lh-Avg-thickness-age-Cor.mtx \
  --surf fsaverage lh \
  --cortex \
  --glmdir lh.gender_age.glmdir
}}}
Notes:
  1. Input is lh.gender_age.thickness.10.mgh
  1. Same FSGD used as with mris_preproc. Maintains subject order!
  1. DODS is specified (it is the default)
  1. Only one contrast is used (lh-Avg-thickness-age-Cor.mtx) but you
  can specify multiple contrasts.
  1. "--cortex" specifies that the analysis only be done in cortex
  (ie, medial wall is zeroed out). Other labels can be used.
  1. The output directory is lh.gender_age.glmdir
  1. Should only take about 1min to run.

Things to do:

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:

{{{
y.fsgd -- copy of input FSGD file (text)
Xg.dat -- design matrix (text)
mask.mgh -- binary mask (surface overlay)
beta.mgh -- all parameter estimates (surface overlay)
rstd.mgh -- residual standard deviation (surface overlay)
sar1.mgh -- residual spatial AR1 (surface overlay)
fwhm.dat -- average FWHM of residual (text)
dof.dat -- degrees of freedom (text)
lh-Avg-thickness-age-Cor -- contrast subdirectory
mri_glmfit.log -- log file (text, send this with bug reports)
}}}

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 .mtx extension). The '''lh-Avg-thickness-age-Cor'''
directory will have the following files:

Study Questions:
  1. What was the DOF for this experiment?
  1. What was the FWHM?
  1. How many frames does beta have? Why?

Look at the contrast subdirectory:

{{{
C.dat -- original contrast matrix (text)
gamma.mgh -- contrast effect size (surface overlay)
F.mgh -- F ratio of contrast (surface overlay)
sig.mgh -- significance, -log10(pvalue), uncorrected (surface overlay)
}}}

View the uncorrected significance map with tksurfer:

{{{
tksurfer fsaverage lh inflated \
  -annot aparc.annot -fthresh 2 \
  -overlay lh.gender_age.glmdir/sig.mgh \
}}}

You should see:

{{attachment:groupana.lat.uncor.th2.gif}}
{{attachment:groupana.med.uncor.th2.gif}}

Notes:
  1. The threshold is set to 2, meaning vertices with p<.01,
  uncorrected, will have color.
  1. Blue means a negative correlation (thickness decreases with age),
  red is positive correlation.
  1. Click on a point in the Precentral Gyrus. What is it's value?
  What does it mean?

Viewing the medial surface, change the overlay threshold to something
very, very low (say, .01):

{{{
View->Configure->Overlay, set "Min" .01
}}}

You should see:

{{attachment:groupana.med.uncor.th0.gif}}

Notes:
  1. Almost all of cortex now has color
  1. The non-cortical areas are still blank (0) because they were
  excluded with --cortex in mri_glmfit above.

== Clusterwise Correction for Multiple Comparisons ==

To perform a cluster-wise correction for multiple comparisons, we will
run a simulation. The simulation is a way to get a measure of the
distribution of the maximum cluster size under the null
hypothesis. This is done by iterating over the following steps:
  1. Synthesize a z map
  1. Smooth z map (using residual FWHM)
  1. Threshold z map (level and sign)
  1. Find clusters in thresholded map
  1. Record area of maximum cluster
  1. Repeat over desired number of iterations (usually > 5000)

In FreeSurfer, this information is stored in a simple text file called
a CSD (Cluster Simulation Data).

Once we have the distribution of the maximum cluster size, we correct
for multiple comparisons by:
  1. Going back to the original data
  1. Thresholding using same level and sign
  1. Finding clusters in thresholded map
  1. For each cluster, p = probability of seeing a maximum cluster
  that size or larger during simulation.

=== Run the simulation ===

All these steps are performed with the mri_glmfit-sim:

{{{
mri_glmfit-sim \
Line 140: Line 278:
  --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:

 ar1.mgh - spatial AR1 coefficients <<BR>>
 beta.mgh - all regression coefficients<<BR>>
 eres.mgh - residual error<<BR>>
 fsgd.X.mat - final design matrix in matlab format <<BR>>
 mri_glmfit.log - execution parameters<<BR>>
 rstd.mgh - residual error stddev (just sqrt of rvar)<<BR>>
 rvar.mgh - residual error variance<<BR>>
 Xg.dat - design matrix in text format <<BR>>
 y.fsgd - fsgd file<<BR>>

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:

 C.dat - copy of contrast matrix<<BR>>
 F.mgh - F-ratio<<BR>>
 gamma.mgh - contrast<<BR>>
 sig.mgh - significance from F-test (actually -log10(p))<<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. Note: 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 simulator synthesizes data under the null hypothesis, analyzes, thresholds, clusters, and then keeps track of the number of times clusters of a given size were encounted. 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:
 (1) perm - perumation, randomly permute rows of X (cf FSL randomise) <<BR>>
 (2) mc-full - full Monte Carlo simulation, replace input with white gaussian noise, smooth, and analyze. <<BR>>
 (3) mc-z - z Monte Carlo simulation. Does not actually do analysis, just assume the output is z-distributed (cf ANFI AlphaSim)<<BR>>
Which one should you use? Permutation makes the fewest assumptions and is fastest, but it requires that the design matrix be orthogonal (e.g., you cannot have continuous variables such as age or IQ). If you cannot use permutation, then a mc-z will be fastest, but the z assumption requires that you have many subjects (on the order of 80). Otherwise, you must run the full MC simulation, which can take quite a while. For both MC simulations, you must supply the smoothness of your data as Full-Width/Half-Max (FWHM) of the residual. This can be obtained from the ResidualFWHM value in y.fsgd in the glm output directory.

The next parameter is ''nsim'' which corresponds to the number of simulations to run. If you want to make inferences at the .01 level then you'll need about 10000 iterations. 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 .01:

{{{
mri_glmfit --y lh.gender_age.thickness.10.mgh \
  --glmdir lh.gender_age.glmdir \
  --fsgd gender_age.txt doss \
  --surf fsaverage lh \
  --fwhm 14.517 --C age.mat \
  --sim mc-full 10000 2 lh.gender_age.glmdir/csd1
}}}

this will create lh.gender_age.glmdir/csd1-age.csd

If you want to split this into multiple runs you could use the following two commands:
{{{
mri_glmfit --y lh.gender_age.thickness.10.mgh \
  --fsgd gender_age.txt doss\
  --fwhm 14.517 \
  --glmdir lh.gender_age.glmdir \
  --surf fsaverage lh \
  --C age.mat \
  --sim mc-full 5000 2 lh.gender_age.glmdir/csd1
mri_glmfit --y lh.gender_age.thickness.10.mgh \
  --fsgd gender_age.txt doss\
  --fwhm 14.517 \
  --glmdir lh.gender_age.glmdir \
  --surf average lh \
  --C age.mat \
  --sim mc-full 5000 2 lh.gender_age.glmdir/csd2
}}}

which will generate csd1-age.csd and csd2-age.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 --in lh.gender_age.glmdir/age/sig.mgh \
  --csd lh.gender_age.glmdir/csd1-age.csd \
  --csd lh.gender_age.glmdir/csd2-age.csd \
  --sum lh.gender_age.glmdir/age/sig.cluster.sum \
  --cwsig lh.gender_age.glmdir/age/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 NVtxs
   1 1.850 134208 391.63 -28.8 34.0 -8.0 0.00100 0.00060 0.00140 542
   2 1.439 18669 0.56 -10.9 9.1 -13.2 0.99960 0.99930 0.99980 1
   3 1.366 148132 13.84 -7.4 11.8 -14.1 0.97440 0.97240 0.97640 25
}}}

CWP stands for ''cluster-wise probability'' which 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. NVtxs is the number of vertexes in the cluster. The output surface map, sig.cluster.mgh, will be a map of these clusters with their CWP. This can be viewed with: 		  --sim mc-z 5 4 mc-z.negative \
  --sim-sign neg \
  --overwrite
}}}

Notes:
  1. Specify the same GLM directory (--glmdir)
  1. The simulation type is Z Monte Carlo (mc-z)
  1. Specify the sign (neg for negative, pos, or abs)
  1. Vertex-wise thershold of 4 (p < .0001)
  1. Only 5 iterations (usually > 5000)
  1. CSD files will have mc-z.negative base
  1. --overwrite deletes previous CSD files

This has been run this data with 1000 iterations (CSD base of
mc-z.neg). This can often take many hours.

=== View CSD File ===

The CSD files are stored in lh.gender_age.glmdir/csd. Look
at one of them:

{{{
less lh.gender_age.glmdir/csd/mc-z.neg4.j001-lh-Avg-thickness-age-Cor.csd
}}}
Or click [[FsTutorial/GroupAnalysisCsd|here]].

Notes:
  1. Rows that begin with '#' are headers or comments.
  1. Each data row is an iteration.
  1. The 3rd column is the maximum cluster size for that iteration.
  1. The CSD file itself is not particularly important, so don't worry
  too much about it.

=== View the Corrected Results ===

In the contrast subdirectory, you will see several new files, all
beginning with mc-z:
{{{
mc-z.neg4.sig.cluster.summary - summary of clusters (text)
mc-z.neg4.sig.cluster.mgh - cluster-wise corrected map (overlay)
mc-z.neg4.sig.ocn.annot - annotation of clusters
}}}

First, look at the cluster summary (or click
[[FsTutorial/GroupAnalysisClusterSummary|here]]):
{{{
less lh.gender_age.glmdir/lh-Avg-thickness-age-Cor.csd/mc-z.neg4.sig.cluster.summary
}}}
Notes:
  1. This is a list of all the clusters that were found (38 of them)
  1. The CWP column is the cluster-wise probability (the number you
  are interested in). It is a simple p (ie, NOT -log10(p)).
  1. Note that ALL clusters, regardless of their significance.
  1. Cluster number 1 has a CWP of p=.001.
  1. Cluster number 21 has a CWP of p=.007.

Load the cluster annotation in tksurfer
Line 278: Line 339:
  -overlay lh.gender_age.glmdir/age/sig.cluster.mgh
}}} 		  -annot lh.gender_age.glmdir/lh-Avg-thickness-age-Cor.csd/mc-z.neg4.sig.ocn.annot \
  -fthresh 2
}}}

You should see:

{{attachment:groupana.lat.clusters.gif}}
{{attachment:groupana.med.clusters.gif}}

Notes:
  1. These are all clusters, regardless of significance.
  1. When you click on a cluster, the label will tell you the
  cluster number (eg, cluster-021).

Change the label mode to "Outline" by hitting the outline button
{{attachment:tksurfer.outline.gif}}. Now load the cluster p-value
overlay with
{{{
  File->LoadOverlay
  "Browse..."
  Select mc-z.neg4.sig.cluster.mgh
  OK
  OK
}}}

You should see:

{{attachment:groupana.lat.cor.th2.gif}}
{{attachment:groupana.med.cor.th2.gif}}

Note:


Things to do:

  1. Find and click on cluster 1. It is blue and has a value of -3
  since this is -log10(.001). The -3 is because the correlation is
  negative.
  1. Find and click on cluster 21. Its value is -1.11 because this
  is -log10(.007). Note that it is transparent because its significance
  is worse than the threshold we set (-fthresh 2, p < .01).
  1. All vertices within a cluster are the same value (the p-value
  of the cluster).
  1. You can changed the cluster-wise threshold by clicking on
  View->Configure->Overlay, and setting the "Min" to your desired
  level. As you do this, some clusters will change transparency.

Back to Tutorial Top

Contents

  1. Introduction
  2. This Data Set
  3. General Linear Model (GLM) DODS Setup
    1. Design Matrix/FSGD File
    2. Contrasts
  4. Analysis
    1. Assemble the Data (mris_preproc)
      1. Previously Cached (qcached) Data
      2. Uncached Data
    2. GLM Analysis (mri_glmfit)
    3. Clusterwise Correction for Multiple Comparisons
      1. Run the simulation
      2. View CSD File
      3. View the Corrected Results

Introduction

This tutorial is designed to introduce you to the "command-line" group analysis stream in FreeSurfer (as opposed to QDEC which is GUI-driven), including correction for multiple comparisons. While this tutorial shows you how to perform a surface-based thickness study, it is important to realize that most of the concepts learned here apply to any group analysis in FreeSurfer, surface or volume, thickness or fMRI.

Here are some useful Group Analysis Links:<<BR>> QDEC Tutorial
FSGD Format
FSGD Examples
DODS vs DOSS
Old Group Analysis Tutorial

This Data Set

The data used for this tutorial is 40 subjects from Randy Buckner's lab. It consists of males and females ages 18 to 93. You can see the demographics here. You will perform an analysis looking for the effect of age on cortical thickness accounting for the effects of gender in the analysis. This is the same data set used in the QDEC Tutorial, and you will get the same result.

General Linear Model (GLM) DODS Setup

Design Matrix/FSGD File

For this example, we will model the thickness as a straight line. A line has two parameters: intercept (or offset) and a slope.

  1. The slope is the change of thickness with age
  2. The intercept/offset is interpreted as the thickness at age=0.

Parameter estimates are also called "regression coefficients".

To account for effects of gender, we will model each sex with its own line, meaning that there will be four linear parameters (also called "betas"):

  1. Intercept for Females
  2. Intercept for Males
  3. Slope for Females
  4. Slope for Males

In FreeSurfer, this type of design is called DODS (for "Different-Offset, Different-Slope").

In FreeSurfer, you can either create your own design matrices, or, if you can specify your design in terms of a FreeSurfer Group Descriptor File (FSGD), FreeSurfer will create them for you. The FSGD file is a simple text file. See this page for the format. The demographics page also has an example FSGD file for this data.

Things to do:

  1. Create an FSGD file for the above design. One (gender_age.fsgd) already exists so that you can continue with the exercises.

Contrasts

A contrast is a vector that embodies the hypothesis we want to test. In this case, we wish to test the change in thickness with age after removing the effects of gender. To do this, create a simple text file with the following numbers: 

0 0 0.5 0.5

Notes:

  1. There is one value for each parameter (so 4 values)
  2. The intercept/offset values are 0 (nuisance)
  3. The slope values are 0.5 so as to average the F and M slopes
  4. A file called lh-Avg-thickness-age-Cor.mtx already exists

Analysis

Assemble the Data (mris_preproc)

Assembling the data simply means:

  1. Resampling each subject's data into a common space, and
  2. Concatenating all the subjects' into a single file
  3. Spatial smoothing (can be done between 1 and 2)

This can be done in two equivalent ways:

Previously Cached (qcached) Data

When you have run recon-all with the -qcache option, recon-all will resample data onto the average subject (fsaverage) and smooth it at various FWHM (full-width/half-max), usually 0, 5, 10, 10, 20, and 25mm. This can speed later processing. The data for this tutorial have been cached, so run: 

mris_preproc --fsgd gender_age.fsgd \
  --cache-in thickness.fwhm10.fsaverage \
  --target fsaverage --hemi lh \
  --out lh.gender_age.thickness.10.mgh

Notes:

  1. Only takes a few seconds because the data have been cached
  2. The FSGD file lists all the subjects, helps keep order.
  3. The independent variable is the thickness smoothed to 10mm FWHM.
  4. The data are for the left hemisphere
  5. The output is lh.gender_age.thickness.10.mgh (which already exists)

Things to do:

  1. Run mri_info on lh.gender_age.thickness.10.mgh to find its dimensions.

Uncached Data

In the case that you have no cached the data, you can perform the same operations manually using the two commands below. There is no need to do this for this tutorial. 

OPTIONAL: THIS WILL TAKE ABOUT 10 MINUTES
mris_preproc --fsgd gender_age.fsgd \
  --target fsaverage --hemi lh \
  --meas thickness \
  --out lh.gender_age.thickness.00.mgh

Notes:

  1. This resamples each subjects data to fsaverage
  2. Output is lh.gender_age.thickness.00.mgh, which is unsmoothed 

OPTIONAL: THIS WILL TAKE ABOUT 5 MINUTES
mri_surf2surf --hemi lh \
  --s fsaverage \
  --sval lh.gender_age.thickness.mgh \
  --fwhm 10 \
  --cortex \
  --tval lh.gender_age.thickness.10B.mgh
  1. This smooths by 10mm FWHM
  2. "--cortex" means only smooth areas in cortex (exclude medial wall). This is automatically done with qcache. You can also specify other labels.
  3. Output is lh.gender_age.thickness.10B.mgh.

GLM Analysis (mri_glmfit)

mri_glmfit \
  --y lh.gender_age.thickness.10.mgh \
  --fsgd gender_age.fsgd dods\
  --C lh-Avg-thickness-age-Cor.mtx \
  --surf fsaverage lh \
  --cortex \
  --glmdir lh.gender_age.glmdir

Notes:

  1. Input is lh.gender_age.thickness.10.mgh
  2. Same FSGD used as with mris_preproc. Maintains subject order!
  3. DODS is specified (it is the default)
  4. Only one contrast is used (lh-Avg-thickness-age-Cor.mtx) but you can specify multiple contrasts.
  5. "--cortex" specifies that the analysis only be done in cortex (ie, medial wall is zeroed out). Other labels can be used.
  6. The output directory is lh.gender_age.glmdir
  7. Should only take about 1min to run.

Things to do:

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: 

y.fsgd -- copy of input FSGD file  (text)
Xg.dat -- design matrix  (text)
mask.mgh -- binary mask (surface overlay)
beta.mgh -- all parameter estimates (surface overlay)
rstd.mgh -- residual standard deviation (surface overlay)
sar1.mgh -- residual spatial AR1  (surface overlay)
fwhm.dat -- average FWHM of residual (text)
dof.dat  -- degrees of freedom (text)
lh-Avg-thickness-age-Cor -- contrast subdirectory
mri_glmfit.log -- log file (text, send this with bug reports)

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 .mtx extension). The lh-Avg-thickness-age-Cor directory will have the following files:

Study Questions:

  1. What was the DOF for this experiment?
  2. What was the FWHM?
  3. How many frames does beta have? Why?

Look at the contrast subdirectory: 

C.dat -- original contrast matrix (text)
gamma.mgh -- contrast effect size (surface overlay)
F.mgh -- F ratio of contrast  (surface overlay)
sig.mgh -- significance, -log10(pvalue), uncorrected (surface overlay)

View the uncorrected significance map with tksurfer: 

tksurfer fsaverage lh inflated \
  -annot aparc.annot -fthresh 2 \
  -overlay lh.gender_age.glmdir/sig.mgh \

You should see:

Notes:

  1. The threshold is set to 2, meaning vertices with p<.01, uncorrected, will have color.

  2. Blue means a negative correlation (thickness decreases with age), red is positive correlation.
  3. Click on a point in the Precentral Gyrus. What is it's value? What does it mean?

Viewing the medial surface, change the overlay threshold to something very, very low (say, .01): 

View->Configure->Overlay, set "Min" .01

You should see:

Notes:

  1. Almost all of cortex now has color
  2. The non-cortical areas are still blank (0) because they were excluded with --cortex in mri_glmfit above.

Clusterwise Correction for Multiple Comparisons

To perform a cluster-wise correction for multiple comparisons, we will run a simulation. The simulation is a way to get a measure of the distribution of the maximum cluster size under the null hypothesis. This is done by iterating over the following steps:

  1. Synthesize a z map
  2. Smooth z map (using residual FWHM)
  3. Threshold z map (level and sign)
  4. Find clusters in thresholded map
  5. Record area of maximum cluster

  6. Repeat over desired number of iterations (usually > 5000)

In FreeSurfer, this information is stored in a simple text file called a CSD (Cluster Simulation Data).

Once we have the distribution of the maximum cluster size, we correct for multiple comparisons by:

  1. Going back to the original data
  2. Thresholding using same level and sign
  3. Finding clusters in thresholded map
  4. For each cluster, p = probability of seeing a maximum cluster that size or larger during simulation.

Run the simulation

All these steps are performed with the mri_glmfit-sim: 

mri_glmfit-sim \
  --glmdir lh.gender_age.glmdir \
  --sim mc-z 5 4 mc-z.negative \
  --sim-sign neg \
  --overwrite

Notes:

  1. Specify the same GLM directory (--glmdir)
  2. The simulation type is Z Monte Carlo (mc-z)
  3. Specify the sign (neg for negative, pos, or abs)

  4. Vertex-wise thershold of 4 (p < .0001)

  5. Only 5 iterations (usually > 5000)

  6. CSD files will have mc-z.negative base
  7. --overwrite deletes previous CSD files

This has been run this data with 1000 iterations (CSD base of mc-z.neg). This can often take many hours.

View CSD File

The CSD files are stored in lh.gender_age.glmdir/csd. Look at one of them: 

less lh.gender_age.glmdir/csd/mc-z.neg4.j001-lh-Avg-thickness-age-Cor.csd

Or click here.

Notes:

  1. Rows that begin with '#' are headers or comments.
  2. Each data row is an iteration.
  3. The 3rd column is the maximum cluster size for that iteration.
  4. The CSD file itself is not particularly important, so don't worry too much about it.

View the Corrected Results

In the contrast subdirectory, you will see several new files, all beginning with mc-z: 

mc-z.neg4.sig.cluster.summary - summary of clusters (text)
mc-z.neg4.sig.cluster.mgh - cluster-wise corrected map (overlay)
mc-z.neg4.sig.ocn.annot - annotation of clusters

First, look at the cluster summary (or click here): 

less lh.gender_age.glmdir/lh-Avg-thickness-age-Cor.csd/mc-z.neg4.sig.cluster.summary

Notes:

  1. This is a list of all the clusters that were found (38 of them)
  2. The CWP column is the cluster-wise probability (the number you are interested in). It is a simple p (ie, NOT -log10(p)).
  3. Note that ALL clusters, regardless of their significance.
  4. Cluster number 1 has a CWP of p=.001.
  5. Cluster number 21 has a CWP of p=.007.

Load the cluster annotation in tksurfer 

tksurfer fsaverage lh inflated \
  -annot  lh.gender_age.glmdir/lh-Avg-thickness-age-Cor.csd/mc-z.neg4.sig.ocn.annot \
  -fthresh 2

You should see:

Notes:

  1. These are all clusters, regardless of significance.
  2. When you click on a cluster, the label will tell you the cluster number (eg, cluster-021).

Change the label mode to "Outline" by hitting the outline button tksurfer.outline.gif. Now load the cluster p-value overlay with 

  File->LoadOverlay
  "Browse..."
  Select  mc-z.neg4.sig.cluster.mgh 
  OK
  OK

You should see:

Note:

Things to do:

  1. Find and click on cluster 1. It is blue and has a value of -3 since this is -log10(.001). The -3 is because the correlation is negative.
  2. Find and click on cluster 21. Its value is -1.11 because this is -log10(.007). Note that it is transparent because its significance

    is worse than the threshold we set (-fthresh 2, p < .01).

  3. All vertices within a cluster are the same value (the p-value of the cluster).
  4. You can changed the cluster-wise threshold by clicking on

    View->Configure->Overlay, and setting the "Min" to your desired level. As you do this, some clusters will change transparency.

FsTutorial/GroupAnalysis_tktools (last edited 2014-01-28 16:08:46 by LouisVinke)