Skip to content

Commit

Permalink
update code and tests
Browse files Browse the repository at this point in the history
  • Loading branch information
Remi-Gau committed Jul 25, 2024
1 parent e1ee976 commit ab138de
Show file tree
Hide file tree
Showing 10 changed files with 138 additions and 23 deletions.
140 changes: 126 additions & 14 deletions src/batches/stats/setBatchFactorialDesign.m
Original file line number Diff line number Diff line change
Expand Up @@ -36,33 +36,107 @@
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);

Check warning on line 49 in src/batches/stats/setBatchFactorialDesign.m

View check run for this annotation

Codecov / codecov/patch

src/batches/stats/setBatchFactorialDesign.m#L48-L49

Added lines #L48 - L49 were not covered by tests
elseif numel(groupBy) == 2
[matlabbatch, contrastsList, groups] = tTestForGroup(matlabbatch, opt, nodeName);

Check warning on line 51 in src/batches/stats/setBatchFactorialDesign.m

View check run for this annotation

Codecov / codecov/patch

src/batches/stats/setBatchFactorialDesign.m#L51

Added line #L51 was not covered by tests
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
% the later is needed to be able to create proper RFX dir name
contrastsList = {};
groups = {};

Check warning on line 136 in src/batches/stats/setBatchFactorialDesign.m

View check run for this annotation

Codecov / codecov/patch

src/batches/stats/setBatchFactorialDesign.m#L135-L136

Added lines #L135 - L136 were not covered by tests

contrasts = getContrastsListForFactorialDesign(opt, nodeName);

Check warning on line 138 in src/batches/stats/setBatchFactorialDesign.m

View check run for this annotation

Codecov / codecov/patch

src/batches/stats/setBatchFactorialDesign.m#L138

Added line #L138 was not covered by tests

% collect all con images from all subjects
for iSub = 1:numel(opt.subjects)
subLabel = opt.subjects{iSub};
Expand All @@ -87,15 +161,23 @@
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
% the later is needed to be able to create proper RFX dir name
contrastsList = {};
groups = {};

Check warning on line 170 in src/batches/stats/setBatchFactorialDesign.m

View check run for this annotation

Codecov / codecov/patch

src/batches/stats/setBatchFactorialDesign.m#L169-L170

Added lines #L169 - L170 were not covered by tests

[BIDS, opt] = getData(opt, opt.dir.preproc);

Check warning on line 172 in src/batches/stats/setBatchFactorialDesign.m

View check run for this annotation

Codecov / codecov/patch

src/batches/stats/setBatchFactorialDesign.m#L172

Added line #L172 was not covered by tests

contrasts = getContrastsListForFactorialDesign(opt, nodeName);

Check warning on line 174 in src/batches/stats/setBatchFactorialDesign.m

View check run for this annotation

Codecov / codecov/patch

src/batches/stats/setBatchFactorialDesign.m#L174

Added line #L174 was not covered by tests

node = opt.model.bm.get_nodes('Name', nodeName);
groupBy = node.GroupBy;

Check warning on line 177 in src/batches/stats/setBatchFactorialDesign.m

View check run for this annotation

Codecov / codecov/patch

src/batches/stats/setBatchFactorialDesign.m#L176-L177

Added lines #L176 - L177 were not covered by tests

groupColumnHdr = setxor(groupBy, {'contrast'});
groupColumnHdr = groupColumnHdr{1};

Check warning on line 180 in src/batches/stats/setBatchFactorialDesign.m

View check run for this annotation

Codecov / codecov/patch

src/batches/stats/setBatchFactorialDesign.m#L179-L180

Added lines #L179 - L180 were not covered by tests

checkColumnParticipantsTsv(BIDS, groupColumnHdr);

Check warning on line 182 in src/batches/stats/setBatchFactorialDesign.m

View check run for this annotation

Codecov / codecov/patch

src/batches/stats/setBatchFactorialDesign.m#L182

Added line #L182 was not covered by tests

Expand Down Expand Up @@ -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);

Expand Down
5 changes: 3 additions & 2 deletions src/bids_model/BidsModel.m
Original file line number Diff line number Diff line change
Expand Up @@ -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', ...
Expand All @@ -106,7 +107,7 @@

opt.verbosity = 0;
if obj.verbose
opt.verbosity = 1;
opt.verbosity = 3;

Check warning on line 110 in src/bids_model/BidsModel.m

View check run for this annotation

Codecov / codecov/patch

src/bids_model/BidsModel.m#L110

Added line #L110 was not covered by tests
end
logger('INFO', msg, 'filename', mfilename(), 'options', opt);

Expand Down
2 changes: 1 addition & 1 deletion src/bids_model/checkGroupBy.m
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion tests/create_3_groups_dataset.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"]

Expand Down
2 changes: 1 addition & 1 deletion tests/data/models/model-bug616_smdl.json
Original file line number Diff line number Diff line change
Expand Up @@ -209,7 +209,7 @@
"Type": "glm",
"X": [
1,
"group"
"Group"
]
},
"Contrasts": [
Expand Down
2 changes: 1 addition & 1 deletion tests/data/models/model-vislocalizer2sampleTTest_smdl.json
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@
"Type": "glm",
"X": [
1,
"group"
"Group"
],
"Options": {
"Mask": {
Expand Down
2 changes: 1 addition & 1 deletion tests/data/models/model-vismotionNoOverWrite_smdl.json
Original file line number Diff line number Diff line change
Expand Up @@ -199,7 +199,7 @@
"Type": "glm",
"X": [
1,
"group"
"Group"
],
"Options": {
"Mask": {
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
2 changes: 1 addition & 1 deletion tests/utils/setOptions.m
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down

0 comments on commit ab138de

Please sign in to comment.