Merge pull request #585 from Remi-Gau/parametric_and_f_test

[ENH] Add parametric modulation to run / subject level GLM

.github/workflows/system_tests_matlab.m

@@ -5,7 +5,7 @@
 
 addpath(fullfile(root_dir, 'spm12'));
 
-cd demos/MoAE;
+cd(fullfile(root_dir, 'demos', 'MoAE'));
 download_moae_ds(true);
 
 cd(fullfile(root_dir, 'manualTests'));

demos/face_repetition/face_rep_02_stats.m

@@ -19,7 +19,7 @@
 opt = face_rep_get_option_results();
 opt.space = 'IXI549Space';
 
-% bidsFFX('specifyAndEstimate', opt);
+bidsFFX('specifyAndEstimate', opt);
 bidsFFX('contrasts', opt);
 
 bidsResults(opt);

demos/face_repetition/face_rep_02_stats_parametric.m

@@ -0,0 +1,28 @@
+% This script will run the FFX and contrasts
+% on the the face repetition dataset from SPM.
+%
+% Results might be a bit different from those in the manual as some
+% default options are slightly different in this pipeline (e.g use of FAST
+% instead of AR(1)...)
+%
+% (C) Copyright 2019 Remi Gau
+
+% TODO
+% - create the contrast similar to the tuto
+
+clear;
+clc;
+
+addpath(fullfile(pwd, '..', '..'));
+cpp_spm();
+
+opt = face_rep_get_option_results();
+opt.space = 'IXI549Space';
+
+opt.model.file = spm_file(opt.model.file, 'basename', 'model-faceRepetitionParametric_smdl');
+opt.model.bm = BidsModel('file', opt.model.file);
+
+bidsFFX('specifyAndEstimate', opt);
+% bidsFFX('contrasts', opt);
+
+% bidsResults(opt);

demos/face_repetition/models/model-faceRepetitionParametric_smdl.json

@@ -0,0 +1,117 @@
+{
+   "Name": "parametric modulation",
+   "BIDSModelVersion": "1.0.0",
+   "Description": "model for face repetition",
+   "Input": {
+      "task": [
+         "facerepetition"
+      ]
+   },
+   "Nodes": [
+      {
+         "Level": "Run",
+         "Name": "run_level",
+         "GroupBy": [
+            "run",
+            "subject"
+         ],
+         "Transformations": {
+            "Description": "merge the familiarity and repetition column to create the trial type column",
+            "Transformer": "cpp_spm",
+            "Instructions": [
+               {
+                  "Name": "Concatenate",
+                  "Input": [
+                     "face_type",
+                     "repetition_type"
+                  ],
+                  "Output": "trial_type"
+               },
+               {
+                  "Name": "Copy",
+                  "Input": "lag",
+                  "Output": "pmod_lag"
+               },
+               {
+                  "Name": "Power",
+                  "Input": "lag",
+                  "Value": 2,
+                  "Output": "pmod_lag_squared"
+               }
+            ]
+         },
+         "Model": {
+            "X": [
+               "trial_type.famous_1",
+               "trial_type.famous_2",
+               "trial_type.unfamiliar_1",
+               "trial_type.unfamiliar_2",
+               "trans_?",
+               "rot_?"
+            ],
+            "HRF": {
+               "Variables": [
+                  "trial_type.famous_1",
+                  "trial_type.famous_2",
+                  "trial_type.unfamiliar_1",
+                  "trial_type.unfamiliar_2"
+               ],
+               "Model": "DoubleGamma"
+            },
+            "Type": "glm",
+            "Options": {
+               "HighPassFilterCutoffHz": 0.0078
+            },
+            "Software": {
+               "SPM": {
+                  "SerialCorrelation": "AR(1)",
+                  "HRFderivatives": "TemporalDispersion"
+               }
+            }
+         },
+         "DummyContrasts": {
+            "Test": "t",
+            "Contrasts": [
+               "trial_type.famous_1",
+               "trial_type.famous_2",
+               "trial_type.unfamiliar_1",
+               "trial_type.unfamiliar_2"
+            ]
+         },
+         "Contrasts": [
+            {
+               "Name": "faces_gt_baseline",
+               "ConditionList": [
+                  "trial_type.famous_1",
+                  "trial_type.famous_2",
+                  "trial_type.unfamiliar_1",
+                  "trial_type.unfamiliar_2"
+               ],
+               "Weights": [
+                  1,
+                  1,
+                  1,
+                  1
+               ],
+               "Test": "t"
+            },
+            {
+               "Name": "faces_lt_baseline",
+               "ConditionList": [
+                  "trial_type.famous_1",
+                  "trial_type.famous_2",
+                  "trial_type.unfamiliar_1",
+                  "trial_type.unfamiliar_2"
+               ],
+               "Weights": [
+                  -1,
+                  -1,
+                  -1,
+                  -1
+               ],
+               "Test": "t"
+            }
+         ]
+      }
+   ]
+}

