diff --git a/CHANGELOG.md b/CHANGELOG.md index 397aff4a..dc3fbed9 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] make `addConfoundsToDesignMatrix` a method of `BidsModel` #1294 by @Remi-Gau * [ENH] add Apptainer definition #1254 by @Remi-Gau and @monique2208 * [ENH] allow to copy anat only on raw datasets #1181 by @Remi-Gau * [ENH] add option to concatenate runs at subject level to facilite running PPI analysis #1133 by @Remi-Gau diff --git a/docs/source/API.rst b/docs/source/API.rst index 0a73f918..8e189bfa 100644 --- a/docs/source/API.rst +++ b/docs/source/API.rst @@ -188,7 +188,6 @@ bids bids_model ========== -.. autofunction:: src.bids_model.addConfoundsToDesignMatrix .. autofunction:: src.bids_model.checkContrast .. autofunction:: src.bids_model.checkGroupBy .. autofunction:: src.bids_model.createDefaultStatsModel diff --git a/src/bids_model/BidsModel.m b/src/bids_model/BidsModel.m index 1def9928..3722a2e1 100644 --- a/src/bids_model/BidsModel.m +++ b/src/bids_model/BidsModel.m @@ -346,6 +346,156 @@ end + function obj = addConfoundsToDesignMatrix(obj, varargin) + % + % Add some typical confounds to the design matrix of bids stat model. + % + % This will update the design matrix of the root node of the model. + % + % Similar to the :func:`nilearn.interfaces.fmriprep.load_confounds` + % + % USAGE:: + % + % bm = bm.addConfoundsToDesignMatrix('strategy', strategy); + % + % + % :param bm: bids stats model. + % :type bm: :obj:`BidsModel` instance or path + % to a ``_smdl.json`` file + % + % :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 = bm.addConfoundsToDesignMatrix('strategy', strategy); + % + % + + % (C) Copyright 2023 bidspm developers + + args = inputParser; + args.CaseSensitive = false; + args.KeepUnmatched = false; + args.FunctionName = 'addConfoundsToDesignMatrix'; + + addParameter(args, 'strategy', defaultStrategy(), @isstruct); + addParameter(args, 'updateName', false, @islogical); + + parse(args, varargin{:}); + + strategy = args.Results.strategy; + strategy = setFieldsStrategy(strategy); + + [~, name] = obj.get_root_node(); + [~, idx] = obj.get_nodes('Name', name); + designMatrix = obj.get_design_matrix('Name', name); + + strategiesToApply = strategy.strategies; + for i = 1:numel(strategiesToApply) + + switch strategiesToApply{i} + + case 'motion' + switch strategy.motion{1} + 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{1} + 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{1} + 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.\n', ... + 'Supported strategies are:%s'], ... + strategiesToApply{i}, ... + bids.internal.create_unordered_list(supportedStrategies()))); + otherwise + logger('WARNING', sprintf('Unknown strategey: "%s".', ... + strategiesToApply{i}), ... + 'filename', mfilename(), ... + 'id', 'unknownStrategy'); + end + end + + designMatrix = cleanDesignMatrix(designMatrix); + + obj.Nodes{idx}.Model.X = designMatrix; + + if args.Results.updateName + obj.Nodes{idx}.Name = appendSuffixToNodeName(obj.Nodes{idx}.Name, strategy); + end + + end + function validateConstrasts(obj) % validate all contrasts spec in the model @@ -425,3 +575,71 @@ function bidsModelError(obj, id, msg) 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) + % remove empty and duplicate + toClean = cellfun(@(x) isempty(x), designMatrix); + designMatrix(toClean) = []; + + if isempty(designMatrix) + return + end + if size(designMatrix, 1) > 1 + designMatrix = designMatrix'; + end + + numeric = cellfun(@(x) isnumeric(x), designMatrix); + tmp = unique(designMatrix(~numeric)); + + designMatrix = cat(2, tmp, designMatrix(numeric)); +end + +function strategy = setFieldsStrategy(strategy) + + tmp = defaultStrategy(); + + strategies = fieldnames(defaultStrategy()); + for i = 1:numel(strategies) + + if ~isfield(strategy, strategies{i}) + strategy.(strategies{i}) = tmp.(strategies{i}); + end + + if ~iscell(strategy.(strategies{i})) + strategy.(strategies{i}) = {strategy.(strategies{i})}; + end + + if ~isempty(strategy.(strategies{i})) && ... + isnumeric(strategy.(strategies{i}){1}) && ... + isnan(strategy.(strategies{i}){1}) + strategy.(strategies{i}){1} = tmp.(strategies{i}); + end + + end + +end diff --git a/src/bids_model/addConfoundsToDesignMatrix.m b/src/bids_model/addConfoundsToDesignMatrix.m deleted file mode 100644 index ac97c4a5..00000000 --- a/src/bids_model/addConfoundsToDesignMatrix.m +++ /dev/null @@ -1,220 +0,0 @@ -function bm = addConfoundsToDesignMatrix(varargin) - % - % Add some typical confounds to the design matrix of bids stat model. - % - % This will update the design matrix of the root node of the model. - % - % Similar to the :func:`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 - % - % :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 - - args = inputParser; - args.CaseSensitive = false; - args.KeepUnmatched = false; - args.FunctionName = 'addConfoundsToDesignMatrix'; - - isBidsModelOrFile = @(x) isa(x, 'BidsModel') || exist(x, 'file') == 2; - addRequired(args, 'bm', isBidsModelOrFile); - - addParameter(args, 'strategy', defaultStrategy(), @isstruct); - addParameter(args, 'updateName', false, @islogical); - - 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{1} - 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{1} - 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{1} - 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.\n', ... - 'Supported strategies are:%s'], ... - strategiesToApply{i}, ... - bids.internal.create_unordered_list(supportedStrategies()))); - otherwise - logger('WARNING', sprintf('Unknown strategey: "%s".', ... - strategiesToApply{i}), ... - 'filename', mfilename(), ... - 'id', 'unknownStrategy'); - end - end - - designMatrix = cleanDesignMatrix(designMatrix); - - 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) - % remove empty and duplicate - toClean = cellfun(@(x) isempty(x), designMatrix); - designMatrix(toClean) = []; - - if isempty(designMatrix) - return - end - if size(designMatrix, 1) > 1 - designMatrix = designMatrix'; - end - - numeric = cellfun(@(x) isnumeric(x), designMatrix); - tmp = unique(designMatrix(~numeric)); - - designMatrix = cat(2, tmp, designMatrix(numeric)); -end - -function strategy = setFieldsStrategy(strategy) - - tmp = defaultStrategy(); - - strategies = fieldnames(defaultStrategy()); - for i = 1:numel(strategies) - - if ~isfield(strategy, strategies{i}) - strategy.(strategies{i}) = tmp.(strategies{i}); - end - - if ~iscell(strategy.(strategies{i})) - strategy.(strategies{i}) = {strategy.(strategies{i})}; - end - - if ~isempty(strategy.(strategies{i})) && ... - isnumeric(strategy.(strategies{i}){1}) && ... - 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 index cf4c7571..0faf25c4 100644 --- a/src/bids_model/createModelFamilies.m +++ b/src/bids_model/createModelFamilies.m @@ -39,7 +39,7 @@ function createModelFamilies(varargin) args = inputParser; args.CaseSensitive = false; args.KeepUnmatched = false; - args.FunctionName = 'addConfoundsToDesignMatrix'; + args.FunctionName = 'createModelFamilies'; isBidsModelOrFile = @(x) isa(x, 'BidsModel') || exist(x, 'file') == 2; @@ -81,7 +81,7 @@ function createModelFamilies(varargin) 'non_steady_state', multiverse.non_steady_state{l}); bm = defaultModel; - bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy, 'updateName', true); + bm = bm.addConfoundsToDesignMatrix('strategy', strategy, 'updateName', true); name = bm.Nodes{idx}.Name; bm.Name = name; diff --git a/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m b/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m index 20ed2436..14c771c9 100644 --- a/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m +++ b/tests/tests_bids_model/test_addConfoundsToDesignMatrix.m @@ -11,32 +11,14 @@ function test_addConfoundsToDesignMatrix_default() bm = BidsModel('init', true); - bm = addConfoundsToDesignMatrix(bm); + bm = bm.addConfoundsToDesignMatrix(); 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', ... - 'model-balloonanalogrisktask_smdl.json'); - - strategy.strategies = {'motion', 'non_steady_state'}; - strategy.motion = 'basic'; - strategy.non_steady_state = true; - - bm = addConfoundsToDesignMatrix(modelFile, 'strategy', strategy); - + bm = bm.addConfoundsToDesignMatrix('strategy', strategy); assertEqual(numel(bm.Nodes{1}.Model.X), 3); assertEqual(ismember({'rot_?', 'trans_?', 'non_steady_state_outlier*'}, ... bm.Nodes{1}.Model.X), ... @@ -51,7 +33,7 @@ function test_addConfoundsToDesignMatrix_with_numerial_and_no_duplicate() strategy.strategies = {'motion'}; strategy.motion = 'basic'; - bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy); + bm = bm.addConfoundsToDesignMatrix('strategy', strategy); assertEqual(numel(bm.Nodes{1}.Model.X), 3); @@ -72,7 +54,7 @@ function test_addConfoundsToDesignMatrix_full() strategy.wm_csf = 'full'; strategy.non_steady_state = true; - bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy); + bm = bm.addConfoundsToDesignMatrix('strategy', strategy); assertEqual(numel(bm.Nodes{1}.Model.X), 6); @@ -95,10 +77,10 @@ function test_addConfoundsToDesignMatrix_update_name() strategy.wm_csf = 'full'; strategy.non_steady_state = true; - bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy); + bm = bm.addConfoundsToDesignMatrix('strategy', strategy); assertEqual(bm.Nodes{1}.Name, nameBefore); - bm = addConfoundsToDesignMatrix(bm, 'strategy', strategy, 'updateName', true); + bm = bm.addConfoundsToDesignMatrix('strategy', strategy, 'updateName', true); assertEqual(bm.Nodes{1}.Name, [nameBefore '_rp-full_scrub-1_tissue-full_nsso-1']); end @@ -112,11 +94,11 @@ function test_addConfoundsToDesignMatrix_warning() strategy.strategies = {'foo'}; - assertWarning(@()addConfoundsToDesignMatrix(bm, 'strategy', strategy), ... - 'addConfoundsToDesignMatrix:unknownStrategy'); + assertWarning(@()bm.addConfoundsToDesignMatrix('strategy', strategy), ... + 'BidsModel:unknownStrategy'); strategy.strategies = {'global_signal'}; - assertWarning(@()addConfoundsToDesignMatrix(bm, 'strategy', strategy), ... - 'addConfoundsToDesignMatrix:notImplemented'); + assertWarning(@()bm.addConfoundsToDesignMatrix('strategy', strategy), ... + 'BidsModel:notImplemented'); end