%% CANlab Robust Regression Walkthrough
% This script demonstrates how to use the CANlab robust regression toolbox to run robust 2nd-level analyses.
% Robust regression is a popular approach for principled down-weighting of outliers in a General Linear Model analysis. Robust regression has particular advantages when conducting maps over many tests, as in many neuroimaging analyses. In such cases, outliers cannot be easily checked and dealt with for each individual test.  The application of robust regression to neuroimaging is described in this paper:
% Wager, T. D., Keller, M. C., Lacey, S. C. & Jonides, J. (2005). Increased sensitivity in neuroimaging analyses using robust regression. Neuroimage. 26:99-113. PMID: 15862210. [PDF]
% Table of Contents
% Theory
% Installation and code
%     Toolbox installation
%     Main functions in the robust regression toolbox
% Sample datasets
% Robust regression with robfit.m
%     Covariates
%     robfit inputs
%     Load the dataset
%     Load the behavioral data
%     Define the  structure with input variables for robfit
%     Create a directory for the results and go to it
% Run robfit
%     Examine the output files
%     Get batch results
% robfit_parcelwise
%% Theory
% • Outliers can violate assumptions, have very large effects on regression coefficients
% • High-variance observations tend to dominate if their observed values are extreme
% • When assumptions cannot be checked at each voxel, automatic procedures for weighting based on outlier status advantageous
% • Robust regression: An automatic procedure for identifying cases that are potential outliers and down-weighting
% • Iterative GLS, estimation procedure, but weights are based on residual value rather than variancePotential way to deal with heteroscedastic variances
%
%
%% Installation and code
%% Toolbox installation
% The robust regression toolbox is available from Github, here. For code, dependencies, and installation, along with more information on these objects and CANlab tools, see canlab.github.io and the CANlab Core Tools github repository.
% In brief, you will need SPM software (for image reading/writing only) and the Robust toolbox on your Matlab path. For full functionality, you'll also need  the CANlab Core Tools repository, including subfolders, on your path.
% Download or "clone" (add to your Github destkop) these two toolboxes and add each to your Matlab path with subfolders:
% https://github.com/canlab/CanlabCore
% https://github.com/canlab/RobustToolbox
% In addition, you will need Statistical Parametric Mapping (SPM, SPM2/5/12) installed and on your Matlab path. This is used for image reading and writing, but not for statistics.
% If you are using robfit_parcelwise (see below), you will also need:
% https://github.com/canlab/Neuroimaging_pattern_masks
% Installation. To install these, place the source code folders in a folder on your hard drive, and type
% >> pathtool
% at the Matlab command prompt.
% Select each of the folders named above, and select “add with subfolders”
% SPM and Matlab signal processing toolboxes required
% The robust regression toolbox is not an SPM toolbox per se.  It uses SPM image manipulation (data I/O) functions, but does not rely on any SPM functions for statistics.  It does, however, use the Matlab Statistics Toolbox, so you will need that as well.
%% Main functions in the robust regression toolbox
% The main functions used in the toolbox are:
% robfit.m
% Runs robust regression for a series of group analyses, creating a directory and saving results for each
% robfit_parcelwise.m
% Runs robust regression for a single group analysis on each of a series of "parcels", regions of interest that are larger than voxels and generally defined to conform to pre-defined anatomical or functional boundaries.
% robustfit.m
% Matlab's internal function, which is run for each test conducted across brain voxels or parcels.
% robseed.m
% Robust regression correlation analysis, correlating each voxel with a seed region you specify.
% robust_results_batch.m
% Threshold maps and generate and save plots and tables of results.
% publish_robust_regression_report.m
% Thresholds maps and generates an HTML report with results
% Note: If you are using object-oriented tools in CANlab Core, the regress( ) method for fmri_data objects also has a robust regression option. This works independent of the robust regression toolbox.
%% Sample datasets
% A sample dataset is  in the "Sample_datasets" folder in CANlab_Core_Tools.
% See canlab.github.io for more sample datasets, and in particular the function load_image_set( ) for some easy, pre-packaged datasets, atlases, and maps.
% Some others are in this repository as well: https://github.com/canlab/CANlab_data_public
% This walkthough will use emotion regulation data in the "Sample_datasets" folder, in this subfolder: "Wager_et_al_2008_Neuron_EmotionReg"
% The dataset is a series of contrast images from N = 30 participants. Each image is a contrast image for [reappraise neg vs. look neg] for one participant.
% These data were published in:
% Wager, T. D., Davidson, M. L., Hughes, B. L., Lindquist, M. A., Ochsner, K. N.. (2008). Prefrontal-subcortical pathways mediating successful emotion regulation. Neuron, 59, 1037-50.
% We will load it with load_image_set( ).
%% Robust regression with robfit.m
% Robfit runs a series of analyses.  Each analysis is a univariate robust regression analysis, regressing data from each voxel in a series of images on a GLM design matrix.
% The typical use case for a single analysis is to enter a series of contrast images, one for each of n participants.
% robfit can handle a series of contrast images, with one set of n images per contrast.  It creates directories called robust0001, robust0002, etc., one per analysis (i.e., per contrast).
% A SETUP.mat file for each analysis contains information about the input images and design matrix, along with some other meta-data.
%% Covariates
% You can run robfit to test the group mean (intercept) only or you can enter one or more covariates.  Each covariate is entered as a predictor in the design matrix.
% As with any multiple regression, the intercept will be tested (with associated t- and p-values) at the zero-point of all covariates, i.e., when they are zero.  Therefore, if you mean-center covariates and use effects codes (1, -1) for categorical variables, the intercept will be interpreted as the group effect for the average participant, assuming input images correspond to participants.
% Robfit will save t-maps and P-maps for each regressor in your design matrix, with the intercept first (rob_tmap_0001.nii) followed by covariate effects if covariates are entered (rob_tmap_0002.nii, etc.)
% Robfit automatically mean-centers continuous covariates, but does not mean-center categorical covariates, and instead codes them with effects codes of 1 and -1.
%
%
%% robfit inputs
% robfit uses a data structure with a specific format that stores information about image files and the second-level experimental design.  This structure has a special variable name here, by convention: EXPT.
% The main 2nd-level analysis function is called robfit.m.  Once you create an EXPT structure, you will pass it into robfit.m.
% Here is a list of the fields in EXPT required to ultimately run analyses with robfit:
% EXPT.subjects       List of subject directory names, in a cell array.{‘subj1’ ‘subj2’ … etc.}
% EXPT.SNPM.P          Cell array; each cell contains names of a set of images (one per subject) to be subjected to a group analysis. Each image list is a string matrix of image file names. Thus, images for multiple 2nd level analyses can be specified at once by entering images for each analysis into separate cells in EXPT.SNPM.P
% EXPT.SNPM.connames  String matrix of names, whose rows correspond to names for the analyses.  Each row (name) corresponds to the cell with the same index in EXPT.SNPM.P.
%                      E.g.,
%                      ['Faces-Houses'
%                     'Houses-Faces']
% EXPT.SNPM.connums   A vector of integers numbering the analyses, and corresponding to cells in EXPT.SNPM.P.  These numbers determine output directory names, which will be called:
%                      robust_0001
%                      robust_0002
%                      …and so forth (one directory per 2nd-level analysis)
% EXPT.mask            Name of mask image containing in-analysis voxels. This need not be in the same space as the images you’re analyzing.
%
% EXPT.cov      Matrix of subjects x covariates for analysis of between subjects effects. One row per subject, one column per predictor. Do not include an intercept here.  It will be added to the front end of the design matrix, as the first predictor.
% The output images, i.e., rob_tmap_0001.img, are numbered so that 0001 is the intercept (overall activation) and 0002 – 0000x are the maps for covariates 1 thru (x-1)
%% Load the dataset