src/batches/stats/setBatchEstimateModel.m

@@ -1,15 +1,24 @@
 function matlabbatch = setBatchEstimateModel(matlabbatch, opt, nodeName, contrastsList)
   %
-  % Short description of what the function does goes here.
+  % Set up the estimate model batch for run/subject or group level GLM
   %
   % USAGE::
   %
-  %   matlabbatch = setBatchEstimateModel(matlabbatch, grpLvlCon)
+  %   matlabbatch = setBatchEstimateModel(matlabbatch, opt)
+  %   matlabbatch = setBatchEstimateModel(matlabbatch, opt, nodeName, contrastsList)
   %
   % :param matlabbatch:
   % :type matlabbatch: structure
-  % :param grpLvlCon:
-  % :type grpLvlCon:
+  %
+  % :param opt:
+  % :type opt: structure
+  %
+  % :param nodeName:
+  % :type nodeName: char
+  %
+  % :param contrastsList:
+  % :type contrastsList: cell string
+  %
   %
   % :returns: - :matlabbatch: (structure)
   %

src/batches/stats/setBatchSubjectLevelGLMSpec.m

@@ -16,7 +16,7 @@
   % :type opt: structure
   %
   % :param subLabel:
-  % :type subLabel: string
+  % :type subLabel: char
   %
   % :returns: - :matlabbatch: (structure)
   %
@@ -228,8 +228,7 @@
 
   % get events file from raw data set and convert it to a onsets.mat file
   % store in the subject level GLM directory
