Dear experts:
Perhaps I am just too bleary-eyed to see the error of my ways, but for whatever reason, I am having problems adding condition derivatives to a first-level glm.
My script is a little complicated. I am iterating through two runs of data, constructing their respective design matrices, concatenating, then adding extra regressors to account for between run variance. My glm constructs, stats maps look reasonable, but I believe the default ‘glover’ hrf is used, rather than the model I’ve specified, as the derivative regressors are missing from the design matrix.
Thanks in advance for your input, the code, built on version 0.10.3, is below! The relevant block is towards the very end of the script.
Best,
JB
import os
import nibabel as nb
import numpy as np
import pandas as pd
import configparser
import matplotlib.pyplot as plt # Import matplotlib
from nilearn import datasets, surface, plotting, image, maskers
from nilearn.image import concat_imgs, mean_img
from nilearn.plotting import view_img, plot_design_matrix, plot_contrast_matrix
from nilearn.reporting import make_glm_report
from nilearn.glm.first_level import make_first_level_design_matrix, FirstLevelModel
from nilearn.maskers import NiftiMasker
from nilearn.image import concat_imgs, mean_img
# Function to read and parse the config file
def read_config(config_file):
config = configparser.ConfigParser()
config.read(config_file)
return config
# Main script
config_file = 'config.ini'
config = read_config(config_file)
# Config variables
study = config['Paths']['study']
codedir = config['Paths']['codedir']
sublist = config['Paths']['sublist']
prepdir = config['Paths']['prepdir']
outdir = config['Paths']['outdir']
space = config['Variables']['space']
mask = config['Variables']['mask']
session = config['Variables']['session']
task = config['Variables']['task']
tr = config.getfloat('Variables', 'tr')
#hrf_model = config['Variables']['hrf_model'].split('+') # Split components
confounds = config['Variables']['confounds'].split(',')
run_mapping = {'1': 'AP', '2': 'PA'}
with open(sublist, "r") as subjects_list:
subjects = subjects_list.read().splitlines()
for s in subjects:
sub = f'sub-{s}'
sub_outdir = os.path.join(outdir, sub)
print(f'Starting 1st level GLM for {sub}')
design_matrices = []
fmri_images = []
scans = []
mask_path = os.path.join(prepdir, sub, session, 'func', f'{sub}_{mask}')
print(mask_path)
for run_identifier in run_mapping:
run_name = run_mapping[run_identifier]
confounds_path = os.path.join(prepdir, sub, session, 'func', f'{sub}_{session}_task-{task}_dir-{run_name}_desc-confounds_timeseries.tsv')
events_path = os.path.join(outdir, sub, f'{sub}_{session}_task-{task}_dir-{run_name}_events.tsv')
fmri_path = os.path.join(prepdir, sub, session, 'func', f'{sub}_{session}_task-{task}_dir-{run_name}_space-{space}_desc-preproc_bold.nii.gz')
if not all(os.path.exists(path) for path in [fmri_path, mask_path, events_path, confounds_path]):
print(f"Warning: Missing one or more files for {run_name} at {sub}")
continue
# Load events, confounds
confounds_df = pd.read_table(confounds_path, usecols=confounds).fillna(0)
events = pd.read_table(events_path).fillna(0)
# Load the fMRI data directly as 4D Nifti image
fmri_data = nb.load(fmri_path)
fmri_images.append(fmri_data)
# Build design matrix
num_scans = fmri_data.shape[-1]
frame_times = np.arange(num_scans) * tr
scans.append(num_scans) # Append number of scans to the scans list
design_matrix = make_first_level_design_matrix(frame_times, events, drift_model=None, add_regs=confounds_df)
design_matrices.append(design_matrix)
# Concatenate all runs' fmri data
if fmri_images:
concatenated_fmri = concat_imgs(fmri_images, auto_resample=True) # Adjust auto_resample as needed
# Concatenate all design matrices
if design_matrices:
design_matrix = pd.concat(design_matrices, ignore_index=True)
# Total scans
n_scans = sum(scans)
print(n_scans)
# Create binary regressors for each run
run_regressors = np.zeros((n_scans, len(scans)))
current_scan_index = 0
for idx, scan_count in enumerate(scans):
run_regressors[current_scan_index:current_scan_index + scan_count, idx] = 1
current_scan_index += scan_count
# Add run regressors to the design matrix
for idx in range(len(scans)):
design_matrix[f'Run_{idx+1}'] = run_regressors[:, idx]
# Specify contrasts
contrast_matrix = np.eye(design_matrix.shape[1])
basic_contrasts = {
column: contrast_matrix[i]
for i, column in enumerate(design_matrix.columns)
}
# Specify contrasts of interest
contrasts = {
"Inhibit-Go": basic_contrasts["Stop Success"] - basic_contrasts["Go Correct"],
"Fail-Go": basic_contrasts["Stop Fail"] - basic_contrasts["Go Correct"],
"Inh-Fail": basic_contrasts["Stop Success"] - basic_contrasts["Stop Fail"],
"Fail-Inh": basic_contrasts["Stop Fail"] - basic_contrasts["Stop Success"],
"effects_of_interest": np.vstack(
(basic_contrasts["Go Correct"], basic_contrasts["Stop Success"], basic_contrasts["Stop Fail"])
),
}
# Fit GLM
print("Fitting a GLM")
fmri_glm = FirstLevelModel(t_r=tr, hrf_model='spm + derivative', smoothing_fwhm=6)
fmri_glm = fmri_glm.fit(concatenated_fmri, design_matrices=design_matrix)
# Inspect design matrix
new_matrix = fmri_glm.design_matrices[0]
print(new_matrix) # Print or plot the design matrix to inspect it
# Plot design matrix
plot_design_matrix(new_matrix)
plt.title(f"Design Matrix for {sub}")
plt.show()