%% First-level model with CANlab regress( )
% The CANlab fmri_data object has a regress( ) method that will run a basic GLM. It returns a data structure with a series of statistic_image objects containing regression parameter estimates ("beta images") and t maps with P values included. These are statistically thresholded for display by default, but can be re-thresholded after running regress( ).
% This lab walks you through the steps.
%
%% Setting up a first-level analysis: Ingredients
% We're setting up for the first-level analysis.  We need five things:
% 1) Brain data. The functional fMRI data, preprocessed and ready for analysis (usually, realigned, warped to MNI space, and smoothed)
% Here, we want "swr*" images, which are Smoothed, Warped, and Realigned
% Other scaling and "denoising" (artifact-removal) steps are sometimes done in preprocessing as well
% 2) Task information. Onset times, durations, and names for each event type (condition) of interest
% In BIDS standard format, these are in a .tsv file
% We'll create cell arrays called ons, durations, and names (by convention, as used in SPM software)
% 3) Nuisance covariates. Any custom nuisance covariates we want to add, including
% "Spike" (indicator) regressors for outliers
% Movement parameters and their transformations
% Any other custom regressors, e.g., global CSF or aCompCorr regressors
% We want to be careful about how we choose these, to avoid multicolinearity with task regressors
% We will save these in one matrix, R (by convention, as used in SPM software)
% For BIDS/fmriprep, these are in *_desc-confounds_regressors.tsv
% 4) Filter. An intended HP filter cutoff in sec
%     - SPM will add the high-pass filter when it does the analysis
% 5) Contrast weights for effects (comparisons across conditions) we care about
%
%% Find and examine the brain data
%% Navigate to files
% It's helpful to start in a folder with data in it, so you don't have to navigate extensively through the SPM file selection browser.

pwd  % where am i now? (in terms of working directory, at least)

% Use cd( ) to go to your sample subject directory. Drag and drop, or copy and paste the path from Terminal.
% cd('/Users/torwager/Dropbox (Dartmouth College)/COURSES/Courses_Dartmouth/2021_3_Spring_fMRI_Class/Labs_and_Lectures_for_Instructors/Labs_and_assignments/Lab_5/Pinel_data_sample_subject_INSTRUCTOR')

ls  % list files here

%% Find the functional run and examine the data
% It's always important to look at your data!  Check: Do the brains look like brains? is the range of data values what you expected? Are there crazy outliers? Does the image space match the template space? Are the images relatively homogenous, or do some appear to be really different, i.e., on the same scale?
% We'll use CANlab object-oriented tools, specifically the fmri_data object, to load the preprocessed data and make a summary plot. Here, we want "swr*" images, which are Smoothed, Warped, and Realigned...and hopefully thus ready for analysis.

preprocessed_image_name = 'swrsub-sid001567_task-pinel_acq-s1p2_run-03_bold.nii';

% Load the images into an fmri_data object
dat = fmri_data(preprocessed_image_name);

% Plot the images. This will output potential outliers too, but this has been pre-prepped here.
plot(dat);

%
%% Load setup file with other variables we need
% The file SPM_onsets_and_nuisance_regressors should have everything we need to build the design in SPM!  It's always a good idea to save prepared variables before they are entered into SPM, etc.
% This .mat file was created by the first_level_spm12 lab, so go back to that if needed. It has the key variables:
% onsets             Cell array of onset times for each condition (10 conditions)
% durations        Cell array of event durations for each condition
% names             Cell array of names for each condition
% R                      Matrix of nuisance regressors
% C                      Structure with contrast weights and names

load SPM_onsets_and_nuisance_regressors

%
%% Construct design matrix
% Now, we've extracted the information in a format we can use to construct a design matrix. Let's try it.
% Note that we need to know the TR for this! It's 2 sec. In BIDS formatted data, this can be found in the ".json sidecar" file, in the "RepetitionTime" field. We also need the scan length in sec, which in this case is 300 sec or 5 min.

TR = 2;
X = onsets2fmridesign(onsets, TR, 300);

% let's plot it
plotDesign(onsets, [], TR);

% This returns the design matrix, X.
% We can look at the Variance Inflation Factors to see whether there are mulitcolinearity issues or not. We'll use the 'plot' option and save the handle to this figure so we can re-activate it and add things later:

fig_han = create_figure('vifs'); % Get handle to current figure

vifs = getvif(X, 0, 'plot');

%% Add nuisance covariates and check VIFs again
% Finally, let's re-calculate the VIFs for task regressors and add them to our original plot:

X3 = [X R];

figure(fig_han)   % re-activate the VIF figure

vifs3 = getvif(X3);
han = plot(vifs3(1:10), 'bo', 'MarkerSize', 8, 'MarkerFaceColor', [.5 .5 1]);

%% Attach design matrix to dat object
% The CANlab regress( ) method

X3 = intercept(X3, 'remove');
dat.X = X3;

%% Pad C with zeros so that we have an entry in C for each regressor

C.weights = [C.weights; zeros(size(X3, 2) - size(C.weights, 1), size(C.weights, 2))];

% add one for the intercept, which will be added automatically in regress().
C.weights(end+1, :) = 0;

%
%
%% Run regression
% We will add some optional inputs to save names in the output structure.

regression_results_ols = regress(dat, 'variable_names', names, ...
    'analysis_name', 'Pinel localizer 1st-level GLM', ...
    'C', C.weights, 'contrast_names', C.names);

