From 300b62cd898e079c5b6d47ccd3739866da9f8802 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Sat, 12 Aug 2023 13:49:53 -0400 Subject: [PATCH 1/6] add confound strategy --- demos/bayes/ds000114_run.m | 54 +++---- src/bids_model/addConfoundsToDesignMatrix.m | 139 ++++++++++++++++++ .../test_addConfoundsToDesignMatrix.m | 93 ++++++++++++ 3 files changed, 259 insertions(+), 27 deletions(-) create mode 100644 src/bids_model/addConfoundsToDesignMatrix.m create mode 100644 tests/tests_bids_model/test_addConfoundsToDesignMatrix.m diff --git a/demos/bayes/ds000114_run.m b/demos/bayes/ds000114_run.m index 8a95fb2f0..0c851aab2 100644 --- a/demos/bayes/ds000114_run.m +++ b/demos/bayes/ds000114_run.m @@ -50,21 +50,21 @@ 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.strategy = {'motion', 'wm_csf', 'scrub', 'non_steady_state'}; +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.strategy = {'motion', 'wm_csf', 'scrub', 'non_steady_state'}; + multiverse.motion = {'none', 'basic'}; + multiverse.scrub = [false, true]; + multiverse.wm_csf = {'none'}; + multiverse.non_steady_state = false; end -create_model_families(models_dir, default_model_file, mutliverse); +create_model_families(models_dir, default_model_file, multiverse); %% Statistics preproc_dir = fullfile(output_dir, 'bidspm-preproc'); @@ -93,38 +93,38 @@ 'verbosity', VERBOSITY); %% -function create_model_families(models_dir, default_model_file, mutliverse) +function create_model_families(models_dir, default_model_file, multiverse) % 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 = fieldnames(multiverse); + idxStrategyToSkip = ~ismember(fieldnames(multiverse), multiverse.strategy); strategyToSkip = strategyToSkip(idxStrategyToSkip); for i = 1:numel(strategyToSkip) - mutliverse.(strategyToSkip{i}) = {''}; + multiverse.(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) + for i = 1:numel(multiverse.motion) + for j = 1:numel(multiverse.scrub) + for k = 1:numel(multiverse.wm_csf) + for l = 1:numel(multiverse.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)); + multiverse.motion{i}, ... + multiverse.scrub(j), ... + multiverse.wm_csf{k}, ... + multiverse.non_steady_state(l)); model.Name = name; model.Nodes.Name = name; design_matrix = model.Nodes.Model.X; - switch mutliverse.motion{i} + switch multiverse.motion{i} case 'none' case 'basic' design_matrix{end + 1} = 'rot_?'; @@ -135,11 +135,11 @@ function create_model_families(models_dir, default_model_file, mutliverse) design_matrix{end + 1} = 'trans_*'; end - if mutliverse.scrub(j) == 1 + if multiverse.scrub(j) == 1 design_matrix{end + 1} = 'motion_outlier*'; %#ok<*AGROW> end - switch mutliverse.wm_csf{k} + switch multiverse.wm_csf{k} case 'none' case 'basic' design_matrix{end + 1} = 'csf'; @@ -149,7 +149,7 @@ function create_model_families(models_dir, default_model_file, mutliverse) design_matrix{end + 1} = 'white_*'; end - if mutliverse.non_steady_state(l) + if multiverse.non_steady_state(l) design_matrix{end + 1} = 'non_steady_state_outlier*'; end diff --git a/src/bids_model/addConfoundsToDesignMatrix.m b/src/bids_model/addConfoundsToDesignMatrix.m new file mode 100644 index 000000000..b1080035e --- /dev/null +++ b/src/bids_model/addConfoundsToDesignMatrix.m @@ -0,0 +1,139 @@ +function bm = addConfoundsToDesignMatrix(varargin) + % + % Add some typical confounds to the design matrix of bids stat model. + % + % Similar to the module nilearn.interfaces.fmriprep.load_confounds. + % + % USAGE:: + % + % bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy); + % + % + % :param bm: bids stats model. + % :type bm: :obj:`BidsModel` instance + % or path to a ``_smdl.json`` file. + % + % :param strategy: structure with the fields: + % ``strategy``, + % ``motion``, ``non_steady_state``, ``scrub``, ``wm_csf`` + % :type strategy: struct + % + + % (C) Copyright 2023 bidspm developers + + args = inputParser; + args.CaseSensitive = false; + args.KeepUnmatched = false; + args.FunctionName = 'addConfoundsToDesignMatrix'; + + isBidsModelOrFile = @(x) isa(x, 'BidsModel') || exist(x, 'file') == 2; + addRequired(args, 'bm', isBidsModelOrFile); + + defaultStrategy.strategies = {'motion', 'non_steady_state'}; + addParameter(args, 'strategy', defaultStrategy, @isstruct); + + parse(args, varargin{:}); + + bm = args.Results.bm; + if ischar(bm) + bm = BidsModel('file', bm); + end + + strategy = args.Results.strategy; + strategy = setFieldsStrategy(strategy); + + [~, name] = bm.get_root_node(); + [~, idx] = bm.get_nodes('Name', name); + designMatrix = bm.get_design_matrix('Name', name); + + strategiesToApply = strategy.strategies; + for i = 1:numel(strategiesToApply) + + switch strategiesToApply{i} + + case 'motion' + switch strategy.motion + case 'none' + case 'basic' + designMatrix{end + 1} = 'rot_?'; %#ok<*AGROW> + designMatrix{end + 1} = 'trans_?'; + case {'power2', 'derivatives' } + notImplemented(mfilename(), ... + sprintf('motion "%s" not implemented.', strategy.motion)); + case 'full' + designMatrix{end + 1} = 'rot_*'; + designMatrix{end + 1} = 'trans_*'; + end + + case 'non_steady_state' + if strategy.non_steady_state + designMatrix{end + 1} = 'non_steady_state_outlier*'; + end + + case 'scrub' + if strategy.scrub == 1 + designMatrix{end + 1} = 'motion_outlier*'; + end + + case 'wm_csf' + switch strategy.wm_csf + case 'none' + case 'basic' + designMatrix{end + 1} = 'csf'; + designMatrix{end + 1} = 'white'; + case 'full' + designMatrix{end + 1} = 'csf_*'; + designMatrix{end + 1} = 'white_*'; + otherwise + notImplemented(mfilename(), ... + sprintf('wm_csf "%s" not implemented.', strategiesToApply{i})); + end + + case {'global_signal', 'compcorstr', 'n_compcorstr'} + notImplemented(mfilename(), ... + sprintf('Strategey "%s" not implemented.', strategiesToApply{i})); + otherwise + logger('WARNING', sprintf('Unknown strategey: "%s".', ... + strategiesToApply{i}), ... + 'filename', mfilename(), ... + 'id', 'unknownStrategy'); + end + end + + designMatrix = cleanDesignMatrix(designMatrix); + + bm.Nodes{idx}.Model.X = designMatrix; + +end + +function designMatrix = cleanDesignMatrix(designMatrix) + % remove empty and duplicate + toClean = cellfun(@(x) isempty(x), designMatrix); + designMatrix(toClean) = []; + + numeric = cellfun(@(x) isnumeric(x), designMatrix); + tmp = unique(designMatrix(~numeric)); + if size(tmp, 1) > 1 + tmp = tmp'; + end + if size(designMatrix, 1) > 1 + designMatrix = designMatrix'; + end + designMatrix = cat(2, tmp, designMatrix(numeric)); +end + +function strategy = setFieldsStrategy(strategy) + + defaultStrategy.motion = 'basic'; + defaultStrategy.scrub = false; + defaultStrategy.wm_csf = 'none'; + defaultStrategy.non_steady_state = true; + + strategies = fieldnames(defaultStrategy); + for i = 1:numel(strategies) + if ~isfield(strategy, strategies{i}) + strategy.(strategies{i}) = defaultStrategy.(strategies{i}); + end + end + +end diff --git a/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m b/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m new file mode 100644 index 000000000..35521d97f --- /dev/null +++ b/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m @@ -0,0 +1,93 @@ +function test_suite = test_addConfoundsToDesignMatrix %#ok<*STOUT> + % (C) Copyright 2023 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 + end + initTestSuite; +end + +function test_addConfoundsToDesignMatrix_from_file() + + modelFile = fullfile(getTestDataDir(), 'models', ... + 'model-balloonanalogrisktask_smdl.json'); + + bm = addConfoundsToDesignMatrix(modelFile); + + assertEqual(numel(bm.Nodes{1}.Model.X), 3); + assertEqual(ismember({'rot_?', 'trans_?', 'non_steady_state_outlier*'}, ... + bm.Nodes{1}.Model.X), ... + true(1, 3)); + +end + +function test_addConfoundsToDesignMatrix_with_numerial_and_no_duplicate() + + bm = BidsModel('init', true); + bm.Nodes{1}.Model.X = {1, 'rot_?'}; + + bm = addConfoundsToDesignMatrix(bm); + + assertEqual(numel(bm.Nodes{1}.Model.X), 4); + + numeric = cellfun(@(x) isnumeric(x), bm.Nodes{1}.Model.X); + assertEqual(ismember({'rot_?', 'trans_?', 'non_steady_state_outlier*'}, ... + bm.Nodes{1}.Model.X(~numeric)), ... + true(1, 3)); + +end + +function test_addConfoundsToDesignMatrix_full() + + bm = BidsModel('init', true); + + strategy.strategies = {'motion', 'wm_csf', 'scrub', 'non_steady_state'}; + strategy.motion = 'full'; + strategy.scrub = true; + strategy.wm_csf = 'full'; + strategy.non_steady_state = true; + + bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy); + + bm.Nodes{1}.Model.X; + + assertEqual(numel(bm.Nodes{1}.Model.X), 6); + + assertEqual(ismember({'rot_*', 'trans_*', ... + 'csf_*', 'white_*', ... + 'motion_outlier*', 'non_steady_state_outlier*'}, ... + bm.Nodes{1}.Model.X), ... + true(1, 6)); +end + +function test_addConfoundsToDesignMatrix_default() + + bm = BidsModel('init', true); + + bm = addConfoundsToDesignMatrix(bm); + + assertEqual(numel(bm.Nodes{1}.Model.X), 3); + assertEqual(ismember({'rot_?', 'trans_?', 'non_steady_state_outlier*'}, ... + bm.Nodes{1}.Model.X), ... + true(1, 3)); + +end + +function test_addConfoundsToDesignMatrix_warning() + + if bids.internal.is_octave() + moxunit_throw_test_skipped_exception('testing warning with octave'); + end + + bm = BidsModel('init', true); + + strategy.strategies = {'foo'}; + + assertWarning(@()addConfoundsToDesignMatrix(bm, 'strategy', strategy), ... + 'addConfoundsToDesignMatrix:unknownStrategy'); + + strategy.strategies = {'global_signal'}; + + assertWarning(@()addConfoundsToDesignMatrix(bm, 'strategy', strategy), ... + 'addConfoundsToDesignMatrix:notImplemented'); +end From 13a71ab808dbabcaae5fb3a1ef41e896847e3de4 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Sat, 12 Aug 2023 18:14:34 -0400 Subject: [PATCH 2/6] create family models --- CHANGELOG.md | 1 + demos/bayes/ds000114_run.m | 82 +------------ demos/bayes/models/default_model.json | 3 + docs/source/dev_doc.rst | 12 ++ src/QA/compileScrubbingStats.m | 2 +- src/bids_model/addConfoundsToDesignMatrix.m | 116 +++++++++++++++--- src/bids_model/createModelFamilies.m | 100 +++++++++++++++ src/defaults/checkOptions.m | 2 +- src/messages/bidspmHelp.m | 17 +-- .../test_addConfoundsToDesignMatrix.m | 43 +++++-- .../test_createModelFamilies.m | 60 +++++++++ 11 files changed, 325 insertions(+), 113 deletions(-) create mode 100644 src/bids_model/createModelFamilies.m create mode 100644 tests/tests_bids_model/test_createModelFamilies.m diff --git a/CHANGELOG.md b/CHANGELOG.md index f4ed76663..3e4e2c652 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +* [ENH] add CLI to run bayesian model selection #1121 be @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 index 0c851aab2..9a66aa35e 100644 --- a/demos/bayes/ds000114_run.m +++ b/demos/bayes/ds000114_run.m @@ -50,21 +50,18 @@ default_model_file = fullfile(models_dir, 'default_model.json'); -multiverse.strategy = {'motion', 'wm_csf', 'scrub', 'non_steady_state'}; multiverse.motion = {'none', 'basic', 'full'}; -multiverse.scrub = [false, true]; +multiverse.scrub = {false, true}; multiverse.wm_csf = {'none', 'basic', 'full'}; -multiverse.non_steady_state = [false, true]; +multiverse.non_steady_state = {false, true}; if TESTING - multiverse.strategy = {'motion', 'wm_csf', 'scrub', 'non_steady_state'}; multiverse.motion = {'none', 'basic'}; - multiverse.scrub = [false, true]; multiverse.wm_csf = {'none'}; - multiverse.non_steady_state = false; + multiverse.non_steady_state = {false}; end -create_model_families(models_dir, default_model_file, multiverse); +createModelFamilies(default_model_file, multiverse, models_dir); %% Statistics preproc_dir = fullfile(output_dir, 'bidspm-preproc'); @@ -91,74 +88,3 @@ 'fwhm', FWHM, ... 'skip_validation', true, ... 'verbosity', VERBOSITY); - -%% -function create_model_families(models_dir, default_model_file, multiverse) - % create models from a default one - % - - % TODO incorporate into bidspm - - % TODO add support for 12 motion regressors - strategyToSkip = fieldnames(multiverse); - idxStrategyToSkip = ~ismember(fieldnames(multiverse), multiverse.strategy); - strategyToSkip = strategyToSkip(idxStrategyToSkip); - for i = 1:numel(strategyToSkip) - multiverse.(strategyToSkip{i}) = {''}; - end - - for i = 1:numel(multiverse.motion) - for j = 1:numel(multiverse.scrub) - for k = 1:numel(multiverse.wm_csf) - for l = 1:numel(multiverse.non_steady_state) - - model = bids.util.jsondecode(default_model_file); - - name = sprintf('rp-%s_scrub-%i_tissue-%s_nsso-%i', ... - multiverse.motion{i}, ... - multiverse.scrub(j), ... - multiverse.wm_csf{k}, ... - multiverse.non_steady_state(l)); - model.Name = name; - model.Nodes.Name = name; - - design_matrix = model.Nodes.Model.X; - - switch multiverse.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 multiverse.scrub(j) == 1 - design_matrix{end + 1} = 'motion_outlier*'; %#ok<*AGROW> - end - - switch multiverse.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 multiverse.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 index 99cb28f6b..a69ccc6bb 100644 --- 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 index b8ed7d075..de0b6f133 100644 --- 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 index 1280c54b0..7a590d4fd 100644 --- 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:: diff --git a/src/bids_model/addConfoundsToDesignMatrix.m b/src/bids_model/addConfoundsToDesignMatrix.m index b1080035e..e12c1f354 100644 --- a/src/bids_model/addConfoundsToDesignMatrix.m +++ b/src/bids_model/addConfoundsToDesignMatrix.m @@ -2,7 +2,10 @@ % % Add some typical confounds to the design matrix of bids stat model. % - % Similar to the module nilearn.interfaces.fmriprep.load_confounds. + % This will update the design matrix of the root node of the model. + % + % Similar to the module + % https://nilearn.github.io/dev/modules/generated/nilearn.interfaces.fmriprep.load_confounds.html % % USAGE:: % @@ -10,13 +13,49 @@ % % % :param bm: bids stats model. - % :type bm: :obj:`BidsModel` instance - % or path to a ``_smdl.json`` file. + % :type bm: :obj:`BidsModel` instance or path to a ``_smdl.json`` file % - % :param strategy: structure with the fields: - % ``strategy``, - % ``motion``, ``non_steady_state``, ``scrub``, ``wm_csf`` % :type strategy: struct + % :param strategy: structure describing the confoudd strategy. + % + % The structure must have the following field: + % + % - ``strategy``: cell array of char with the strategies to apply. + % + % The structure may have the following field: + % + % - ``motion``: motion regressors strategy + % - ``scrub``: scrubbing strategy + % - ``wm_csf``: white matter and cerebrospinal fluid regressors strategy + % - ``non_steady_state``: non steady state regressors strategy + % + % See the nilearn documentation (mentioned above) + % for more information on the possible values those strategies can take. + % + % :type updateName: logical + % :param updateName: Append the name of the root node + % with a string describing the counfounds added. + % + % ``rp-{motion}_scrub-{scrub}_tissue-{wm_csf}_nsso-{non_steady_state}`` + % + % default = ``false`` + % + % + % :rtype: :obj:`BidsModel` instance + % :return: bids stats model with the confounds added. + % + % EXAMPLE: + % + % .. code-block:: matlab + % + % + % strategy.strategies = {'motion', 'wm_csf', 'scrub', 'non_steady_state'}; + % strategy.motion = 'full'; + % strategy.scrub = true; + % strategy.non_steady_state = true; + % + % bm = addConfoundsToDesignMatrix(path_to_statsmodel_file, 'strategy', strategy); + % % % (C) Copyright 2023 bidspm developers @@ -29,8 +68,8 @@ isBidsModelOrFile = @(x) isa(x, 'BidsModel') || exist(x, 'file') == 2; addRequired(args, 'bm', isBidsModelOrFile); - defaultStrategy.strategies = {'motion', 'non_steady_state'}; - addParameter(args, 'strategy', defaultStrategy, @isstruct); + addParameter(args, 'strategy', defaultStrategy(), @isstruct); + addParameter(args, 'updateName', false, @islogical); parse(args, varargin{:}); @@ -52,7 +91,7 @@ switch strategiesToApply{i} case 'motion' - switch strategy.motion + switch strategy.motion{1} case 'none' case 'basic' designMatrix{end + 1} = 'rot_?'; %#ok<*AGROW> @@ -66,17 +105,17 @@ end case 'non_steady_state' - if strategy.non_steady_state + if strategy.non_steady_state{1} designMatrix{end + 1} = 'non_steady_state_outlier*'; end case 'scrub' - if strategy.scrub == 1 + if strategy.scrub{1} designMatrix{end + 1} = 'motion_outlier*'; end case 'wm_csf' - switch strategy.wm_csf + switch strategy.wm_csf{1} case 'none' case 'basic' designMatrix{end + 1} = 'csf'; @@ -91,7 +130,10 @@ case {'global_signal', 'compcorstr', 'n_compcorstr'} notImplemented(mfilename(), ... - sprintf('Strategey "%s" not implemented.', strategiesToApply{i})); + sprintf(['Strategey "%s" not implemented.\n', ... + 'Supported strategies are:%s'], ... + strategiesToApply{i}, ... + bids.internal.create_unordered_list(supportedStrategies()))); otherwise logger('WARNING', sprintf('Unknown strategey: "%s".', ... strategiesToApply{i}), ... @@ -104,6 +146,35 @@ bm.Nodes{idx}.Model.X = designMatrix; + if args.Results.updateName + bm.Nodes{idx}.Name = appendSuffixToNodeName(bm.Nodes{idx}.Name, strategy); + end + +end + +function name = appendSuffixToNodeName(name, strategy) + if ~isempty(name) + name = [name, '_']; + end + suffix = sprintf('rp-%s_scrub-%i_tissue-%s_nsso-%i', ... + strategy.motion{1}, ... + strategy.scrub{1}, ... + strategy.wm_csf{1}, ... + strategy.non_steady_state{1}); + + name = [name suffix]; +end + +function value = supportedStrategies() + value = {'motion', 'non_steady_state', 'wm_csf', 'scrub'}; +end + +function value = defaultStrategy() + value.strategies = {}; + value.motion = 'none'; + value.scrub = false; + value.wm_csf = 'none'; + value.non_steady_state = false; end function designMatrix = cleanDesignMatrix(designMatrix) @@ -124,16 +195,23 @@ function strategy = setFieldsStrategy(strategy) - defaultStrategy.motion = 'basic'; - defaultStrategy.scrub = false; - defaultStrategy.wm_csf = 'none'; - defaultStrategy.non_steady_state = true; + tmp = defaultStrategy(); - strategies = fieldnames(defaultStrategy); + strategies = fieldnames(defaultStrategy()); for i = 1:numel(strategies) + if ~isfield(strategy, strategies{i}) - strategy.(strategies{i}) = defaultStrategy.(strategies{i}); + strategy.(strategies{i}) = tmp.(strategies{i}); + end + + if ~iscell(strategy.(strategies{i})) + strategy.(strategies{i}) = {strategy.(strategies{i})}; + end + + if isnan(strategy.(strategies{i}){1}) + strategy.(strategies{i}){1} = tmp.(strategies{i}); end + end end diff --git a/src/bids_model/createModelFamilies.m b/src/bids_model/createModelFamilies.m new file mode 100644 index 000000000..9b103ffb7 --- /dev/null +++ b/src/bids_model/createModelFamilies.m @@ -0,0 +1,100 @@ +function createModelFamilies(varargin) + % + % Create a family of models from a default one. + % + % USAGE:: + % + % createModelFamilies(defaultModel, multiverse, outputDir) + % + % + % :param defaultModel: bids stats model that serves as template. + % :type defaultModel: :obj:`BidsModel` instance or path to a ``_smdl.json`` file + % + % :param multiverse: Structure to describe the multiverse of models. + % + % Each field of the structure is a dimension of the multiverse. + % Possible dimensions are: + % + % - ``motion``: motion regressors strategy + % - ``scrub``: scrubbing strategy + % - ``wm_csf``: white matter and cerebrospinal fluid regressors strategy + % - ``non_steady_state``: non steady state regressors strategy + % :type multiverse: struct + % + % EXAMPLE: + % + % .. code-block:: matlab + % + % multiverse.motion = {'none', 'basic', 'full'}; + % multiverse.scrub = {false, true}; + % multiverse.wm_csf = {'none', 'basic', 'full'}; + % multiverse.non_steady_state = {false, true}; + % + % createModelFamilies(path_to_statsmodel_file, multiverse, output_path); + % + % + + % (C) Copyright 2023 bidspm developers + + args = inputParser; + args.CaseSensitive = false; + args.KeepUnmatched = false; + args.FunctionName = 'addConfoundsToDesignMatrix'; + + isBidsModelOrFile = @(x) isa(x, 'BidsModel') || exist(x, 'file') == 2; + + addRequired(args, 'defaultModel', isBidsModelOrFile); + addRequired(args, 'multiverse', @isstruct); + addRequired(args, 'outputDir', @isdir); + + parse(args, varargin{:}); + + defaultModel = args.Results.defaultModel; + multiverse = args.Results.multiverse; + outputDir = args.Results.outputDir; + + if ischar(defaultModel) + defaultModel = BidsModel('file', defaultModel); + end + + [~, name] = defaultModel.get_root_node(); + [~, idx] = defaultModel.get_nodes('Name', name); + defaultModel.Nodes{idx}.Name = ''; + + missingDimensionIdx = ~ismember(supportedDimensions(), fieldnames(multiverse)); + missingDimension = supportedDimensions(); + missingDimension = missingDimension(missingDimensionIdx); + + for i = 1:numel(missingDimension) + multiverse.(missingDimension{i}) = {nan}; + end + + for i = 1:numel(multiverse.motion) + for j = 1:numel(multiverse.scrub) + for k = 1:numel(multiverse.wm_csf) + for l = 1:numel(multiverse.non_steady_state) + + strategy = struct('strategies', {fieldnames(multiverse)}, ... + 'motion', multiverse.motion{i}, ... + 'scrub', multiverse.scrub{j}, ... + 'wm_csf', multiverse.wm_csf{k}, ... + 'non_steady_state', multiverse.non_steady_state{l}); + + bm = defaultModel; + bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy, 'updateName', true); + + name = bm.Nodes{idx}.Name; + bm.Name = name; + output_file = fullfile(outputDir, ['model_' name '_smdl.json']); + + bm.write(output_file); + + end + end + end + end +end + +function value = supportedDimensions() + value = {'motion', 'non_steady_state', 'wm_csf', 'scrub'}; +end diff --git a/src/defaults/checkOptions.m b/src/defaults/checkOptions.m index 7e4995832..c3660ec0d 100644 --- a/src/defaults/checkOptions.m +++ b/src/defaults/checkOptions.m @@ -111,7 +111,7 @@ % ``p(grayMatter) + p(whiteMatter) + p(CSF) > threshold`` % will be included in the mask. % - ``opt.skullstrip.do = true`` - - % Set to ``true`` to skip skullstripping. + % Set to ``true`` to skip skullstripping. % % - ``opt.stc.skip = false`` - % Boolean flag to skip slice time correction or not. diff --git a/src/messages/bidspmHelp.m b/src/messages/bidspmHelp.m index ee882c2b6..f2b16d030 100644 --- a/src/messages/bidspmHelp.m +++ b/src/messages/bidspmHelp.m @@ -46,7 +46,7 @@ function bidspmHelp() % - ``'default_model'``: creates a default BIDS stats model % - ``'create_roi'``: creates ROIs from a given atlas % - ``'stats'``: runs model specification / estimation, - % contrast computation, display results + % contrast computation, display results % - ``'contrasts'``: runs contrast computation, display results % - ``'results'``: displays results % - ``'bms'``: performs bayesian model selection @@ -188,13 +188,16 @@ function bidspmHelp() % :type space: cell string % % :param roi_atlas: Can be any of: - % - ``'visfatlas'`` - % - ``'anatomy_toobox'`` - % - ``'neuromorphometrics'`` - % - ``'hcpex'`` - % - ``'wang'`` - % - ``'glasser'`` + + % - ``'visfatlas'`` + % - ``'anatomy_toobox'`` + % - ``'neuromorphometrics'`` + % - ``'hcpex'`` + % - ``'wang'`` + % - ``'glasser'`` + % Defaults to ``'neuromorphometrics'`` + % :type roi_atlas: char % % :param roi_name: Name of the roi to create. If the ROI does not exist in the atlas, diff --git a/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m b/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m index 35521d97f..2e3c830b5 100644 --- a/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m +++ b/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m @@ -12,7 +12,11 @@ function test_addConfoundsToDesignMatrix_from_file() modelFile = fullfile(getTestDataDir(), 'models', ... 'model-balloonanalogrisktask_smdl.json'); - bm = addConfoundsToDesignMatrix(modelFile); + strategy.strategies = {'motion', 'non_steady_state'}; + strategy.motion = 'basic'; + strategy.non_steady_state = true; + + bm = addConfoundsToDesignMatrix(modelFile, 'strategy', strategy); assertEqual(numel(bm.Nodes{1}.Model.X), 3); assertEqual(ismember({'rot_?', 'trans_?', 'non_steady_state_outlier*'}, ... @@ -26,14 +30,16 @@ function test_addConfoundsToDesignMatrix_with_numerial_and_no_duplicate() bm = BidsModel('init', true); bm.Nodes{1}.Model.X = {1, 'rot_?'}; - bm = addConfoundsToDesignMatrix(bm); + strategy.strategies = {'motion'}; + strategy.motion = 'basic'; + bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy); - assertEqual(numel(bm.Nodes{1}.Model.X), 4); + assertEqual(numel(bm.Nodes{1}.Model.X), 3); numeric = cellfun(@(x) isnumeric(x), bm.Nodes{1}.Model.X); - assertEqual(ismember({'rot_?', 'trans_?', 'non_steady_state_outlier*'}, ... + assertEqual(ismember({'rot_?', 'trans_?'}, ... bm.Nodes{1}.Model.X(~numeric)), ... - true(1, 3)); + true(1, 2)); end @@ -49,8 +55,6 @@ function test_addConfoundsToDesignMatrix_full() bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy); - bm.Nodes{1}.Model.X; - assertEqual(numel(bm.Nodes{1}.Model.X), 6); assertEqual(ismember({'rot_*', 'trans_*', ... @@ -60,12 +64,37 @@ function test_addConfoundsToDesignMatrix_full() true(1, 6)); end +function test_addConfoundsToDesignMatrix_update_name() + + bm = BidsModel('init', true); + + nameBefore = bm.Nodes{1}.Name; + + strategy.strategies = {'motion', 'wm_csf', 'scrub', 'non_steady_state'}; + strategy.motion = 'full'; + strategy.scrub = true; + strategy.wm_csf = 'full'; + strategy.non_steady_state = true; + + bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy); + assertEqual(bm.Nodes{1}.Name, nameBefore); + + bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy, 'updateName', true); + assertEqual(bm.Nodes{1}.Name, [nameBefore '_rp-full_scrub-1_tissue-full_nsso-1']); +end + function test_addConfoundsToDesignMatrix_default() bm = BidsModel('init', true); bm = addConfoundsToDesignMatrix(bm); + assertEqual(numel(bm.Nodes{1}.Model.X), 0); + strategy.strategies = {'motion', 'non_steady_state'}; + strategy.motion = 'basic'; + strategy.non_steady_state = true; + + bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy); assertEqual(numel(bm.Nodes{1}.Model.X), 3); assertEqual(ismember({'rot_?', 'trans_?', 'non_steady_state_outlier*'}, ... bm.Nodes{1}.Model.X), ... diff --git a/tests/tests_bids_model/test_createModelFamilies.m b/tests/tests_bids_model/test_createModelFamilies.m new file mode 100644 index 000000000..9635b5951 --- /dev/null +++ b/tests/tests_bids_model/test_createModelFamilies.m @@ -0,0 +1,60 @@ +function test_suite = test_createModelFamilies %#ok<*STOUT> + % (C) Copyright 2023 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 + end + initTestSuite; +end + +function test_createModelFamilies_ones_dimension() + + defaultModel = BidsModel('init', true); + + outputDir = tempName; + + multiverse.motion = {'full'}; + multiverse.scrub = {true, false}; + + createModelFamilies(defaultModel, multiverse, outputDir); + + files = spm_select('FPList', outputDir, '.*smdl.json'); + assertEqual(size(files, 1), 2); + +end + +function test_createModelFamilies_basic() + + defaultModel = BidsModel('init', true); + + outputDir = tempName; + + multiverse.motion = {'none'}; + multiverse.scrub = {false}; + multiverse.wm_csf = {'none'}; + multiverse.non_steady_state = {false}; + + createModelFamilies(defaultModel, multiverse, outputDir); + + files = spm_select('FPList', outputDir, '.*smdl.json'); + assertEqual(size(files, 1), 1); + +end + +function test_createModelFamilies_all() + + defaultModel = BidsModel('init', true); + + outputDir = tempName; + + multiverse.motion = {'none', 'basic', 'full'}; + multiverse.scrub = {false, true}; + multiverse.wm_csf = {'none', 'basic', 'full'}; + multiverse.non_steady_state = {false, true}; + + createModelFamilies(defaultModel, multiverse, outputDir); + + files = spm_select('FPList', outputDir, '.*smdl.json'); + assertEqual(size(files, 1), 3 * 2 * 3 * 2); + +end From 09d5e6170ee870d9dfb9c5ed0c839c01c375ef9a Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Sat, 12 Aug 2023 18:18:42 -0400 Subject: [PATCH 3/6] update changelog --- CHANGELOG.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3e4e2c652..56e6abf3b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,7 +24,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added -* [ENH] add CLI to run bayesian model selection #1121 be @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 From 36151fbc033be68547e236584d0d936170196b78 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Sat, 12 Aug 2023 18:23:28 -0400 Subject: [PATCH 4/6] fix doc --- src/messages/bidspmHelp.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/messages/bidspmHelp.m b/src/messages/bidspmHelp.m index f2b16d030..e48c5aba0 100644 --- a/src/messages/bidspmHelp.m +++ b/src/messages/bidspmHelp.m @@ -188,16 +188,16 @@ function bidspmHelp() % :type space: cell string % % :param roi_atlas: Can be any of: - + % % - ``'visfatlas'`` % - ``'anatomy_toobox'`` % - ``'neuromorphometrics'`` % - ``'hcpex'`` % - ``'wang'`` % - ``'glasser'`` - + % % Defaults to ``'neuromorphometrics'`` - + % % :type roi_atlas: char % % :param roi_name: Name of the roi to create. If the ROI does not exist in the atlas, From 69fb49e4c50c0af7326cedabb479fc9f8d6b8515 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Sat, 12 Aug 2023 19:23:34 -0400 Subject: [PATCH 5/6] fix bugs --- demos/bayes/ds000114_run.m | 19 ++++++---- src/bids_model/addConfoundsToDesignMatrix.m | 4 +- .../test_addConfoundsToDesignMatrix.m | 38 +++++++++---------- 3 files changed, 33 insertions(+), 28 deletions(-) diff --git a/demos/bayes/ds000114_run.m b/demos/bayes/ds000114_run.m index 9a66aa35e..3ad4a7423 100644 --- 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; % 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 @@ -56,9 +56,9 @@ multiverse.non_steady_state = {false, true}; if TESTING - multiverse.motion = {'none', 'basic'}; - multiverse.wm_csf = {'none'}; - multiverse.non_steady_state = {false}; + multiverse.motion = {'basic', 'full'}; + multiverse.scrub = {false, true}; + multiverse.non_steady_state = {true}; end createModelFamilies(default_model_file, multiverse, models_dir); @@ -69,6 +69,9 @@ %% 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', ... diff --git a/src/bids_model/addConfoundsToDesignMatrix.m b/src/bids_model/addConfoundsToDesignMatrix.m index e12c1f354..fb9fb01f1 100644 --- a/src/bids_model/addConfoundsToDesignMatrix.m +++ b/src/bids_model/addConfoundsToDesignMatrix.m @@ -208,7 +208,9 @@ strategy.(strategies{i}) = {strategy.(strategies{i})}; end - if isnan(strategy.(strategies{i}){1}) + if ~isempty(strategy.(strategies{i})) && ... + isnumeric(strategy.(strategies{i}){1}) && ... + isnan(strategy.(strategies{i}){1}) strategy.(strategies{i}){1} = tmp.(strategies{i}); end diff --git a/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m b/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m index 2e3c830b5..20ed2436f 100644 --- a/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m +++ b/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m @@ -7,6 +7,25 @@ initTestSuite; end +function test_addConfoundsToDesignMatrix_default() + + bm = BidsModel('init', true); + + bm = addConfoundsToDesignMatrix(bm); + assertEqual(numel(bm.Nodes{1}.Model.X), 0); + + strategy.strategies = {'motion', 'non_steady_state'}; + strategy.motion = 'basic'; + strategy.non_steady_state = true; + + bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy); + assertEqual(numel(bm.Nodes{1}.Model.X), 3); + assertEqual(ismember({'rot_?', 'trans_?', 'non_steady_state_outlier*'}, ... + bm.Nodes{1}.Model.X), ... + true(1, 3)); + +end + function test_addConfoundsToDesignMatrix_from_file() modelFile = fullfile(getTestDataDir(), 'models', ... @@ -83,25 +102,6 @@ function test_addConfoundsToDesignMatrix_update_name() assertEqual(bm.Nodes{1}.Name, [nameBefore '_rp-full_scrub-1_tissue-full_nsso-1']); end -function test_addConfoundsToDesignMatrix_default() - - bm = BidsModel('init', true); - - bm = addConfoundsToDesignMatrix(bm); - assertEqual(numel(bm.Nodes{1}.Model.X), 0); - - strategy.strategies = {'motion', 'non_steady_state'}; - strategy.motion = 'basic'; - strategy.non_steady_state = true; - - bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy); - assertEqual(numel(bm.Nodes{1}.Model.X), 3); - assertEqual(ismember({'rot_?', 'trans_?', 'non_steady_state_outlier*'}, ... - bm.Nodes{1}.Model.X), ... - true(1, 3)); - -end - function test_addConfoundsToDesignMatrix_warning() if bids.internal.is_octave() From f66f858e09d12af28db3eccd1264fcc83b527ba3 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Sat, 12 Aug 2023 20:27:13 -0400 Subject: [PATCH 6/6] skip warning --- demos/bayes/ds000114_run.m | 3 ++- tests/tests_bids/test_getAnatFilename.m | 4 ++++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/demos/bayes/ds000114_run.m b/demos/bayes/ds000114_run.m index 3ad4a7423..cad84817c 100644 --- a/demos/bayes/ds000114_run.m +++ b/demos/bayes/ds000114_run.m @@ -11,7 +11,7 @@ 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; @@ -30,6 +30,7 @@ participant_label = {'[0-9]*'}; %#ok<*NASGU> if TESTING participant_label = {'^0[12]$'}; + participant_label = {'[0-9]*'}; end %% Smooth diff --git a/tests/tests_bids/test_getAnatFilename.m b/tests/tests_bids/test_getAnatFilename.m index 7063e1705..5e6c61bbd 100644 --- a/tests/tests_bids/test_getAnatFilename.m +++ b/tests/tests_bids/test_getAnatFilename.m @@ -31,6 +31,10 @@ function test_getAnatFilename_return_several() assertEqual(numel(anatImage), 3); warning('ON'); + if bids.internal.is_octave() + moxunit_throw_test_skipped_exception('Octave:mixed-string-concat warning thrown'); + end + assertWarning(@()getAnatFilename(BIDS, opt, subLabel, nbImgToReturn), ... 'getAnatFilename:severalAnatFile');