[ENH] Add several confound strategies to add to bids stats model (#1126)

* add confound strategy * create family models * update changelog * fix doc * fix bugs * skip warning
cpp-lln-lab · Aug 13, 2023 · 03c1bae · 03c1bae
1 parent 35397c0
commit 03c1bae
Show file tree

Hide file tree

Showing 12 changed files with 552 additions and 97 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -24,6 +24,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+* [ENH] Add several confound strategies to add to bids stats model and a function to create family of models #1126 by @Remi-Gau
+* [ENH] add CLI to run bayesian model selection #1121 by @Remi-Gau
 * [ENH] support label of activations with all atlases [1100](https://github.com/cpp-lln-lab/bidspm/pull/1100) by [Remi-Gau](https://github.com/Remi-Gau)
 * [ENH] add support for session level models #1116 be @Remi-Gau
 

diff --git a/demos/bayes/ds000114_run.m b/demos/bayes/ds000114_run.m
@@ -6,17 +6,17 @@
 addpath(fullfile(pwd, '..', '..'));
 bidspm();
 
+VERBOSITY = 2;
+
+FWHM = 8;
+
 % set to false to not re run the smoothing
-SMOOTH = true;
+SMOOTH = false;
 
 % set to false to not re run the model specification
 FIRST_LEVEL = true;
 
-VERBOSITY = 1;
-
-FWHM = 8;
-
-% to run on fewer subjects
+% set to true to run on fewer subjects and fewer models
 TESTING = true;
 
 % The directory where the data are located
@@ -30,6 +30,7 @@
 participant_label = {'[0-9]*'}; %#ok<*NASGU>
 if TESTING
   participant_label = {'^0[12]$'};
+  participant_label = {'[0-9]*'};
 end
 
 %% Smooth
@@ -50,28 +51,28 @@
 
 default_model_file = fullfile(models_dir, 'default_model.json');
 
-mutliverse.strategy = {'motion', 'wm_csf', 'scrub', 'non_steady_state'};
-mutliverse.motion = {'none', 'basic', 'full'};
-mutliverse.scrub = [false, true];
-mutliverse.wm_csf = {'none', 'basic', 'full'};
-mutliverse.non_steady_state = [false, true];
+multiverse.motion = {'none', 'basic', 'full'};
+multiverse.scrub = {false, true};
+multiverse.wm_csf = {'none', 'basic', 'full'};
+multiverse.non_steady_state = {false, true};
 
 if TESTING
-  mutliverse.strategy = {'motion', 'wm_csf', 'scrub', 'non_steady_state'};
-  mutliverse.motion = {'none', 'basic'};
-  mutliverse.scrub = [false, true];
-  mutliverse.wm_csf = {'none'};
-  mutliverse.non_steady_state = false;
+  multiverse.motion = {'basic', 'full'};
+  multiverse.scrub = {false, true};
+  multiverse.non_steady_state = {true};
 end
 
-create_model_families(models_dir, default_model_file, mutliverse);
+createModelFamilies(default_model_file, multiverse, models_dir);
 
 %% Statistics
 preproc_dir = fullfile(output_dir, 'bidspm-preproc');
 
 %% Subject level analysis
 if FIRST_LEVEL
 
+  % Silence this warning as this dataset has not been slice time corrected.
+  warning('OFF', 'setBatchSubjectLevelGLMSpec:noSliceTimingInfoForGlm');
+
   bidspm(bids_dir, output_dir, 'subject', ...
          'participant_label', participant_label, ...
          'action', 'specify_only', ...
@@ -91,74 +92,3 @@
        'fwhm', FWHM, ...
        'skip_validation', true, ...
        'verbosity', VERBOSITY);
-
-%%
-function create_model_families(models_dir, default_model_file, mutliverse)
-  % create models from a default one
-  %
-
-  % TODO incorporate into bidspm
-
-  % TODO add support for 12 motion regressors
-  strategyToSkip = fieldnames(mutliverse);
-  idxStrategyToSkip = ~ismember(fieldnames(mutliverse), mutliverse.strategy);
-  strategyToSkip = strategyToSkip(idxStrategyToSkip);
-  for i = 1:numel(strategyToSkip)
-    mutliverse.(strategyToSkip{i}) = {''};
-  end
-
-  for i = 1:numel(mutliverse.motion)
-    for j = 1:numel(mutliverse.scrub)
-      for k = 1:numel(mutliverse.wm_csf)
-        for l = 1:numel(mutliverse.non_steady_state)
-
-          model = bids.util.jsondecode(default_model_file);
-
-          name = sprintf('rp-%s_scrub-%i_tissue-%s_nsso-%i', ...
-                         mutliverse.motion{i}, ...
-                         mutliverse.scrub(j), ...
-                         mutliverse.wm_csf{k}, ...
-                         mutliverse.non_steady_state(l));
-          model.Name = name;
-          model.Nodes.Name = name;
-
-          design_matrix = model.Nodes.Model.X;
-
-          switch mutliverse.motion{i}
-            case 'none'
-            case 'basic'
-              design_matrix{end + 1} = 'rot_?';
-              design_matrix{end + 1} = 'trans_?';
-            case 12
-            case 'full'
-              design_matrix{end + 1} = 'rot_*';
-              design_matrix{end + 1} = 'trans_*';
-          end
-
-          if mutliverse.scrub(j) == 1
-            design_matrix{end + 1} = 'motion_outlier*'; %#ok<*AGROW>
-          end
-
-          switch mutliverse.wm_csf{k}
-            case 'none'
-            case 'basic'
-              design_matrix{end + 1} = 'csf';
-              design_matrix{end + 1} = 'white';
-            case 'full'
-              design_matrix{end + 1} = 'csf_*';
-              design_matrix{end + 1} = 'white_*';
-          end
-
-          if mutliverse.non_steady_state(l)
-            design_matrix{end + 1} = 'non_steady_state_outlier*';
-          end
-
-          model.Nodes.Model.X = design_matrix;
-
-          output_file = fullfile(models_dir, ['model_' name '_smdl.json']);
-          bids.util.jsonencode(output_file, model);
-        end
-      end
-    end
-  end
-end
diff --git a/demos/bayes/models/default_model.json b/demos/bayes/models/default_model.json
@@ -34,6 +34,9 @@
         "Options": {
           "HighPassFilterCutoffHz": 0.008,
           "Mask": {
+            "ses": [
+              "test"
+            ],
             "desc": [
               "brain"
             ],

diff --git a/docs/source/dev_doc.rst b/docs/source/dev_doc.rst
@@ -274,6 +274,8 @@ QA
 .. autofunction:: src.QA.anatQA
 .. _censoring:
 .. autofunction:: src.QA.censoring
+.. _compileScrubbingStats:
+.. autofunction:: src.QA.compileScrubbingStats
 .. _computeDesignEfficiency:
 .. autofunction:: src.QA.computeDesignEfficiency
 .. _computeFDandRMS:
@@ -342,12 +344,16 @@ bids
 
 bids_model
 ==========
+.. _addConfoundsToDesignMatrix:
+.. autofunction:: src.bids_model.addConfoundsToDesignMatrix
 .. _checkContrast:
 .. autofunction:: src.bids_model.checkContrast
 .. _checkGroupBy:
 .. autofunction:: src.bids_model.checkGroupBy
 .. _createDefaultStatsModel:
 .. autofunction:: src.bids_model.createDefaultStatsModel
+.. _createModelFamilies:
+.. autofunction:: src.bids_model.createModelFamilies
 .. _getContrastsFromParentNode:
 .. autofunction:: src.bids_model.getContrastsFromParentNode
 .. _getContrastsList:
@@ -368,6 +374,8 @@ cli
 ===
 .. _baseInputParser:
 .. autofunction:: src.cli.baseInputParser
+.. _cliBayesModel:
+.. autofunction:: src.cli.cliBayesModel
 .. _cliCopy:
 .. autofunction:: src.cli.cliCopy
 .. _cliCreateRoi:
@@ -384,6 +392,8 @@ cli
 .. autofunction:: src.cli.getBidsFilterFile
 .. _getOptionsFromCliArgument:
 .. autofunction:: src.cli.getOptionsFromCliArgument
+.. _inputParserForBayesModel:
+.. autofunction:: src.cli.inputParserForBayesModel
 .. _inputParserForCopy:
 .. autofunction:: src.cli.inputParserForCopy
 .. _inputParserForCreateModel:
@@ -686,6 +696,8 @@ utils
 .. autofunction:: src.utils.createDataDictionary
 .. _deregexify:
 .. autofunction:: src.utils.deregexify
+.. _displayArguments:
+.. autofunction:: src.utils.displayArguments
 .. _getDist2surf:
 .. autofunction:: src.utils.getDist2surf
 .. _getFuncVoxelDims:

diff --git a/src/QA/compileScrubbingStats.m b/src/QA/compileScrubbingStats.m
@@ -1,6 +1,6 @@
 function compileScrubbingStats(statsFolder)
   %
-  % Make a list of *_desc-confounds_timeseries.json
+  % Make a list of ``*_desc-confounds_timeseries.json``
   % and compile their results in a single tsv.
   %
   % EXAMPLE::