% This script uses object-oriented tools, which we cover in another tutorial. 
% For now, it's just a convenient way to get a set of image filenames.

[data_obj, subject_names, image_names] = load_image_set('emotionreg');

% data_obj is an fmri_data object, with its own methods and properties
% We'll use the names

%
%% Load the behavioral data


beh_data_file = which('Wager_2008_emotionreg_behavioral_data.txt')
if isempty(beh_data_file), error('Sample behavioral data not found. '); end
beh = readtable(beh_data_file);

%% Define the  structure with input variables for robfit
% Here we fill in fields of a structure, called EXPT, that contain what we need to run the analysis. We'll do this for only a sngle contrast, though we could add more cells to the fields to add others as well.

% The fields of EXPT should be:
%
% EXPT.SNPM.P           Cell vector. Each cell specifies images for one analysis. Each cell contains a string matrix with image names for the analysis. Image files can be 3-D or 4-D images.
% EXPT.SNPM.connames    Cell vector. Each cell contains a string with the analysis name (e.g., contrast name) for this contrast
% EXPT.SNPM.connums     Vector of contrast numbers. Determines folder names (e.g., 1 = robust0001 
% EXPT.cov              [n x k] matrix of n observations (must match number of images) empty for 1-sample ttest, or containing covariates
% EXPT.mask             Optional mask file image name, for voxels to include

EXPT.subjects = subject_names;
EXPT.SNPM.P{1} = image_names{1};  % string matrix of images names
% Note: In the sample data file, 30 images for 30 subjects are saved in 1
% 4-D .nii file. So we only need to list the name of this one file. 
% If you had 3-D image files, this would be a string matrix with 30 rows

EXPT.SNPM.connums = 1;
EXPT.SNPM.connames = 'Reapp_vs_Look';
EXPT.mask = which('gray_matter_mask.img');
EXPT.cov = beh.Y_Reappraisal_Success;   

% Because we have one column in EXPT.cov, our 2nd-level design matrix 
% will contain two images (e.g., two t-maps/p-maps): 
% One for the intercept (always first) and one for the covariate entered.

% Let's look at the input structure we have created:

EXPT

%
%% Create a directory for the results and go to it
% This directory will contain subdirectories called robust00XX for each cell containing images entered in EXPT.SNPM.P.
% Make sure you are starting in a directory that you can create a folder for the sample results in before you run this.

basedir = fullfile(pwd, 'Robust_regression_sample_results_dir');

if ~exist(basedir, 'dir'), mkdir(basedir), end
cd(basedir)

% Save the EXPT setup file for later reference:
save EXPT EXPT

%
%% Run robfit

EXPT = robfit(EXPT);

% Try help robfit for more options.

% e.g., you can specify that you want to run only some of the contrasts 
% in EXPT.SNPM.P.  You can specify that you want to create OLS images 
% side by side with robust images to compare.  
% And you can specify a mask filename.  
% The mask image should contain 1’s and 0’s, with 1’s in voxels you want to analyze.

%
% Note: The robfit command will also attempt to run publish_robust_regression_report for each set of images (each analysis) it runs. This generates an HTML report in the results directory and opens it. Thus, if this command runs correctly, it will bring up a report (or multiple reports).
%
%% Examine the output files
%

ls -lt

% published_output             A folder with HTML reports from publish_robust_regression_report
% irls-ols_p_0002.nii           P-value image for the covariate, difference between robust and OLS
% irls-ols_z_0002.nii           Z-value image for the covariate, difference between robust and OLS
% ols_p_0002.nii                 P-value image for the covariate, OLS regression
% ols_tmap_0002.nii           t-value image for the covariate, OLS regression
% ols_beta_0002.nii            beta (slope) image for the covariate, OLS regression
% rob_p_0002.nii                P-value image for the covariate, robust regression
% rob_tmap_0002.nii         t-value image for the covariate, robust regression
% rob_beta_0002.nii          beta (slope) image for the covariate, robust regression
% irls-ols_p_0001.nii          P-value image for the intercept, difference between robust and OLS
% irls-ols_z_0001.nii          The images below are the same as those above, but for the intercept
% ols_p_0001.nii
% ols_tmap_0001.nii
% ols_beta_0001.nii
% rob_p_0001.nii
% rob_tmap_0001.nii
% rob_beta_0001.nii
% weights.nii                    Robust regression weights for each image (i.e., participant)
% mask.nii                        mask of voxels included in the analysis
% nsubjects.nii                Image with integer values for number of subjects with valid data
% SETUP.mat                 Metadata file
% If you are using CANlab object-oriented tools to threshold and view the results maps, you can load and combine the t- and P-maps into t-statistic image objects for each regressor.
% [trob, names, mask_obj, nsubjects, weights, SETUP] = robust_reg_load_files_to_objects(pwd);
%
%
%% Get batch results
% Results should have been saved by publish_robust_regression_report as
% HTML. To re-run batch results, go to the robust regression folder and
% run the batch script.
% You can edit this script to customize it as well.

cd(basedir)
cd('robust0001')
robust_results_batch

% Or,. run publish_robust_regression_report( ) to create a report and save an HTML file in the robust0001 folder.
% publish_robust_regression_report
% See the object oriented help for more details on how to run robust analyses, and how to extract and plot data from significant regions.
%
%% robfit_parcelwise
% Voxel-wise results are often under-powered with moderate samples, and it may be advantageous to test averages across a series of parcels instead.  If the parcel boundaries match the underlying boundaries of activated regions, this could substantially increase power.
% robfit_parcelwise extracts mean values for each of a series of regions -- with a default CANlab atlas parcellation that covers the brain -- and runs robust regressions on each parcel
% Corrections (i.e., False Discovery Rate) are performed across parcels rather than voxels
% For convenience in viewing results, parcel-by-parcel results are returned in voxel-wise maps where all values in the parcel have the same statistic values.
% robfit_parcelwise uses the object-oriented tools, so you need images in an fmri_data object to run it.  We have created one above, in the variable data_obj.  See canlab.github.io for more details.
%

OUT = robfit_parcelwise(data_obj);

%  When you run robfit_parcelwise, it generates additional plots of the mask, data, and weights by default.
% It also saves a structure with t-statitic objects, results tables, region-class objects with extracted data, and more:

OUT

% Parcel-by-parcel betas (slopes), t-values, P-values, image counts:
% betas: [489×1 double]
% tscores: [489×1 double]
% pvalues: [489×1 double]
% nsubjects: [489×1 double]
% maskvol: [489×1 double]
% Parcel-by-parcel robust regression weights, degrees of freedom, and data  (parcel averages)
% weights: [489×30 double]
% dfe: [489×1 double]
% datmatrix: [30×489 double]
% False discovery rate threshold and significant parcels
% pthr_FDRq05: 0.0156
% pthr_FDRq05_descrip: 'Highest P-value significant at FDR q < 0.05; threshold at p <= pthr'
% sig_q05: [489×1 logical]
% cohens_d_fdr05: 0.4206
% cohens_d_fdr05_descrip: 'Min Cohen's d detectable at FDR q < 0.05'
% Quality control metrics and outliers for input images:
% group_metrics: [1×1 struct]
% individual_metrics: [1×1 struct]
% global_gm_wm_csf_values: [30×3 double]
% outliers_corr: [30×1 logical]
% outliers_uncorr: [30×1 logical]
%  ind_quality_dat: [30×9 table]
% t-maps, slopes, mask, and subject counts in voxel-wise fmri_data object format:
% t_obj: [1×1 statistic_image]
% beta_obj: [1×1 fmri_data]
% mask: [1×1 fmri_data]
% nsubjects_obj: [1×1 fmri_data]
% resultstable: [1×8 table]
% Region-class objects with extracted data:
% region_objects: {[1×10 region]}
% Table objects with significant effects (FDR):
% contrast_tables_FDR05: {[10×8 table]}
%
% CANlab object-oriented tools, or other toolboxes, can be used to visualize the results.
% Let's view the results and do some post-run interactive visualization:

o2 = canlab_results_fmridisplay(OUT.t_obj, 'montagetype', 'full', 'noverbose');

% o2 = removeblobs(o2);
% o2 = addblobs(o2, OUT.t_obj);


disp(OUT.contrast_tables_FDR05{1})

% Extract data from all significant regions to make a plot of individual subject post-hoc values.  Note that this induces a voxel-selection bias ("circularity") that inflates the apparent effect sizes in these regions. But it can be useful for visualizing individual subjects' data and distributions. You can also use the same code for regions of interest defined a priori, which are free of this bias.
% Region-class object stored in OUT.region_objects{1} is saved in the "r" variable for convenience.
% The extract_data( ) method extracts the data from the original contrast images.
% Regions are already auto-labeled with a standard atlas, and names stored in the .shorttitle field.

r = extract_data(OUT.region_objects{1}, data_obj);
dat_matrix = cat(2, r.dat);
create_figure('individual subjects');
barplot_columns(dat_matrix, 'names', {r.shorttitle});

%
%
%
%   Adjust for global intensity and global components and re-run:
%   imgs2 = imgs.rescale('l2norm_images');
%   OUT = robfit_parcelwise(imgs2);
%
%   Run adding global WM and CSF signal covariates:
%   OUT = robfit_parcelwise(imgs2, 'csf_wm_covs');
%
%   run with global covariates and removing outliers based on Mahalanobis distance:
%   OUT = robfit_parcelwise(imgs2, 'csf_wm_covs', 'remove_outliers');