How to correctly model and define constrasts in SPM12

1

Dear community,

I am analyzing some fMRI data acquired while participants underwent 4 separate sessions. In sessions 1 & 3 as well in 2 & 4 they did the same task, respectively. The task consisted of 2 conditions (M or I). So in session 1 they did “M”, in session 2 “I”, in session 3 “M” and in session 4 “I”.

I now want to estimate some contrasts of interest by combining the effects of sessions 1 + 3 and 2 + 4, respectively.

I use the following code but this creates contrasts that do not separate / model separately the effects in 1 + 3 or 2 + 4 but rather combines all sessions’ effect together. How can I correctly model that?

filedirectory = fullfile(source_directory, 'SPM.mat');
matlabbatch{1}.spm.stats.con.spmmat = {filedirectory};

matlabbatch{1}.spm.stats.con.consess{1}.tcon.name = 'PriorPress';
matlabbatch{1}.spm.stats.con.consess{1}.tcon.weights = [1 0];
matlabbatch{1}.spm.stats.con.consess{1}.tcon.sessrep = 'both';

matlabbatch{1}.spm.stats.con.consess{2}.tcon.name = 'RespondPeriod';
matlabbatch{1}.spm.stats.con.consess{2}.tcon.weights = [0 1];
matlabbatch{1}.spm.stats.con.consess{2}.tcon.sessrep = 'both';

matlabbatch{1}.spm.stats.con.consess{3}.tcon.name = 'PriorPress VS Respond';
matlabbatch{1}.spm.stats.con.consess{3}.tcon.weights = [1 -1];
matlabbatch{1}.spm.stats.con.consess{3}.tcon.sessrep = 'both';
matlabbatch{1}.spm.stats.con.delete = 0;

and this produces:

Screenshot 2022-11-22 at 14.17.38
Screenshot 2022-11-22 at 14.17.381046×802 96.9 KB

I could do manually define it like:

matlabbatch{1}.spm.stats.con.consess{1}.tcon.name = 'M_all VS I_all';
matlabbatch{1}.spm.stats.con.consess{1}.tcon.weights = [1 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 -1];
matlabbatch{1}.spm.stats.con.consess{1}.tcon.sessrep = 'none';

but is there a more automatic way?

The design matrix looks like this:

Screenshot 2022-11-22 at 14.12.43
Screenshot 2022-11-22 at 14.12.431510×1788 168 KB

2

Not aware of an automatic way to do it.

What I usually do is create a logical vector by looking in the column names of the design matrix.

Adding a small script based on your example and a slightly modified version of this function:

Should hopefully make your life a tad easier.


filedirectory = fullfile(source_directory, 'SPM.mat');

load(filedirectory);

[cdtName, regIdx, status] = getRegressorIdx('M', SPM)

matlabbatch{1}.spm.stats.con.consess{1}.tcon.name = 'M_all';
matlabbatch{1}.spm.stats.con.consess{1}.tcon.weights = regIdx;
matlabbatch{1}.spm.stats.con.consess{1}.tcon.sessrep = 'none';

function  [cdtName, regIdx, status] = getRegressorIdx(cdtName, SPM)
    %
    % Gets from the SPM structure the regressors index
    % corresponding to the a condition convolved with the
    % canonical HRF. This can also look for non convolved conditions to identify a confound regressor.
    %
    % Throws a warning if there is no regressor for that condition.
    %
    % USAGE::
    %
    %   [cdtName, regIdx, status] = getRegressorIdx(cdtName, SPM)
    %
    % :param cdtName: name of the condition to look for
    % :type cdtName: char or cellstr
    %
    % :param SPM: content of SPM.mat
    % :type SPM: structure
    %
    %
    % :returns:
    %
    % - :cdtName: (char) name of the condition stripped of any eventual
    %             ``'trial_type.'`` prefix
    % - :regIdx:  (logical) vector of the columns of the design matrix
    %             containing the regressor of interest
    % - :status:  (logical) is ``false`` if no regressor was found for that
    %             condition
    %
    %
  
    % (C) Copyright 2022 bidspm developers
  
    % in case the condition is of something like trial_type.foo
    tokens = regexp(cdtName, '\.', 'split');
    if numel(tokens) > 1
      cdtName = tokens{2};
    end
  
    % construct regexp pattern
    convolvedWithCanonicalHRF = ['^.* (' cdtName '\*bf\(1\))$'];
    nonConvolved = ['^.* ' cdtName '$'];
  
    pattern = [convolvedWithCanonicalHRF '|' nonConvolved];
  
    regIdx = regexp(SPM.xX.name', pattern, 'match');
    regIdx = ~cellfun('isempty', regIdx);
  
    status = checkRegressorFound(regIdx, cdtName);
  
  end
  
  function status = checkRegressorFound(regIdx, cdtName)
  
    status = true;
    regIdx = find(regIdx);
  
    if all(~regIdx)
      status = false;
      msg = sprintf('No regressor found for condition ''%s''', cdtName);
      warning(msg);
    end
  
  end
3

This is nice thanks for sharing!

1 Like