[ENH] allow to run substeps of bms (#1145)

* allow to run substeps of bms

* update changelog

CHANGELOG.md

@@ -24,7 +24,10 @@ 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 option to concatenate runs at subject level to facilite running PPI analysis #1133 by @Remi-Gau
+* [ENH] allow to run substeps of substeps of the bayesian model selection #1145 by @Remi-Gau
+* [ENH] add quality control for GLM using the MACS toolbox to give a goodness of fit several other information criteria (AIC, BIC) #1135  by @Remi-Gau
+* [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
@@ -77,7 +80,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * [ENH] make it possible to update bidspm from any folder in [948](https://github.com/cpp-lln-lab/bidspm/pull/948) by [Remi-Gau](https://github.com/Remi-Gau)
 * [ENH] use filtering of layout in [944](https://github.com/cpp-lln-lab/bidspm/pull/944) by [Remi-Gau](https://github.com/Remi-Gau)
 * [ENH] make reports more silent in [939](https://github.com/cpp-lln-lab/bidspm/pull/939) by [Remi-Gau](https://github.com/Remi-Gau)
-* [ENH] Add CC0 license by default to all outputs in [898](https://github.com/cpp-lln-lab/bidspm/pull/898) by [Remi-Gau](https://github.com/Remi-Gau)
+* [ENH] add CC0 license by default to all outputs in [898](https://github.com/cpp-lln-lab/bidspm/pull/898) by [Remi-Gau](https://github.com/Remi-Gau)
 
 ### Changed
 
@@ -90,7 +93,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 
-* [FIX] Fix vismotion demo in [1070](https://github.com/cpp-lln-lab/bidspm/pull/1070) by [Remi-Gau](https://github.com/Remi-Gau)
+* [FIX] fix vismotion demo in [1070](https://github.com/cpp-lln-lab/bidspm/pull/1070) by [Remi-Gau](https://github.com/Remi-Gau)
 * [FIX] Exclude from GLM specification events with onsets longer than the run duration in [1060](https://github.com/cpp-lln-lab/bidspm/pull/1060) by [Remi-Gau](https://github.com/Remi-Gau)
 * [FIX] remove dummies from preproc dataset and not raw dataset when using CLI in [1057](https://github.com/cpp-lln-lab/bidspm/pull/1057) by [Remi-Gau](https://github.com/Remi-Gau)
 * [FIX] skip smoothing when running bidspm prepoc in dryRun in [1054](https://github.com/cpp-lln-lab/bidspm/pull/1054) by [Remi-Gau](https://github.com/Remi-Gau)

bidspm.m

@@ -85,7 +85,7 @@
     case {'stats', 'contrasts', 'results', 'specify_only'}
       cliStats(varargin{2:end});
 
-    case {'bms'}
+    case {'bms', 'bms-posterior', 'bms-bms'}
       cliBayesModel(varargin{2:end});
 
     case 'meaning_of_life'
@@ -375,7 +375,9 @@ function update()
            'contrasts'; ...
            'results'; ...
            'specify_only', ...
-           'bms'};
+           'bms'; ...
+           'bms-posterior'; ...
+           'bms-bms'};
 
 end
 

demos/bayes/Makefile

@@ -1,11 +1,21 @@
 .PHONY: clean data data_ds000001
 
-data_ds000114:
+data_install:
 	mkdir -p inputs
 	cd inputs && datalad install ///openneuro/ds000114
 	cd inputs && datalad install ///openneuro-derivatives/ds000114-fmriprep
 	cd inputs/ds000114-fmriprep && datalad get sub-*/anat/*MNI152NLin2009cAsym*desc-preproc*.nii.gz -J 12
+
+data_ds000114: data_install
+	cd inputs/ds000114-fmriprep && datalad get sub-*/anat/*MNI152NLin2009cAsym*desc-preproc*.nii.gz -J 12
 	cd inputs/ds000114-fmriprep && datalad get sub-*/ses-*/func/*tsv -J 12
 	cd inputs/ds000114-fmriprep && datalad get sub-*/ses-*/func/*json -J 12
 	cd inputs/ds000114-fmriprep && datalad get sub-*/ses-*/func/*overtverbgeneration*MNI152NLin2009cAsym*_mask.nii.gz -J 12
 	cd inputs/ds000114-fmriprep && datalad get sub-*/ses-*/func/*overtverbgeneration*MNI152NLin2009cAsym*desc-preproc*bold.nii.gz -J 12
+
+data_ds000114_testing: data_install
+	cd inputs/ds000114-fmriprep && datalad get sub-0[12]/anat/*MNI152NLin2009cAsym*desc-preproc*.nii.gz -J 12
+	cd inputs/ds000114-fmriprep && datalad get sub-0[12]/ses-*/func/*tsv -J 12
+	cd inputs/ds000114-fmriprep && datalad get sub-0[12]/ses-*/func/*json -J 12
+	cd inputs/ds000114-fmriprep && datalad get sub-0[12]/ses-*/func/*overtverbgeneration*MNI152NLin2009cAsym*_mask.nii.gz -J 12
+	cd inputs/ds000114-fmriprep && datalad get sub-0[12]/ses-*/func/*overtverbgeneration*MNI152NLin2009cAsym*desc-preproc*bold.nii.gz -J 12

demos/bayes/ds000114_run.m

@@ -6,7 +6,7 @@
 addpath(fullfile(pwd, '..', '..'));
 bidspm();
 
-VERBOSITY = 2;
+VERBOSITY = 0;
 
 FWHM = 8;
 
@@ -16,6 +16,9 @@
 % set to false to not re run the model specification
 FIRST_LEVEL = true;
 
+% set to false to not compute cross-validated log model evidence
+CVLME = true;
+
 % set to true to run on fewer subjects and fewer models
 TESTING = true;
 
@@ -46,7 +49,7 @@
          'options', opt); %#ok<*UNRCH>
 end
 
-%% create models from a default one
+%% create models family from a default one
 
 default_model_file = fullfile(models_dir, 'default_model.json');
 
@@ -66,7 +69,7 @@
 %% Statistics
 preproc_dir = fullfile(output_dir, 'bidspm-preproc');
 
-%% Subject level analysis
+% Subject level analysis
 if FIRST_LEVEL
 
   % Silence this warning as this dataset has not been slice time corrected.
@@ -84,8 +87,52 @@
 
 end
 
+% Run bayesian model selection
+% 1. MA_model_space:    defines a model space
+% 2. MA_cvLME_auto:     computes cross-validated log model evidence
+% 3. MS_PPs_group_auto: calculate posterior probabilities from cvLMEs
+% 4. MS_BMS_group_auto: perform cross-validated Bayesian model selection
+% 5. MS_SMM_BMS:        generate selected models maps from BMS
+if CVLME
+
+  bidspm(bids_dir, output_dir, 'subject', ...
+         'action', 'bms', ...
+         'participant_label', participant_label, ...
+         'models_dir', models_dir, ...
+         'fwhm', FWHM, ...
+         'skip_validation', true, ...
+         'verbosity', VERBOSITY);
+
+end
+
+%% Redefine the model space
+% note that it must be a subset of the one defined previously
+clear multiverse;
+delete(fullfile(pwd, 'models', 'model*.json'));
+
+multiverse.motion = {'basic', 'full'};
+multiverse.scrub = {false};
+multiverse.non_steady_state = {true};
+
+createModelFamilies(default_model_file, multiverse, models_dir);
+
+%%
+% Runs a new bayesian model selection
+% but rely on CVLME estimated previously
+% 1. MA_model_space:    defines a model space
+% 2. MS_PPs_group_auto: calculate posterior probabilities from cvLMEs
+bidspm(bids_dir, output_dir, 'subject', ...
+       'action', 'bms-posterior', ...
+       'participant_label', participant_label, ...
+       'models_dir', models_dir, ...
+       'fwhm', FWHM, ...
+       'skip_validation', true, ...
+       'verbosity', VERBOSITY);
+
+% 3. MS_BMS_group_auto: perform cross-validated Bayesian model selection
+% 4. MS_SMM_BMS:        generate selected models maps from BMS
 bidspm(bids_dir, output_dir, 'subject', ...
-       'action', 'bms', ...
+       'action', 'bms-bms', ...
        'participant_label', participant_label, ...
        'models_dir', models_dir, ...
        'fwhm', FWHM, ...

src/cli/cliBayesModel.m

@@ -18,13 +18,21 @@ function cliBayesModel(varargin)
 
   validate(args);
 
+  action = args.Results.action;
   opt = getOptionsFromCliArgument(args);
   opt.pipeline.type = 'stats';
   opt.pipeline.isBms = true;
   opt = checkOptions(opt);
 
   saveOptions(opt);
 
-  bidsModelSelection(opt, 'action', 'all');
+  switch action
+    case 'bms'
+      bidsModelSelection(opt, 'action', 'all');
+    case 'bms-posterior'
+      bidsModelSelection(opt, 'action', 'posterior');
+    case 'bms-bms'
+      bidsModelSelection(opt, 'action', 'BMS');
+  end
 
 end

src/cli/getOptionsFromCliArgument.m

@@ -70,7 +70,7 @@
 
   opt = optionsStats(opt, args, action);
 
-  if ismember(lower(action), {'bms'})
+  if ismember(lower(action), {'bms', 'bms-posterior', 'bms-bms'})
     opt.toolbox.MACS.model.dir = args.Results.models_dir;
   end
 
@@ -140,7 +140,10 @@
            'contrasts'; ...
            'results'; ...
            'specify_only'; ...
-           'bms'};
+           'bms'; ...
+           'bms-posterior'; ...
+           'bms-bms' ...
+          };
 end
 
 function opt = getOptions(args)

src/cli/inputParserForBayesModel.m

@@ -10,7 +10,6 @@
   isLogical = @(x) islogical(x) && numel(x) == 1;
   isPositiveScalar = @(x) isnumeric(x) && numel(x) == 1 && x >= 0;
   isFolder = @(x) isdir(x);
-  isEmptyOrCellstr = @(x) isempty(x) || iscellstr(x);  %#ok<*ISCLSTR>
 
   addParameter(args, 'models_dir', pwd, isFolder);
 

src/messages/bidspmHelp.m

@@ -388,6 +388,16 @@ function bidspmHelp()
   %          'dry_run', false, ...
   %          'skip_validation', false)
   %
+  % :param action:       Any of: ``{'bms', 'bms-posterior', 'bms-bms'}``
+  %                      `'bms'` will performm all steps for the baeysian
+  %                      model selection. If `'bms'` has been performed
+  %                      than ``'bms-posterior'`` can``'bms-bms'``
+  %                      be performed one aftr the other set a new model
+  %                      space and perform bayesian model selection on it.
+  %                      See the help section of :func:`bidsModelSelection`
+  %                      for more details.
+  % :type  action:       char
+  %
   % :param models_dir:  A path to a dir can be passed as well.
   %                     In this case all ``*_smdl.json`` files will be used
   %                     and looped over.

tests/tests_slow/tests_workflows/stats/test_bidsModelSelection.m

@@ -1,8 +1,5 @@
 function test_suite = test_bidsModelSelection %#ok<*STOUT>
-  %
-
   % (C) Copyright 2022 bidspm developers
-
   try % assignment of 'localfunctions' is necessary in Matlab >= 2016
     test_functions = localfunctions(); %#ok<*NASGU>
   catch % no problem; early Matlab versions can use initTestSuite fine