Robust regression toolbox

The robust regression toolbox is a tool for analyzing group activation data (mainly). It replaces the “2nd level analysis” in SPM. It can model both group results and group differences and/or brain-behavior correlations in the same model.

There is a more extensive help document in the robust_toolbox folder in matlab_code. The main function is robfit.

The main functions for getting and visualizing results are:

mediation_brain_results
robust_results_threshold
robust_results3

For more details on results see: visualizing_results

Robust nonparametric toolbox

A second set of functions runs robust regression embedded within a permutation testing framework. While this can take a long time to run, the advantage is that it can provide family-wise error rate correction for multiple comparisons based on both height and extent, and within regions of interest (ROIs) or a search mask. A major motivation for creating this separate set of functions was to be able to provide this kind of correction.

The main function is robust_reg_nonparam.m Though it is not as straightforward to use as robfit, a Powerpoint file called Robust_nonparam_flowchart.ppt, in the robust regression directory in the “nonparam” subfolder. provides some notes on the functions and what they do. Additional help is contained in the help files of each function.

One important note is that a mask can include several separate ROIs, which will be tested separately by robust_reg_nonparam as long as their voxels are not contiguous. Correction across all ROIs (map-wise) will also be done automatically.

A second important note is that there is one substantive difference between this and robfit, apart from the correction. The nonparam function pools weights across voxels within a contiguous ROI, and thus regularizes them. Because it pools the weights globally, funny weights by chance will affect the null statistics less. But it also means that this function may work best for reasonably small ROIs (though there may be a number of them) to avoid pooling weights over large areas.

An example of its usage is below.

Example setup and usage for robfit

Below is a complete example of setting up an analysis for a series of contrast images from first-level (individual subject) SPM analyses. The inputs are sets of contrast images from SPM, one image per subject.

%% Get names of 1st level SPM contrast images
 
load('/Volumes/Cort/Imaging/CORT/Analysis-NBACK/Cort1523/SPM.mat')
load('/Volumes/Cort/Imaging/CORT/Analysis-NBACK/EXPT.mat')
 
ncons = length(SPM.xCon);
for i = 1:ncons
imgs = filenames(['*/con_000' num2str(i) '.img'], 'char', 'absolute')
EXPT.SNPM.P{i} = imgs; EXPT.SNPM.connames{i} = SPM.xCon(i).name; EXPT.SNPM.connums(i) = i;
end
 
EXPT.SNPM.connames = char(EXPT.SNPM.connames{:});
 
%% save them in EXPT.SNPM
 
cd('/Volumes/Cort/Imaging/CORT/Analysis-NBACK')
save EXPT -append EXPT
 
%% Here we would add the 2nd-level covariate scores for
% each contrast in EXPT.cov.  They may be different for
% different contrasts, depending on what makes sense.
% A covariate that should be present in all the contrasts in 
% this dataset is the group ID
 
EXPT.cov = EXPT.group; EXPT.covname = 'Trier vs. Control';
 
%% Now let's run robust regression for each contrast image
 
mkdir group_robust_triervscontrol_cov
cd group_robust_triervscontrol_cov/
 
% define the mask
mask = which('brainmask.nii')
 
% run the analyses
EXPT = robfit(EXPT, 1:length(EXPT.SNPM.connums), 0, mask);

Example setup and usage for robust_reg_nonparam

The code below sets up and runs a mask-wise FWE correction, and SVC correction for any separate, contiguous regions within the mask. A mask file is required, and the mask should be resliced to the same space as the data images! See scn_map_image for how to reslice.

One important note is that a mask can include several separate ROIs, which will be tested separately by robust_reg_nonparam as long as their voxels are not contiguous. Correction across all ROIs (map-wise) will also be done automatically.

This code runs the permutations and saves a .mat file with the data structure R (for results)

SETUP.X is a design matrix with 5 columns (in this example). 5000 is the number of iterations. mask is a variable containing a mask image name SETUP.files is a list of image filenames 'names' refer to the names for each column of X. Each of these columns will be tested as a separate effect.

R = robust_reg_nonparam(SETUP.X, 5000,'mask', mask, 'tail', 'onetail', ...
      'data', SETUP.files, 'file', 'nparam_painmatrix_limbic_onetail', ...
  'names',{'Intercept' 'Placebo' 'Order' 'Study' 'Placebo_x_Study'});
This code summarizes the results in tables and figures

You can either load the R structure from the file that is automatically saved on disk, or you can use the existing one from the workspace. The first part creates the table and updates R with new SVC information. The second part creates the figures and returns contiguous clusters for other use (visualization, etc.)

R = robust_nonparam_results(R);
 
% In this example, we are showing results for regressor #2, the first one after the intercept, named 'Placebo'
[A,R,cl_extent,dat_extent,cl_extent_pos,cl_extent_neg] = robust_nonparam_displayregions(R, 2 , 't','svc');
The code below sets up a mask and runs the function

There are a few other things that are required of masks for this function that can trip you up. This code deals with two common ones: The mask must be in the space of the image data, and it must not include voxels without actual data. Also, the mask should have 1/0 values.

cd('/Volumes/Macintosh HD2/Combined_sample_2007/robust0002')
load SETUP
 
cd ../masks
 
disp(mask)
 
% reslice the mask to the data
% I found it was flipped L/R (-2 was x voxel size for the mask, but +2 for the data image)
scn_map_image(mask, SETUP.files(1, :), 'write', 'resliced_mask.img');
 
% in addition, you must exclude voxels from the mask that do not have valid
% values in any image.  the code below does this.
% It also forces the mask to have binary (1 or 0) values.
mask_image('resliced_mask.img', SETUP.files(1, :), 'resliced_mask.img', 'binary');
 
 
R = robust_reg_nonparam(SETUP.X, 5000,'mask', 'resliced_mask.img', 'tail', 'onetail', ...
      'data', SETUP.files, 'file', 'nparam_sACC_SMA_HHvLL', ...
  'names',{'Intercept' 'Placebo' 'Order' 'Study' 'Placebo_x_Study'});
 
% This would show results for the placebo regressor
[A,R,cl_extent,dat_extent,cl_extent_pos,cl_extent_neg] = robust_nonparam_displayregions(R, 2 , 't','svc');
help/robust_regression/robust_regression.txt · Last modified: 2016/08/19 18:00 (external edit)
CC Attribution-Noncommercial-Share Alike 4.0 International
Driven by DokuWiki Recent changes RSS feed Valid CSS Valid XHTML 1.0