NeuroElf wiki

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== MDM::ComputeGLM - compute a random-effects suitable GLM (regression) ======

===== Motivation =====
To combine several subjects' data into one compound analysis, it is custom to compute random-effects statistics on the second level. The MDM file format of BrainVoyager QX references the required time-course files (VTC or MTC) and the respective information containing the model (onsets in [[xff - PRT format|PRT]] files or regressors in [[xff - SDM format|SDM]] files).

This method then allows to perform the regression on a session-by-session basis and combine the results in one compound GLM object.

===== Requirements =====
This function is currently implemented for BrainVoyager QX's GLM format only. The following files must be available:
  * VTC (or MTC) files for all runs of all subjects (e.g. the ones automatically created by spm5_preprojobs, for how to name files see below)
  * [[xff - PRT format|PRT]] (or the derived [[xff - SDM format|SDM]]) files for all runs (if subjects share this information, only one file is needed for each uniquely specified model!)
  * optionally the realignment parameter files (supported formats: TXT for SPM rp*.txt files, [[xff - SDM format|SDM]] for BV's motion correction output)

**Note: The combination of within-subjects regressors (e.g. same condition across several runs) is done by identifying the subject ID as the particle of the time-course filename before the first underscore character!** This means that VTC (or MTC) filenames **must be named** as follows:

  * SUBJ1_RUNorTASK1.vtc
  * SUBJ1_RUNorTASK2.vtc
  * SUBJ1_ ...
  * SUBJ2_RUNorTASK1.vtc
  * SUBJ2_RUNorTASK2.vtc
  * ...

The number of runs per subject can be different (although, for the sake of having the error on beta estimates coming from the same distribution, it is suggested to use roughly the same number of volumes--and usually runs--for all subjects), and also the names **after** the first underscore need not match across subjects (in fact, it needn't be unique, which of course is highly recommended for the user to be able to distinguish time-course files!). The following is thus a valid naming scheme:

<file>jd3512_r48217.vtc
jd3512_r48218.vtc
jd3512_r48219.vtc
pk5199_r89112.vtc
pk5199_r89113.vtc</file>

===== Reference / mdm.Help('ComputeGLM') =====
<file> MDM::ComputeGLM  - compute a GLM from an MDM file
 
 FORMAT:       glm = mdm.ComputeGLM([options])
 
 Input fields:
 
       options     optional 1x1 struct with fields
        .globsigs  add global signals as confound, one of
                   0 - none
                   1 - entire dataset (above threshold/within mask)
                   2 - two (one per hemisphere, split at BV Z=128)
                   3 or more, perform PCA of time courses and first N
                   xff object(s), extract average time course from masks
        .ithresh   intensity threshold, default: 100
        .loadglm   boolean flag, load GLM file named in .outfile
        .mask      optional masking, default: no mask (for now only VTC)
        .motpars   motion parameters (Sx1 cell array with sdm/txt files)
        .motparsd  also add diff of motion parameters (default: false)
        .motparsq  also add squared motion parameters (default: false)
        .ndcreg    if set > 0, perform deconvolution (only with PRTs!)
        .orthconf  orthogonalize confounds (and motion parameters, true)
        .outfile   output filename of GLM file, default: no saving
        .ppicond   list of regressors (or differences) to interact
        .ppirob    perform robust regression on VOI timecourse and remove
                   outliers from timecourse/model (threshold, default: 0)
        .ppivoi    VOI object used to extract time-course from
        .ppivoiidx intra-VOI-object index (default: 1)
        .prtpnorm  normalize parameters of PRT.Conds (true)
        .redo      selected subjects will be overwritten (default: false)
        .regdiff   flag, regress first discreet derivatives (diff) instead
        .restcond  remove rest condition (rest cond. name, default: '')
        .robust    perform robust instead of OLS regression
        .savesdms  token, if not empty, save on-the-fly SDMs (e.g. '.sdm')
        .showsdms  token, passed to SDM::ShowDesign (if valid)
        .shuflab   PRT labels (conditions names) to shuffle
        .shuflabm  minimum number of onsets per label (1x1 or 1xL)
        .sngtrial  single-trial GLM (only with PRTs, default: false)
        .sngtskip  condition list to skip during single-trial conversion
        .subsel    cell array with subject IDs to work on
        .tfilter   add filter regressors to SDMs (cut-off in secs)
        .tfilttype temporal filter type (either 'dct' or {'fourier'})
        .transio   boolean flag, if true, save GLM and use transio
        .vweight   combine runs/studies variance-weighted (default: false)
        .xconfound just as motpars, but without restriction on number
 
 Output fields:
 
       glm         GLM object
 
 Note: if RFX flag in MDM is set to true, predictor separation will be
       set to "Subjects".
       if .outfile is given, GLM is saved after each subject for robust
       regression models (to allow later continuation of crashed/broken
       computation using the .subsel field)
       for .ppi to work, the model filenames must be PRTs!
       all additional fields for the call to PRT::CreateSDM are supported!</file>

===== Options =====

==== .globsigs ====
Setting this field with either
  * a number (other than 0) will add as many "global signals" to each design matrix before regression, whereas for 1 this will be the overall mean, for 2 the mean represents left and right hemisphere separately, and for 3 or higher a PCA is performed and as many components as requested added to each design matrix
  * a single VMR object or a cell array of VMR objects will add as many "global signals" (from all non-0 VMR voxels, sampled at the VTC coordinates) as requested

==== .ithresh ====
This field controls any explicit masking based on threshold (mean timecourse value). If left empty/unset, no threshold masking occurs (full datasets).

==== .loadglm ====
If set to ''true'' (and a valid GLM filename is provided in ''.outfile''!), an already existing GLM file is loaded and then additional subjects added (or replaced) as configured.

==== .mask ====
Either a [[xff - MSK format|MSK object]] fitting the VTC spatial layout or a binary mask sufficing the same condition can be supplied to compute the GLM only within those voxels that are non zero/not false in the mask.

==== .motpars ====
If set to a cell array of filenames with size ''mdm.NrOfStudies''-by-1, the regression will contain any motion parameters found in the referenced files. If set to ''true'' any files found under ''mdm.RunTimeVars.MotionParameters'' will be used instead.

==== .motparsd ====
If set to true any motion parameters will be added as is and their first derivative will also be part of the model.

==== .motparsq ====
If set to true any motion parameters will be added as is and their squared (and mean-removed) version will also be part of the model.

==== .ndcreg ====
If set to an integer number greater than 0, [[xff - PRT format|PRTs]] will be converted to design matrices using deconvolution designs (stick functions for each deconvolution/volume regressor; field is passed to [[prt.CreateSDM|PRT::CreateSDM]] under the same fieldname).

==== .orthconf ====
If set to true, confounds (motion parameters, filters, global timecourses, etc.) will be orthogonalized against one another (which doesn't change the outcome of the regression equation as far as the betas of interest is concerned, but can improve the behavior of the underlying functions used, when the design matrix is inverted before computing the betas).

==== .outfile ====
Can be set to a filename under which the GLM will be stored at the end of the computation (and, for robust regression and enabled transio after each subject). If ''.loadglm'' is also set to true and file already exists, the GLM will be loaded first, and only selected subjects from the MDM that are not already computed will be added at the end (unless ''.redo'' is set to true, in which case those subjects' data, if it is already part of the GLM, will be replaced). See also ''.loadglm'' and ''.redo'' options.

==== .ppicond ====
Can be set to a cell array with condition and/or contrast combinations, which will be used by adding an interaction term between an extracted timecourse (from a PPI seed region) and a (contrast of) regressor(s). To specify a contrast, condition names need to be put into a string like this:
  * '''condition_A + condition_B''' (sum of two conditions)
  * '''condition_A > condition_B''' (difference of two conditions)
  * '''condition_A + condition_B > condition_C + condition_D''' (difference of sums, can be used for a condition interaction term)
See also ''.ppivoi'' and ''.ppivoiidx'' options. 

==== .ppirob ====
**This option is experimental and not fully tested!**

Can be used to remove outliers of the PPI timecourse (which then improves RFX properties of beta estimates associated with the PPI regressors).

==== .ppivoi ====
Must be set to a valid [[xff - VOI format|VOI object]] for PPI to be used.

==== .ppivoiidx ====
Can be set to an index to select a specific VOI from the VOI object/file in option ''.ppivoi'' (default is 1). Timecourses will be extracted from this selected VOI.

==== .prtpnorm ====
If set to true, normalize and parameters in PRT files before computing the parametric conditional regressors (is passed as ''.pnorm'' into the [[prt.CreateSDM|PRT::CreateSDM]] method).

==== .redo ====
If set to ''true'' and ''.loadglm'' is also set to ''true'' **and** an existing GLM file is provided in the ''.outfile'' option, any selected subjects (see ''.subsel'' option) will be overwritten in the GLM file.

==== .regdiff ====
**This option is experimental and not fully tested!**

If set to ''true'' design matrices and data are first passed through ''diff'' before regression occurs (in effect regression the first derivative of data on the first derivative of the design).

==== .restcond ====
This option is passed as ''.rcond'' into [[prt.CreateSDM|PRT::CreateSDM]] calls and defaults to empty (don't remove any condition as rest).

==== .robust ====
If set to ''true'' performs robust regression on all design matrices (**very slow!**) using the [[fitrobustbisquare_img]] function.

==== .savesdms ====
Can be used to store created/altered [[xff - SDM format|SDM files]], and must be set to a non-empty string containing '''.sdm''' at its end.

==== .showsdms ====
If set to either '''image''', '''ortho''', or '''plot''' calls [[sdm.ShowDesign|SDM::ShowDesign]] before each run's regression starts (visualizing the design matrix of the currently regressed run in the desired way).

==== .shuflab ====
**This option is experimental and not fully tested!**

Allows to re-generate the GLM under the null hypothesis of equal effect of stimuli across conditions by shuffling the labels (onsets across conditions).

==== .shuflabm ====
Sub-option for ''.shuflab'', ensures that each condition with onsets retains at least this number of onsets.

==== .sngtrial ====
If set to ''true'' (and the MDM object references PRT objects), calls [[prt.ConvertToSingleTrial|PRT::ConvertToSingleTrial]] prior to converting protocols into design matrices. Uses the ''.sngtskip'' option.

==== .sngtskip ====
Instructs the [[prt.ConvertToSingleTrial|PRT::ConvertToSingleTrial]] method which conditions are to be left unchanged (e.g. onsets of cues, fixation points, active baseline stimuli, etc.)

==== .subsel ====
Cell array of subject IDs, sub-selects subjects from those referenced by the MDM file (see [[#requirements|Requirements section]] above for more info).

==== .tfilter ====
Sets the temporal filtering cutoff in seconds, so that filtering regressors will be added to the design matrix.

==== .tfilttype ====
Sets the temporal filtering type, one of '''dct''' (for Discrete Cosine Transform) or '''fourier''' (Fourier basis set).

==== .transio ====
If set to ''true'' and ''.outfile'' is set to a valid GLM file, will use [[@transio]] access to store data into the output file. This is necessary in case the output file is very large (bigger than 2GB, for instance), so that no out-of-memory error occurs.

==== .vweight ====
If set to ''true'' beta weights of identically named conditions across runs of the same subject will not simply be combined as a sum but will be weighted by the square root of the variance estimate of the design (reciprocal value of the diagonal values of the inverted design matrix), which is roughly equivalent to computing a single regression over all of one subject's runs.

==== .xconfound ====
Can be set to a cell array with filenames of size ''mdm.NrOfStudies''-by-1, in which case additional regressors of no interest (confounds) are read from those files and added to the design matrix (e.g. heart-rate data, eye-tracking output, etc.). The files' content must match the number of volumes of the timecourse files!

===== Usage examples =====
Starting from scratch, the user has to do the following steps:
  * locate the filenames of the VTC (or MTC and SSM/TSM) files to be regressed
  * locate the filenames of the [[xff - PRT format|PRT]] (or [[xff - SDM format|SDM]]) files to be used as regressors (model)
  * optionally locate the filenames of the realignment parameter files
  * make settings in the MDM structure
  * pass optional settings to the ComputeGLM method call

<code matlab mdm_computeglm_sample.m>% locate filenames...
cd /Projects/MAN1/Imaging/Subjects
vtcs = findfiles([pwd '/M*/func'], '*.vtc', 'depth=1');
prts = findfiles([pwd '/M*/model'], '*.prt', 'depth=1');
if numel(vtcs) ~= numel(prts)
    error('Number of VTCs and PRTs must match!');
end

% locate realignment parameter files
rps = findfiles([pwd '/M*/func/r*'], 'rp*.txt', 'depth=1');

% create MDM
mdm = xff('new:mdm');
mdm.RFX_GLM = 1;
mdm.PSCTransformation = 1;
mdm.NrOfStudies = numel(vtcs);
mdm.XTC_RTC = [vtcs(:), prts(:)];

% optionally save MDM, .SaveAs without argument opens a dialog!
% mdm.SaveAs;

% ComputeGLM options
opts = struct( ...
    'motpars',   {rps(:)}, ...
    'restcond',  'rest', ...
    'robust',    true, ...
    'tfilter',   160, ...
    'tfilttype', 'fourier');

% compute GLM (robust takes a **LONG TIME** ! e.g. for each 200-Volume
% run with 15 regressors between 30 and 60 minutes, depending on the
% hardware used)
glm = mdm.ComputeGLM(opts);

% save GLM
glm.SaveAs;</code>