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)
2) Task information. Onset times, durations, and names for each event type (condition) of interest
3) Nuisance covariates. Any custom nuisance covariates we want to add, including
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)
ans = '/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'
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
SPM_onsets_and_nuisance_regressors.mat c1sub-sid001567_acq-MPRAGE_T1w.nii c2sub-sid001567_acq-MPRAGE_T1w.nii c3sub-sid001567_acq-MPRAGE_T1w.nii c4sub-sid001567_acq-MPRAGE_T1w.nii c5sub-sid001567_acq-MPRAGE_T1w.nii canlab_robust_reg_model1 mwc1sub-sid001567_acq-MPRAGE_T1w.nii mwc2sub-sid001567_acq-MPRAGE_T1w.nii mwc3sub-sid001567_acq-MPRAGE_T1w.nii rp_sub-sid001567_task-pinel_acq-s1p2_run-03_bold.txt rsub-sid001567_task-pinel_acq-s1p2_run-03_bold.mat rsub-sid001567_task-pinel_acq-s1p2_run-03_bold.nii rsub-sid001567_task-pinel_acq-s1p2_run-03_bold.nii.gz rsub-sid001567_task-pinel_acq-s1p2_run-03_bold_reorient.mat spm_model1 sub-sid001567_acq-MPRAGE_T1w.nii sub-sid001567_acq-MPRAGE_T1w.nii.gz sub-sid001567_acq-MPRAGE_T1w_reorient.mat sub-sid001567_acq-MPRAGE_T1w_seg8.mat sub-sid001567_task-pinel_acq-s1p2_run-03_bold.json sub-sid001567_task-pinel_acq-s1p2_run-03_bold.mat sub-sid001567_task-pinel_acq-s1p2_run-03_bold.nii sub-sid001567_task-pinel_acq-s1p2_run-03_bold.nii.gz sub-sid001567_task-pinel_acq-s1p2_run-03_events.tsv swrsub-sid001567_task-pinel_acq-s1p2_run-03_bold.nii wrsub-sid001567_task-pinel_acq-s1p2_run-03_bold.nii wrsub-sid001567_task-pinel_acq-s1p2_run-03_bold.nii.gz y_sub-sid001567_acq-MPRAGE_T1w.nii

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);
Using default mask: /Users/torwager/Documents/GitHub/CanlabCore/CanlabCore/canlab_canonical_brains/Canonical_brains_surfaces/brainmask_canlab.nii loading mask. mapping volumes. checking that dimensions and voxel sizes of volumes are the same. Pre-allocating data array. Needed: 145720800 bytes Loading image number: 150 Elapsed time is 2.552240 seconds. Image names entered, but fullpath attribute is empty. Getting path info. .fullpath should have image name for each image column in .dat Attempting to expand image filenames in case image list is unexpanded 4-D images Number of unique values in dataset: 493 Bit rate: 8.95 bits Warning: Number of unique values in dataset is low, indicating possible restriction of bit rate. For comparison, Int16 has 65,536 unique values
% Plot the images. This will output potential outliers too, but this has been pre-prepped here.
plot(dat);
Calculating mahalanobis distances to identify extreme-valued images ...based on union of corr...Retained 2 components for mahalanobis distance Expected 50% of points within 50% normal ellipsoid, found 22.67% Expected 7.50 outside 95% ellipsoid, found 6 Potential outliers based on mahalanobis distance: Bonferroni corrected: 1 images Cases 1 Uncorrected: 6 images Cases 1 2 3 4 5 6 ...and cov...Retained 2 components for mahalanobis distance Expected 50% of points within 50% normal ellipsoid, found 18.67% Expected 7.50 outside 95% ellipsoid, found 6 Potential outliers based on mahalanobis distance: Bonferroni corrected: 1 images Cases 1 Uncorrected: 6 images Cases 1 2 3 4 5 6 done. Outliers: Outliers after p-value correction: Image numbers: 1 Image numbers, uncorrected: 1 2 3 4 5 6
Grouping contiguous voxels: 1 regions
Grouping contiguous voxels: 1 regions
Grouping contiguous voxels: 1 regions

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:
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);