%% Notes on a few key points:
% Output structure (regression_results_ols) contains beta and t images in statistic_image objects
% If contrasts are entered, output structure contains contrast images and contrast t images in statistic_image objects
% T images are thresholded, but beta/contrast images are not. These can be re-thresholded with the threshold( ) method
% statistic_image objects can be visuallized with any CANlab display methods, e.g., surface( ), montage( ), and orthviews( )
% Use the get_wh_image( ) method to select specific beta, t, contrast, contrast t images from statistic_image objects
% Enter condition names and contrast names as cell vectors of strings Labels are saved in the .image_labels property in output objects
% Some diagnostics, including VIFs and leverages, are automatically checked and returned
% There is an option for robust regression
% For first-level models (only), the grandmeanscale option is recommended to increase homogeneity in scale across participants
%% Show some results for selected individual conditions

% Thresholded t map for the first condition
t = get_wh_image(regression_results_ols.t, 1);

create_figure('montage'); axis off
display_obj = montage(t);

figure; surface(t);

% Montages for thresholded t images, condition [2 4 6 8]

t = get_wh_image(regression_results_ols.t, [2 4 6 8]);

create_figure('montage2'); axis off
display_obj = montage(t);

%% Show  results for contrasts

t = get_wh_image(regression_results_ols.con_t, 1:3);

create_figure('montage3'); axis off
display_obj = montage(t);

t = get_wh_image(regression_results_ols.con_t, 4:6);

create_figure('montage4'); axis off
display_obj = montage(t);

% Create and print a table for the first contrast
% at p < 0.001 and k = 50 contiguous voxels
t = get_wh_image(regression_results_ols.con_t, 1);
t = threshold(t, .001, 'unc', 'k', 50);
table(t);

%% All the CANlab display tools are available
% Here is a diagram showing some examples:
%
%% Save the results structure for later use
% Now we can save the results. You should do this in the Matlab Command Window, for your particular path, as it depends on where your local files are stored. Be careful, because save may overwrite existing files!
% cd('/Users/torwager/Dropbox (Dartmouth College)/COURSES/Courses_Dartmouth/2021_3_Spring_fMRI_Class/Labs_and_Lectures_for_Instructors/Labs_and_assignments/Lab_5/Pinel_data_sample_subject_INSTRUCTOR/')
% save canlab_1st_level_ols_regression_results regression_results_ols
%%
%
%% Run robust regression
% This runs robust regression at each voxel. It will take a while!

regression_results_robust = regress(dat, 'variable_names', names, ...
    'analysis_name', 'Pinel localizer 1st-level robust GLM', 'robust', ...
    'C', C.weights, 'contrast_names', C.names);

% Note: Here it may tell us that we have design multicolinearity, but we should consider whether this is because of colinearity among nuisance regressors only (which is fine, as long as we don't want to interpret their individual coefficients), or whether it's related to multicolinearity with the regressors of interest. Only the latter is a problem.
%% Show  results for contrasts

t = get_wh_image(regression_results_robust.con_t, 1:3);

create_figure('montage3'); axis off
display_obj = montage(t);

t = get_wh_image(regression_results_robust.con_t, 4:6);

create_figure('montage4'); axis off
display_obj = montage(t);

%% Save the results structure for later use
% Now we can save the results. You should do this in the Matlab Command Window, for your particular path, as it depends on where your local files are stored. Be careful, because save may overwrite existing files!
% cd('/Users/torwager/Dropbox (Dartmouth College)/COURSES/Courses_Dartmouth/2021_3_Spring_fMRI_Class/Labs_and_Lectures_for_Instructors/Labs_and_assignments/Lab_5/Pinel_data_sample_subject_INSTRUCTOR/')
% save canlab_1st_level_robust_regression_results regression_results_robust
%%
%
%% Key ideas: From 1st-level to 2nd-level (group) analysis
%  if you're running 1st-level models from multiple subjects and are going to use the output for group analysis, here is a workflow, and some things to keep in mind:
% 1. Create a script that loops through all subjects and runs the models
% Use the 'grandmeanscale' input keyword to standardize the overall scale across subjects
% Enter regressor and contrast names so that you have a record of what is what for later
% Save the results structure from each subject.
% When you're ready to run it, publish this script to create an HTML report with all the first-level model output for your records, so you can review it later.
% 2. Create a script that loops through all subjects and aggregates the contrast images
% You'll typically analyze contrast images for each contrast separately.
% For this, you'll need to gather images for contrast 1 from all subjects and aggregate them into an fmri_data object, so that there's one image per subject and all images reflect the same contrast (hypothesis). Ditto for contrast 2, 3, etc.
% Create a second script that loops through the subjects again and loads each results structure. Use the get_wh_image( ) method to select a contrast of interest and group all the subjects.
% Potential aggregation issues: The list of voxels and their brain coordinates must match up exactly across subjects. If some subjects are in different brain spaces / not coregistered to a common space (e.g., MNI sampled to standard voxels based on a common brain mask), you'll need to resample them. The resample_space( ) method can help. This should be rare / not an issue if you use fmri_data( ) with the default options, which saves data in a standard brain mask. However, if you have incomplete coverage of voxels for some subjects, you might end up with  some mismatches. Empty voxels are squeezed out using the remove_empty( ) method to save space. Before resampling, try replace_empty( ) on each subjects' image object to see if it's just an issue with empty voxels. The easiest way to do this is with the image_math( ) method, using 'cat' for concatenate. If obj1 is the first subject and obj2 the second:
% objcat = image_math(obj1, obj2, 'cat');
% 3. Create a script that loops through contrasts and performs 2nd-level (group) analysis
% As before, save the output for later, and publish to create a report. Date-stamping your reports and saving old ones (if you change things) is a good idea.