-  filter = struct( ...
-                  'sub',  subLabel, ...
+  filter = struct('sub',  subLabel, ...
                   'task', task, ...
                   'ses', session, ...
                   'run', run, ...

src/subject_level/convertOnsetTsvToMat.m

@@ -2,11 +2,6 @@
   %
   % Converts an events.tsv file to an onset file suitable for SPM subject level analysis.
   %
-  % The function extracts from the events.tsv file the trials (with type, onsets, and durations)
-  % of the conditions of interest as requested in the model.json.
-  % It then stores them in the GLM folder in a .mat file
-  % that can be fed directly in an SPM GLM batch.
-  %
   % USAGE::
   %
   %   fullpathOnsetFilename = convertOnsetTsvToMat(opt, tsvFile)
@@ -17,6 +12,26 @@
   % :param tsvFile:
   % :type tsvFile: string
   %
+  % Use a BIDS stats model specified in a JSON file to:
+  %
+  % - loads events.tsv and apply the Node.Transformations to its content
+  %
+  % - extract the trials (onsets, durations)
+  %   of the conditions that should be convolved
+  %   as requested from Node.HRF.Variables field
+  %
+  % It then stores them in in a .mat file
+  % that can be fed directly in an SPM GLM batch as 'Multiple conditions'
+  %
+  % Parametric modulation can be specified via columns in the TSV
+  % file starting with ``pmod_``.
+  % These columns can be created via the use of Node.Transformations.
+  % Only polynomial 1 are supported.
+  % More complex modulation should be precomputed via the Transformations.
+  %
+  % if ``opt.glm.useDummyRegressor`` is set to ``true``, any missing condition
+  % will be replaced by a DummyRegressor.
+  %
   % :returns: :fullpathOnsetFilename: (string) name of the output ``.mat`` file.
   %
   % See also: createAndReturnOnsetFile, bids.transformers
@@ -52,13 +67,17 @@
   names = {};
   onsets = {};
   durations = {};
+  pmod = struct('name', {''}, 'param', {}, 'poly', {});
 
   if ~isfield(opt.model, 'bm')
     opt.model.bm = BidsModel('file', opt.model.file);
   end
+  % TODO get / apply transformers from a specific node
   transformers = opt.model.bm.getBidsTransformers();
   tsvContent = bids.transformers(transformers, tsvContent);
 
+  conditionIdx = 1;
+
   for iCond = 1:numel(variablesToConvolve)
 
     trialTypeNotFound = false;
@@ -74,22 +93,23 @@
       trialTypes = tsvContent.(tokens{1});
       conditionName = strjoin(tokens(2:end), '.');
 
-      idx = find(strcmp(conditionName, trialTypes));
+      rows = find(strcmp(conditionName, trialTypes));
 
       printToScreen(sprintf('   Condition %s: %i trials found.\n', ...
                             conditionName, ...
-                            numel(idx)), ...
+                            numel(rows)), ...
                     opt);
 
-      if ~isempty(idx)
+      if ~isempty(rows)
 
         if ~ismember(variablesToConvolve{iCond}, designMatrix)
           continue
         end
 
-        names{1, end + 1} = conditionName;
-        onsets{1, end + 1} = tsvContent.onset(idx)'; %#ok<*AGROW,*NASGU>
-        durations{1, end + 1} = tsvContent.duration(idx)';
+        names{1, conditionIdx} = conditionName;
+        onsets{1, conditionIdx} = tsvContent.onset(rows)'; %#ok<*AGROW,*NASGU>
+        durations{1, conditionIdx} = tsvContent.duration(rows)';
+        pmod = parametricModulation(pmod, tsvContent, rows, conditionIdx);
 
       else
 
@@ -124,6 +144,8 @@
 
     end
 
+    conditionIdx = conditionIdx + 1;
+
   end
 
   %% save the onsets as a matfile
@@ -136,11 +158,38 @@
   fullpathOnsetFilename = fullfile(pth, bf.filename);
 
   save(fullpathOnsetFilename, ...
-       'names', 'onsets', 'durations', ...
+       'names', 'onsets', 'durations', 'pmod', ...
        '-v7');
 
 end
 
+function pmod = parametricModulation(pmod, tsvContent, rows, conditionIdx)
+  % parametric modulation (pmod)
+  %
+  % skipped if parametric modulation is 1 for all onsets
+  %
+  % coerces NaNs into 1
+
+  fields = fieldnames(tsvContent);
+  pmodIdx = ~cellfun('isempty', regexp(fields, '^pmod_.*', 'match'));
+  pmodIdx = find(pmodIdx);
+
+  for iMod = 1:numel(pmodIdx)
+
+    thisMod = fields{pmodIdx(iMod)};
+
+    amplitude = tsvContent.(thisMod)(rows);
+    amplitude(isnan(amplitude)) = 1;
+
+    if ~all(amplitude == 1)
+      pmod(1, conditionIdx).name{iMod}  = strrep(thisMod, 'pmod_', '');
+      pmod(end).param{iMod} = amplitude;
+      pmod(end).poly{iMod}  = 1;
+    end
+
+  end
+end
+
 function [names, onsets, durations] = addDummyRegressor(names, onsets, durations)
 
   names{1, end + 1} = 'dummyRegressor';

src/subject_level/createAndReturnOnsetFile.m

@@ -8,19 +8,22 @@
   %
   % USAGE::
   %
-  %   onsetFilename = createAndReturnOnsetFile(opt, subLabel, tsvFile, funcFWHM)
+  %   onsetFilename = createAndReturnOnsetFile(opt, subLabel, tsvFile)
   %
   % :param opt:
   % :type opt: structure
+  %
   % :param subLabel:
-  % :type subLabel: string
+  % :type subLabel: char
+  %
   % :param tsvFile: fullpath name of the tsv file.
-  % :type tsvFile: string
+  % :type tsvFile: char
   %
-  % :returns: :onsetFilename: (string) fullpath name of the file created.
+  % :returns: :onsetFilename: (path) fullpath name of the file created.
   %
   % See also: convertOnsetTsvToMat
   %
+  %
   % (C) Copyright 2019 CPP_SPM developers
 
   if iscell(tsvFile)