From ab138dedf9f46461615000eba83a04802b47a58e Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Thu, 25 Jul 2024 11:03:44 +0200 Subject: [PATCH] update code and tests --- src/batches/stats/setBatchFactorialDesign.m | 140 ++++++++++++++++-- src/bids_model/BidsModel.m | 5 +- src/bids_model/checkGroupBy.m | 2 +- tests/create_3_groups_dataset.py | 2 +- tests/data/models/model-bug616_smdl.json | 2 +- .../model-vislocalizer2sampleTTest_smdl.json | 2 +- .../model-vismotionNoOverWrite_smdl.json | 2 +- .../model-vismotionOneWayANOVA_smdl.json} | 0 .../test_setBatchFactorialDesign_3_groups.m | 4 +- tests/utils/setOptions.m | 2 +- 10 files changed, 138 insertions(+), 23 deletions(-) rename tests/data/{3_groups/models/model-vismotion_smdl.json => models/model-vismotionOneWayANOVA_smdl.json} (100%) diff --git a/src/batches/stats/setBatchFactorialDesign.m b/src/batches/stats/setBatchFactorialDesign.m index 1fe7e981c..2ceaf7f2b 100644 --- a/src/batches/stats/setBatchFactorialDesign.m +++ b/src/batches/stats/setBatchFactorialDesign.m @@ -36,26 +36,98 @@ return end - [BIDS, opt] = getData(opt, opt.dir.preproc); - printBatchName(sprintf('specify group level fmri model for node "%s"', nodeName), opt); - contrasts = getContrastsListForFactorialDesign(opt, nodeName); - % now we fetch the contrast for each subject and allocate them in the batch % - first case is we pool over all subjects % - second case is we pool over all the subject of a group defined in the % participants.tsv of the raw dataset - if strcmp(glmType, 'one_sample_t_test') - if numel(groupBy) == 1 && strcmpi(groupBy, {'contrast'}) - [matlabbatch, contrastsList, groups] = tTestAcrossSubject(opt, contrasts); - elseif numel(groupBy) == 2 - [matlabbatch, contrastsList, groups] = tTestForGroup(opt, contrasts); - end + switch glmType + case 'one_sample_t_test' + + if numel(groupBy) == 1 && strcmpi(groupBy, {'contrast'}) + [matlabbatch, contrastsList, groups] = tTestAcrossSubject(matlabbatch, opt, nodeName); + elseif numel(groupBy) == 2 + [matlabbatch, contrastsList, groups] = tTestForGroup(matlabbatch, opt, nodeName); + end + + case 'one_way_anova' + + % to keep track of all the contrast we used + % and to which group each batch corresponds to + % the later is needed to be able to create proper RFX dir name + contrastsList = {}; + groups = {}; + + [BIDS, opt] = getData(opt, opt.dir.preproc); + + contrasts = getContrastsListForFactorialDesign(opt, nodeName); + + designMatrix = opt.model.bm.get_design_matrix('Name', nodeName); + designMatrix = cellfun(@(x) num2str(x), designMatrix, 'uniformoutput', false); + + groupColumnHdr = setxor(designMatrix, {'1'}); + groupColumnHdr = groupColumnHdr{1}; + + checkColumnParticipantsTsv(BIDS, groupColumnHdr); + + availableGroups = unique(BIDS.raw.participants.content.(groupColumnHdr)); + + rfxDir = getRFXdir(opt, nodeName, contrasts{1}, '1WayANOVA'); + overwriteDir(rfxDir, opt); + + assert(exist(fullfile(rfxDir, 'SPM.mat'), 'file') == 0); + + matlabbatch = returnOneWayAnovaBatch(matlabbatch, rfxDir); + + mask = getInclusiveMask(opt, nodeName); + matlabbatch{end}.spm.stats.factorial_design.masking.em = {mask}; + + icell(1).scans = {}; + + for iGroup = 1:numel(availableGroups) + + thisGroup = availableGroups{iGroup}; + + % grab subjects label from participants.tsv in raw + % and only keep those that are part of the requested subjects + % + % Note that this will lead to different results + % depending on the requested subejcts + % + tsv = BIDS.raw.participants.content; + subjectsInGroup = strcmp(tsv.(groupColumnHdr), thisGroup); + subjectsLabel = regexprep(tsv.participant_id(subjectsInGroup), '^sub-', ''); + subjectsLabel = intersect(subjectsLabel, opt.subjects); + + % collect all con images from all subjects + for iSub = 1:numel(subjectsLabel) + subLabel = subjectsLabel{iSub}; + conImages{iSub} = findSubjectConImage(opt, subLabel, contrasts); + end + + for iCon = 1:numel(contrasts) + + contrastName = contrasts{iCon}; + + contrastsList{end + 1, 1} = contrastName; + groups{end + 1, 1} = thisGroup; + + msg = sprintf(' Group contrast "%s" for group "%s"', contrastName, thisGroup); + logger('INFO', msg, 'options', opt, 'filename', mfilename()); + + icell = allocateSubjectsContrasts(opt, subjectsLabel, conImages, iCon); + + matlabbatch{end}.spm.stats.factorial_design.des.anova.icell(iGroup).scans = icell.scans; + + end + + end + end end -function [matlabbatch, contrastsList, groups] = tTestAcrossSubject(opt, contrasts) +function [matlabbatch, contrastsList, groups] = tTestAcrossSubject(matlabbatch, opt, nodeName) % to keep track of all the contrast we used % and to which group each batch corresponds to @@ -63,6 +135,8 @@ contrastsList = {}; groups = {}; + contrasts = getContrastsListForFactorialDesign(opt, nodeName); + % collect all con images from all subjects for iSub = 1:numel(opt.subjects) subLabel = opt.subjects{iSub}; @@ -87,7 +161,7 @@ end end -function [matlabbatch, contrastsList, groups] = tTestForGroup(opt, contrasts) +function [matlabbatch, contrastsList, groups] = tTestForGroup(matlabbatch, opt, nodeName) % to keep track of all the contrast we used % and to which group each batch corresponds to @@ -95,7 +169,15 @@ contrastsList = {}; groups = {}; + [BIDS, opt] = getData(opt, opt.dir.preproc); + + contrasts = getContrastsListForFactorialDesign(opt, nodeName); + + node = opt.model.bm.get_nodes('Name', nodeName); + groupBy = node.GroupBy; + groupColumnHdr = setxor(groupBy, {'contrast'}); + groupColumnHdr = groupColumnHdr{1}; checkColumnParticipantsTsv(BIDS, groupColumnHdr); @@ -226,15 +308,45 @@ end +function matlabbatch = returnOneWayAnovaBatch(matlabbatch, directory) + + factorialDesign.dir = {directory}; + + factorialDesign.des.anova.icell(1).scans = {}; + + % Assumes groups are independent + factorialDesign.des.anova.dept = 0; + % 1: Assumes that the variance is not the same across groups + % 0: There is no difference in the variance between groups + factorialDesign.des.anova.variance = 1; + factorialDesign.des.anova.gmsca = 0; + factorialDesign.des.anova.ancova = 0; + + factorialDesign.cov = struct('c', {}, 'cname', {}, 'iCFI', {}, 'iCC', {}); + + factorialDesign.multi_cov = struct('files', {}, 'iCFI', {}, 'iCC', {}); + + factorialDesign.masking.tm.tm_none = 1; + factorialDesign.masking.im = 1; + factorialDesign.masking.em = {''}; + + factorialDesign.globalc.g_omit = 1; + factorialDesign.globalm.gmsca.gmsca_no = 1; + factorialDesign.globalm.glonorm = 1; + + matlabbatch{end + 1}.spm.stats.factorial_design = factorialDesign; + +end + function [status, groupBy, glmType] = checks(opt, nodeName) thisNode = opt.model.bm.get_nodes('Name', nodeName); commonMsg = sprintf('for the dataset level node: "%s"', nodeName); - status = checkGroupBy(thisNode); - participants = bids.util.tsvread(fullfile(opt.dir.raw, 'participants.tsv')); + columns = fieldnames(participants); + status = checkGroupBy(thisNode, columns); [glmType, ~, groupBy] = groupLevelGlmType(opt, nodeName, participants); diff --git a/src/bids_model/BidsModel.m b/src/bids_model/BidsModel.m index 3722a2e1c..16b930b4e 100644 --- a/src/bids_model/BidsModel.m +++ b/src/bids_model/BidsModel.m @@ -90,7 +90,8 @@ [model, nodeName] = obj.getDefaultModel(varargin{:}); - if ~isfield(model.Options, 'HighPassFilterCutoffHz') || ... + if ~isfield(model, 'Options') || ... + ~isfield(model.Options, 'HighPassFilterCutoffHz') || ... isempty(model.Options.HighPassFilterCutoffHz) msg = sprintf(['No high-pass filter for Node "%s"\n', ... @@ -106,7 +107,7 @@ opt.verbosity = 0; if obj.verbose - opt.verbosity = 1; + opt.verbosity = 3; end logger('INFO', msg, 'filename', mfilename(), 'options', opt); diff --git a/src/bids_model/checkGroupBy.m b/src/bids_model/checkGroupBy.m index 9482a16e3..765f6496c 100644 --- a/src/bids_model/checkGroupBy.m +++ b/src/bids_model/checkGroupBy.m @@ -61,7 +61,7 @@ status = true; elseif numel(groupBy) == 2 && iscellstr(extraVar) && numel(extraVar) > 0 for i = 1:numel(extraVar) - if all(ismember(lower(groupBy), {'contrast', extraVar{i}})) + if all(ismember(groupBy, {'contrast', extraVar{i}})) status = true; break end diff --git a/tests/create_3_groups_dataset.py b/tests/create_3_groups_dataset.py index ea52607de..c42274851 100644 --- a/tests/create_3_groups_dataset.py +++ b/tests/create_3_groups_dataset.py @@ -15,7 +15,7 @@ ROOT_DIR = Path(__file__).parent.parent START_DIR = Path(__file__).parent SUBJECTS = ["ctrl01", "blind01", "relative01", "ctrl02", "blind02", "relative02"] -SUBJECTS = ["ctrl01", "blind01", "relative01"] +# SUBJECTS = ["ctrl01", "blind01", "relative01"] SESSIONS = ["01", "02"] SESSIONS = ["01"] diff --git a/tests/data/models/model-bug616_smdl.json b/tests/data/models/model-bug616_smdl.json index 84dba40ac..46817dba1 100644 --- a/tests/data/models/model-bug616_smdl.json +++ b/tests/data/models/model-bug616_smdl.json @@ -209,7 +209,7 @@ "Type": "glm", "X": [ 1, - "group" + "Group" ] }, "Contrasts": [ diff --git a/tests/data/models/model-vislocalizer2sampleTTest_smdl.json b/tests/data/models/model-vislocalizer2sampleTTest_smdl.json index 5fd294226..b87b3ef8d 100644 --- a/tests/data/models/model-vislocalizer2sampleTTest_smdl.json +++ b/tests/data/models/model-vislocalizer2sampleTTest_smdl.json @@ -86,7 +86,7 @@ "Type": "glm", "X": [ 1, - "group" + "Group" ], "Options": { "Mask": { diff --git a/tests/data/models/model-vismotionNoOverWrite_smdl.json b/tests/data/models/model-vismotionNoOverWrite_smdl.json index 1aa9a9ff6..6ed049f33 100644 --- a/tests/data/models/model-vismotionNoOverWrite_smdl.json +++ b/tests/data/models/model-vismotionNoOverWrite_smdl.json @@ -199,7 +199,7 @@ "Type": "glm", "X": [ 1, - "group" + "Group" ], "Options": { "Mask": { diff --git a/tests/data/3_groups/models/model-vismotion_smdl.json b/tests/data/models/model-vismotionOneWayANOVA_smdl.json similarity index 100% rename from tests/data/3_groups/models/model-vismotion_smdl.json rename to tests/data/models/model-vismotionOneWayANOVA_smdl.json diff --git a/tests/tests_slow/tests_batches/stats/test_setBatchFactorialDesign_3_groups.m b/tests/tests_slow/tests_batches/stats/test_setBatchFactorialDesign_3_groups.m index 2220c850e..cae5cae0a 100644 --- a/tests/tests_slow/tests_batches/stats/test_setBatchFactorialDesign_3_groups.m +++ b/tests/tests_slow/tests_batches/stats/test_setBatchFactorialDesign_3_groups.m @@ -14,7 +14,9 @@ function test_setBatchFactorialDesign_between_3_groups() opt.verbosity = 3; matlabbatch = {}; - matlabbatch = setBatchFactorialDesign(matlabbatch, opt, 'between_groups'); + [matlabbatch, contrastsList, groups] = setBatchFactorialDesign(matlabbatch, ... + opt, ... + 'between_groups'); matlabbatch; end diff --git a/tests/utils/setOptions.m b/tests/utils/setOptions.m index 2259ab216..543f98a01 100644 --- a/tests/utils/setOptions.m +++ b/tests/utils/setOptions.m @@ -120,7 +120,7 @@ opt.dir.preproc = fullfile(dataDir, 'derivatives', 'bidspm-preproc'); opt.dir.stats = fullfile(dataDir, 'derivatives', 'bidspm-stats'); - opt.model.file = fullfile(dataDir, 'models', 'model-vismotion_smdl.json'); + opt.model.file = fullfile(getTestDataDir(), 'models', 'model-vismotionOneWayANOVA_smdl.json'); else