From 64b14da52be73373dca1187b1016ade9032f41d4 Mon Sep 17 00:00:00 2001 From: Remi Gau Date: Tue, 30 Jul 2024 15:37:44 +0200 Subject: [PATCH] [ENH] implement one-way ANOVA at group level to compare groups (#1296) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * start adding a 3 group dataset * [ENH] use sub commands for python CLI (#1292) * start adding sub commands * preproc subcommand * deal with create roi parser * deal with default_model sub command * deal with bms action * deal with stats actions * start adapting cli * add test for command building * add test for CLI for preproc, smooth... * add test for command building for stats * add test for command building for bms * rename function * start switching CLI * several fixes * more fix * keep switchingµ * fix * fix * linti * fixes * fixes * minor fixes * additional fixes * fix output ROI * visualize output in circle ci * add temp script for boutiques and bids model graph * start adding tweaks * try 3 groups * rename file * update code and tests * add multi group demo * fixes * adapt example * fix contrasts * deal with results * update doc * minor changes * tmp * refactor * refactor * refactor * several fixes * several fixes * fix and refactor * fix and refactor * fix and refactor * fix models * fix * fix * fix fast tests * fixes * fix and implement F test * run group level node by node * add F test to demo * octave fix --- .flake8 | 1 + .github/workflows/run_tests_cli.yml | 14 +- CHANGELOG.md | 3 +- WIP/find_data_set_with_3_groups.py | 43 +++ demos/openneuro/Makefile | 10 + demos/openneuro/ds000114_run.m | 2 + demos/openneuro/ds003397_run.m | 61 ++++ .../openneuro/models/model-ds003397_smdl.json | 218 ++++++++++++ docs/source/API.rst | 17 +- docs/source/demos/openneuro.md | 16 + docs/source/stats/bids_stats_model.md | 25 +- .../examples/model_oneWayAnova_smdl.json | 125 +++++++ lib/bids-matlab | 2 +- miss_hit.cfg | 4 +- .../setBatchLesionAbnormalitiesDetection.m | 2 - src/batches/stats/setBatchFactorialDesign.m | 314 +++++++++++++----- ...etBatchFactorialDesignGlobalCalcAndNorm.m} | 4 +- .../stats/setBatchGroupLevelContrasts.m | 108 ++++-- .../stats/setBatchSubjectLevelContrasts.m | 13 +- .../stats/setBatchSubjectLevelGLMSpec.m | 12 +- src/batches/stats/setBatchTwoSampleTTest.m | 210 ------------ src/bids/checkColumnParticipantsTsv.m | 21 -- src/bids/getAvailableGroups.m | 25 ++ src/bids_model/BidsModel.m | 192 ++++++++++- src/bids_model/checkContrast.m | 3 +- src/bids_model/checkGroupBy.m | 71 ---- src/bids_model/getContrastsList.m | 8 +- .../getContrastsListForFactorialDesign.m | 27 +- src/bids_model/getDummyContrastsList.m | 12 +- src/bids_model/getInclusiveMask.m | 8 +- src/messages/noSPMmat.m | 45 --- src/stats/group_level/getRFXdir.m | 7 +- src/stats/group_level/groupLevelGlmType.m | 25 -- .../subject_level/convertOnsetTsvToMat.m | 8 +- .../createAndReturnCounfoundMatFile.m | 6 +- src/stats/subject_level/specifyContrasts.m | 34 +- .../subject_level/specifyDummyContrasts.m | 7 +- src/stats/utils/checkSpmMat.m | 19 ++ src/stats/utils/getContrastNb.m | 2 + .../removeIntercept.m | 0 src/workflows/bidsChangeSuffix.m | 4 +- src/workflows/bidsRename.m | 1 - src/workflows/lesion/bidsLesionOverlapMap.m | 2 - .../preproc/bidsWholeBrainFuncMask.m | 2 - src/workflows/roi/bidsRoiBasedGLM.m | 10 +- src/workflows/stats/bidsConcatBetaTmaps.m | 2 +- src/workflows/stats/bidsFFX.m | 30 +- src/workflows/stats/bidsModelSelection.m | 9 +- src/workflows/stats/bidsRFX.m | 27 +- src/workflows/stats/bidsResults.m | 67 ++-- tests/.gitignore | 3 + tests/Makefile | 1 + tests/create_3_groups_dataset.py | 89 +++++ tests/create_dummy_dataset.py | 113 +------ .../data/bidspm-raw/dataset_description.json | 2 +- tests/data/models/model-bug616_smdl.json | 2 +- .../model-vislocalizer2sampleTTest_smdl.json | 11 +- ...vislocalizerSeveralDatasetLevels_smdl.json | 14 + .../model-vismotionNoOverWrite_smdl.json | 34 +- .../model-vismotionOneWayANOVA_smdl.json | 202 +++++++++++ ...del-vismotionSeveralDatasetLevel_smdl.json | 14 + .../preproc/test_setBatchGenerateT1map.m | 2 +- .../stats/test_setBatchGroupLevelContrasts.m | 4 +- tests/tests_bids/test_fileFilterForBold.m | 14 +- .../test_checkColumnParticipantsTsv.m | 23 -- tests/tests_bids_model/test_checkGroupBy.m | 57 ---- .../test_getOptionsFromModel.m | 8 +- .../tests_bids_model/test_groupLevelGlmType.m | 47 +++ tests/tests_bids_model/test_validateGroupBy.m | 76 +++++ .../stats/test_setBatchFactorialDesign.m | 6 + .../test_setBatchFactorialDesign_3_groups.m | 23 ++ .../tests_workflows/stats/test_bidsRFX.m | 128 ------- .../stats/test_bidsRFX_2_groups.m | 121 +++++++ .../stats/test_bidsRFX_3_groups.m | 93 ++++++ .../stats/test_bidsRFX_chain_several_nodes.m | 77 +++++ .../tests_workflows/stats/test_bidsResults.m | 44 --- .../test_renameSegmentParameter.m | 2 +- .../test_renameUnwarpParameter.m | 2 +- .../subject_level/test_specifyContrasts.m | 159 +++++---- .../test_specifyDummyContrasts_Session.m | 3 + .../unit_tests/test_getInclusiveMask.m | 4 +- tests/tests_unit/test_warnings.m | 25 -- .../stats/test_bidsRFX_3_groups.m | 64 ++++ .../tests_workflows/stats/test_bidsResults.m | 88 +++++ tests/utils.py | 107 ++++++ tests/utils/createDummyData.m | 1 + tests/utils/setOptions.m | 12 + 87 files changed, 2417 insertions(+), 1136 deletions(-) create mode 100644 WIP/find_data_set_with_3_groups.py create mode 100644 demos/openneuro/ds003397_run.m create mode 100644 demos/openneuro/models/model-ds003397_smdl.json create mode 100644 docs/source/stats/examples/model_oneWayAnova_smdl.json rename src/batches/stats/{setBatchFatorialDesignGlobalCalcAndNorm.m => setBatchFactorialDesignGlobalCalcAndNorm.m} (54%) delete mode 100644 src/batches/stats/setBatchTwoSampleTTest.m delete mode 100644 src/bids/checkColumnParticipantsTsv.m create mode 100644 src/bids/getAvailableGroups.m delete mode 100644 src/bids_model/checkGroupBy.m delete mode 100644 src/messages/noSPMmat.m delete mode 100644 src/stats/group_level/groupLevelGlmType.m create mode 100644 src/stats/utils/checkSpmMat.m rename src/stats/{subject_level => utils}/removeIntercept.m (100%) create mode 100644 tests/create_3_groups_dataset.py create mode 100644 tests/data/models/model-vismotionOneWayANOVA_smdl.json delete mode 100644 tests/tests_bids/unit_tests/test_checkColumnParticipantsTsv.m delete mode 100644 tests/tests_bids_model/test_checkGroupBy.m create mode 100644 tests/tests_bids_model/test_groupLevelGlmType.m create mode 100644 tests/tests_bids_model/test_validateGroupBy.m create mode 100644 tests/tests_slow/tests_batches/stats/test_setBatchFactorialDesign_3_groups.m create mode 100644 tests/tests_slow/tests_workflows/stats/test_bidsRFX_2_groups.m create mode 100644 tests/tests_slow/tests_workflows/stats/test_bidsRFX_3_groups.m create mode 100644 tests/tests_slow/tests_workflows/stats/test_bidsRFX_chain_several_nodes.m create mode 100644 tests/tests_workflows/stats/test_bidsRFX_3_groups.m create mode 100644 tests/tests_workflows/stats/test_bidsResults.m create mode 100644 tests/utils.py diff --git a/.flake8 b/.flake8 index 9a0a812a1..70be08b06 100644 --- a/.flake8 +++ b/.flake8 @@ -9,6 +9,7 @@ exclude = tests/* _version.py demos/* + WIP/* count = True show-source = True statistics = True diff --git a/.github/workflows/run_tests_cli.yml b/.github/workflows/run_tests_cli.yml index a1a4a3988..eb3b15512 100644 --- a/.github/workflows/run_tests_cli.yml +++ b/.github/workflows/run_tests_cli.yml @@ -60,10 +60,10 @@ jobs: coverage run --source src -m pytest coverage xml - - name: Code coverage - uses: codecov/codecov-action@v4 - with: - file: coverage.xml - flags: cli - name: codecov-cli - fail_ci_if_error: false + # - name: Code coverage + # uses: codecov/codecov-action@v4 + # with: + # file: coverage.xml + # flags: cli + # name: codecov-cli + # fail_ci_if_error: false diff --git a/CHANGELOG.md b/CHANGELOG.md index 84ca5bd38..3e526c589 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +* [ENH] add support for one-way ANOVA across groups at the group level #1296 by @Remi-Gau +* [ENH] allow for 2 sample T-Test, within group T-Test and one-way ANOVA to ne more flexible with respect to what praticipants.tsv column to use to allocate subjects in each group #1296 by @Remi-Gau * [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 @@ -43,7 +45,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed - * [ENH] align specification of F contrasts on the BIDS stats model: they should now be specified as a 2D matrix and not a 1D vector. #1276 @Remi-Gau * [DOC] change theme and structure of the documentation #1256 @Remi-Gau * [REF] Refactor and update CLI in #1096 @Remi-Gau diff --git a/WIP/find_data_set_with_3_groups.py b/WIP/find_data_set_with_3_groups.py new file mode 100644 index 000000000..af91f393d --- /dev/null +++ b/WIP/find_data_set_with_3_groups.py @@ -0,0 +1,43 @@ +from pathlib import Path + +import pandas as pd +from pandas.errors import ParserError + +openneuro = Path("/home/remi/datalad/datasets.datalad.org/openneuro") + +for ds in openneuro.iterdir(): + file = ds / "participants.tsv" + + if file.exists(): + + try: + df = pd.read_csv(file, sep="\t") + except ParserError: + ... + + if len(df) < 10: + continue + if "group" in df.columns: + col = "group" + elif "Group" in df.columns: + col = "Group" + else: + continue + if len(df[col].value_counts()) < 3: + continue + if "age" in df.columns and df["age"].mean() < 18: + continue + bold_files = list(ds.glob("**/*bold*")) + + has_func = len(bold_files) > 0 + if not has_func: + continue + + tasks = {x.name.split("task-")[1].split("_")[0] for x in bold_files} + if len(tasks) == 1 and "rest" in tasks: + continue + + print() + print(file) + print(df[col].value_counts()) + print(tasks) diff --git a/demos/openneuro/Makefile b/demos/openneuro/Makefile index c563da7e8..f714be174 100644 --- a/demos/openneuro/Makefile +++ b/demos/openneuro/Makefile @@ -74,3 +74,13 @@ data_ds002799: datalad get -d inputs/ds002799 inputs/ds002799/derivatives/fmriprep/sub-292/*/func/*tsv datalad get -d inputs/ds002799 inputs/ds002799/derivatives/fmriprep/sub-30[27]/*/func/*MNI152NLin2009cAsym* -J 2 datalad get -d inputs/ds002799 inputs/ds002799/derivatives/fmriprep/sub-30[27]/*/func/*tsv -J 2 + +data_ds003397: + mkdir -p inputs + cd inputs && datalad install ///openneuro/ds003397 + cd inputs && datalad install https://github.com/OpenNeuroDerivatives/ds003397-fmriprep + cd inputs/ds003397-fmriprep && datalad get sub-[01][1267]/anat/*MNI152NLin2009cAsym*T1w.nii.gz -J 12 + cd inputs/ds003397-fmriprep && datalad get sub-[01][1267]/anat/*mask*.nii.gz -J 12 + cd inputs/ds003397-fmriprep && datalad get sub-[01][1267]/func/*time*tsv -J 12 + cd inputs/ds003397-fmriprep && datalad get sub-[01][1267]/func/*json -J 12 + cd inputs/ds003397-fmriprep && datalad get sub-[01][1267]/func/*MNI152NLin2009cAsym*desc-preproc*bold.nii.gz -J 12 diff --git a/demos/openneuro/ds000114_run.m b/demos/openneuro/ds000114_run.m index a5533b5e0..6f9df24b0 100644 --- a/demos/openneuro/ds000114_run.m +++ b/demos/openneuro/ds000114_run.m @@ -1,3 +1,5 @@ +% Demo to compare activations across sessions. + % (C) Copyright 2023 bidspm developers clear; diff --git a/demos/openneuro/ds003397_run.m b/demos/openneuro/ds003397_run.m new file mode 100644 index 000000000..aa793a88f --- /dev/null +++ b/demos/openneuro/ds003397_run.m @@ -0,0 +1,61 @@ +% Run a one-way ANOVA across group +% +% Only a few subjects are run because of the large size of each run. + +% (C) Copyright 2024 bidspm developers + +clear; +clc; + +addpath(fullfile(pwd, '..', '..')); +bidspm(); + +task = 'checkerboard'; + +% The directory where the data are located +root_dir = fileparts(mfilename('fullpath')); +bids_dir = fullfile(root_dir, 'inputs', 'ds003397'); +fmriprep_dir = fullfile(root_dir, 'inputs', 'ds003397-fmriprep'); +output_dir = fullfile(root_dir, 'outputs', 'ds003397', 'derivatives'); + +space = {'MNI152NLin2009cAsym'}; +participant_label = {'01', '06', '11', '12'}; + +% Copy (& smooth) +FWHM = 0; +bidspm(fmriprep_dir, output_dir, 'subject', ... + 'participant_label', participant_label, ... + 'action', 'smooth', ... + 'task', task, ... + 'space', space, ... + 'fwhm', FWHM, ... + 'verbosity', 3); + +%% Stats +preproc_dir = fullfile(output_dir, 'bidspm-preproc'); + +model_file = fullfile(root_dir, ... + 'models', ... + 'model-ds003397_smdl.json'); + +bidspm(bids_dir, output_dir, 'subject', ... + 'participant_label', participant_label, ... + 'action', 'stats', ... + 'preproc_dir', preproc_dir, ... + 'model_file', model_file, ... + 'roi_atlas', 'hcpex', ... + 'space', space, ... + 'fwhm', 0, ... + 'skip_validation', true, ... + 'verbosity', 3); + +bidspm(bids_dir, output_dir, 'dataset', ... + 'participant_label', participant_label, ... + 'action', 'stats', ... + 'preproc_dir', preproc_dir, ... + 'model_file', model_file, ... + 'roi_atlas', 'hcpex', ... + 'space', space, ... + 'fwhm', 0, ... + 'skip_validation', true, ... + 'verbosity', 3); diff --git a/demos/openneuro/models/model-ds003397_smdl.json b/demos/openneuro/models/model-ds003397_smdl.json new file mode 100644 index 000000000..cb5ef1606 --- /dev/null +++ b/demos/openneuro/models/model-ds003397_smdl.json @@ -0,0 +1,218 @@ +{ + "Name": "1_way_ANOVA", + "BIDSModelVersion": "1.0.0", + "Input": { + "task": [ + "checkerboard" + ], + "space": [ + "MNI152NLin2009cAsym" + ] + }, + "Nodes": [ + { + "Level": "Run", + "Name": "run_level", + "GroupBy": [ + "run", + "subject" + ], + "Model": { + "Type": "glm", + "X": [ + "trial_type.flashing checkerboard", + "trans_?", + "rot_?" + ], + "HRF": { + "Variables": [ + "trial_type.flashing checkerboard" + ], + "Model": "spm" + } + }, + "Contrasts": [ + { + "Name": "flashing checkerboard", + "ConditionList": [ + "trial_type.flashing checkerboard" + ], + "Weights": [ + 1 + ], + "Test": "t" + } + ] + }, + { + "Level": "Subject", + "Name": "subject_level", + "GroupBy": [ + "contrast", + "subject" + ], + "Model": { + "Type": "glm", + "X": [ + 1 + ], + "Software": { + "bidspm": { + "Results": [ + { + "name": [ + "flashing checkerboard" + ], + "p": 0.05, + "MC": "FWE", + "png": true, + "binary": false, + "nidm": false, + "montage": { + "do": true, + "slices": [ + -4, + 0, + 4, + 8, + 16 + ], + "background": { + "suffix": "T1w", + "desc": "preproc", + "modality": "anat" + } + } + } + ] + } + } + }, + "DummyContrasts": { + "Test": "t" + } + }, + { + "Level": "Dataset", + "Name": "dataset_level", + "Description": "average across all subjects", + "GroupBy": [ + "contrast" + ], + "Model": { + "Type": "glm", + "X": [ + 1 + ] + } + }, + { + "Level": "Dataset", + "Name": "between_groups", + "Description": "one way anova", + "GroupBy": [ + "contrast" + ], + "Model": { + "Type": "glm", + "X": [ + 1, + "group" + ], + "Software": { + "bidspm": { + "Results": [ + { + "name": [ + "B > I", + "average across groups" + ], + "p": 0.01, + "MC": "FWE", + "png": true, + "binary": false, + "nidm": false, + "montage": { + "do": false + } + } + ] + } + } + }, + "Contrasts": [ + { + "Name": "B > I", + "ConditionList": [ + "group.B", + "group.I" + ], + "Weights": [ + 1, + -1 + ], + "Test": "t" + }, + { + "Name": "average across groups", + "ConditionList": [ + "group.B", + "group.I", + "group.BI" + ], + "Weights": [ + 1, + 1, + 1 + ], + "Test": "t" + }, + { + "Name": "some F test", + "ConditionList": [ + "group.B", + "group.BI", + "group.I" + ], + "Weights": [ + [ + 1, + 0, + 0 + ], + [ + 0, + 1, + 0 + ], + [ + 0, + 0, + 1 + ] + ], + "Test": "F" + } + ] + } + ], + "Edges": [ + { + "Source": "run_level", + "Destination": "subject_level" + }, + { + "Source": "subject_level", + "Destination": "dataset_level" + }, + { + "Source": "subject_level", + "Destination": "between_groups", + "Filter": { + "contrast": [ + "flashing checkerboard" + ] + } + } + ] +} diff --git a/docs/source/API.rst b/docs/source/API.rst index 8e189bfa1..f03955971 100644 --- a/docs/source/API.rst +++ b/docs/source/API.rst @@ -93,8 +93,8 @@ stats .. autofunction:: src.batches.stats.setBatchContrasts .. autofunction:: src.batches.stats.setBatchEstimateModel .. autofunction:: src.batches.stats.setBatchFactorialDesign +.. autofunction:: src.batches.stats.setBatchFactorialDesignGlobalCalcAndNorm .. autofunction:: src.batches.stats.setBatchFactorialDesignImplicitMasking -.. autofunction:: src.batches.stats.setBatchFatorialDesignGlobalCalcAndNorm .. autofunction:: src.batches.stats.setBatchGoodnessOfFit .. autofunction:: src.batches.stats.setBatchGroupLevelContrasts .. autofunction:: src.batches.stats.setBatchGroupLevelResults @@ -167,13 +167,13 @@ bids ==== .. autofunction:: src.bids.addStcToQuery .. autofunction:: src.bids.buildIndividualSpaceRoiFilename -.. autofunction:: src.bids.checkColumnParticipantsTsv .. autofunction:: src.bids.checkFmriprep .. autofunction:: src.bids.fileFilterForBold .. autofunction:: src.bids.generatedBy .. autofunction:: src.bids.getAnatFilename .. autofunction:: src.bids.getAndCheckRepetitionTime .. autofunction:: src.bids.getAndCheckSliceOrder +.. autofunction:: src.bids.getAvailableGroups .. autofunction:: src.bids.getBoldFilename .. autofunction:: src.bids.getInfo .. autofunction:: src.bids.getMeanFuncFilename @@ -189,7 +189,6 @@ bids bids_model ========== .. autofunction:: src.bids_model.checkContrast -.. autofunction:: src.bids_model.checkGroupBy .. autofunction:: src.bids_model.createDefaultStatsModel .. autofunction:: src.bids_model.createModelFamilies .. autofunction:: src.bids_model.getContrastsFromParentNode @@ -199,6 +198,12 @@ bids_model .. autofunction:: src.bids_model.getDummyContrastsList .. autofunction:: src.bids_model.getInclusiveMask +bidspm +====== + +data +---- + cli === .. autofunction:: src.cli.baseInputParser @@ -221,6 +226,7 @@ cli constants ========= +.. autofunction:: src.constants.bidsAppsActions .. autofunction:: src.constants.lowLevelActions defaults @@ -263,7 +269,6 @@ messages .. autofunction:: src.messages.errorHandling .. autofunction:: src.messages.logger .. autofunction:: src.messages.noRoiFound -.. autofunction:: src.messages.noSPMmat .. autofunction:: src.messages.notImplemented .. autofunction:: src.messages.printAvailableContrasts .. autofunction:: src.messages.printBatchName @@ -340,7 +345,6 @@ subject_level .. autofunction:: src.stats.subject_level.getSessionForRegressorNb .. autofunction:: src.stats.subject_level.newContrast .. autofunction:: src.stats.subject_level.orderAndPadCounfoundMatFile -.. autofunction:: src.stats.subject_level.removeIntercept .. autofunction:: src.stats.subject_level.reorderCounfounds .. autofunction:: src.stats.subject_level.sanitizeConfounds .. autofunction:: src.stats.subject_level.saveRoiGlmSummaryTable @@ -355,10 +359,10 @@ group_level .. autofunction:: src.stats.group_level.computeCumulativeFwhm .. autofunction:: src.stats.group_level.findSubjectConImage .. autofunction:: src.stats.group_level.getRFXdir -.. autofunction:: src.stats.group_level.groupLevelGlmType utils ----- +.. autofunction:: src.stats.utils.checkSpmMat .. autofunction:: src.stats.utils.createGlmDirName .. autofunction:: src.stats.utils.designMatrixFigureName .. autofunction:: src.stats.utils.endsWithRunNumber @@ -367,6 +371,7 @@ utils .. autofunction:: src.stats.utils.getRegressorIdx .. autofunction:: src.stats.utils.labelActivations .. autofunction:: src.stats.utils.labelSpmSessWithBidsSesAndRun +.. autofunction:: src.stats.utils.removeIntercept .. autofunction:: src.stats.utils.returnContrastImageFile .. autofunction:: src.stats.utils.returnNumberScrubbedTimePoints .. autofunction:: src.stats.utils.validateContrasts diff --git a/docs/source/demos/openneuro.md b/docs/source/demos/openneuro.md index b1b704715..00026f14e 100644 --- a/docs/source/demos/openneuro.md +++ b/docs/source/demos/openneuro.md @@ -144,3 +144,19 @@ transformers cannot yet be appled to confounds - resting state and task - several sessions - fmriprep data + +## ds003379 + +- [Source dataset](https://openneuro.org/datasets/ds003379) + +### Features + +- checkerboard localizer +- 3 groups +- fmriprep data + +### Scripts + +```{eval-rst} +.. autoscript:: ds003379_run +``` diff --git a/docs/source/stats/bids_stats_model.md b/docs/source/stats/bids_stats_model.md index dc7e631af..d3a610c7a 100644 --- a/docs/source/stats/bids_stats_model.md +++ b/docs/source/stats/bids_stats_model.md @@ -487,30 +487,43 @@ Subject level contrast averaging beta across runs At the moment only, the only type of models that are supported are: -- one sample t-test: averaging across all subjects +### one sample t-test: averaging across all subjects ```{literalinclude} ./examples/model-datasetLevel_smdl.json :language: json ``` -- one sample t-test: averaging across all subjects of a specific group +### one sample t-test: averaging across all subjects of a specific group ```{literalinclude} ./examples/model_withinGroup_smdl.json :language: json :emphasize-lines: 5-8 ``` -- 2 samples t-test: comparing 2 groups +### 2 samples t-test: comparing 2 groups -At the moment this can only be based on how participants are allocated to a -group based on a `group` or `Group` column in the `participants.tsv` of in the -raw dataset. +Participants are allocated to a group based on a the corresponding column +in the `participants.tsv` of in the raw dataset (`Group` in this case). ```{literalinclude} ./examples/model_betweenGroups_smdl.json :language: json :emphasize-lines: 8-28 ``` +### one way ANOVA: comparing several groups + +Participants are allocated to a group based on a the corresponding column +in the `participants.tsv` of in the raw dataset (`Group` in this case). + +```{literalinclude} ./examples/model_oneWayAnova_smdl.json +:language: json +:emphasize-lines: 64-108 +``` + +```{note} +Contrast must be explicitly passed the `Edges.Filter.contrast`. +``` + ### Method section It is possible to write a draft of method section based on a BIDS statistical model. diff --git a/docs/source/stats/examples/model_oneWayAnova_smdl.json b/docs/source/stats/examples/model_oneWayAnova_smdl.json new file mode 100644 index 000000000..afc164f93 --- /dev/null +++ b/docs/source/stats/examples/model_oneWayAnova_smdl.json @@ -0,0 +1,125 @@ +{ + "Name": "1_way_ANOVA", + "BIDSModelVersion": "1.0.0", + "Input": { + "task": [ + "checkerboard" + ], + "space": [ + "MNI152NLin2009cAsym" + ] + }, + "Nodes": [ + { + "Level": "Run", + "Name": "run_level", + "GroupBy": [ + "run", + "subject" + ], + "Model": { + "Type": "glm", + "X": [ + "trial_type.flashing checkerboard", + "trans_?", + "rot_?" + ], + "HRF": { + "Variables": [ + "trial_type.flashing checkerboard" + ], + "Model": "spm" + } + }, + "Contrasts": [ + { + "Name": "flashing checkerboard", + "ConditionList": [ + "trial_type.flashing checkerboard" + ], + "Weights": [ + 1 + ], + "Test": "t" + } + ] + }, + { + "Level": "Subject", + "Name": "subject_level", + "GroupBy": [ + "contrast", + "subject" + ], + "Model": { + "Type": "glm", + "X": [ + 1 + ] + }, + "DummyContrasts": { + "Test": "t" + } + }, + { + "Level": "Dataset", + "Name": "between_groups", + "Description": "one way anova", + "GroupBy": [ + "contrast" + ], + "Model": { + "Type": "glm", + "X": [ + 1, + "group" + ] + }, + "Contrasts": [ + { + "Name": "B > I", + "ConditionList": [ + "group.B", + "group.I", + "group.BI" + ], + "Weights": [ + 1, + -1, + 0 + ], + "Test": "t" + }, + { + "Name": "average across groups", + "ConditionList": [ + "group.B", + "group.I", + "group.BI" + ], + "Weights": [ + 1, + 1, + 1 + ], + "Test": "t" + } + ] + } + ], + "Edges": [ + { + "Source": "run_level", + "Destination": "subject_level" + }, + { + "Source": "subject_level", + "Destination": "between_groups", + "Filter": { + "contrast": [ + "flashing checkerboard" + ] + } + } + ] +} diff --git a/lib/bids-matlab b/lib/bids-matlab index 5833faeb4..c575a111e 160000 --- a/lib/bids-matlab +++ b/lib/bids-matlab @@ -1 +1 @@ -Subproject commit 5833faeb4f236cccd472ad2281f2b97bb3710dd6 +Subproject commit c575a111edfd5fd188f587099990fe9fb0f66866 diff --git a/miss_hit.cfg b/miss_hit.cfg index 706fee817..cd8e82beb 100644 --- a/miss_hit.cfg +++ b/miss_hit.cfg @@ -37,6 +37,6 @@ tab_width: 2 # metrics limit for the code quality (https://florianschanda.github.io/miss_hit/metrics.html) metric "cnest": limit 5 -metric "file_length": limit 750 -metric "cyc": limit 19 +metric "file_length": limit 1000 +metric "cyc": limit 22 metric "parameters": limit 7 diff --git a/src/batches/lesion/setBatchLesionAbnormalitiesDetection.m b/src/batches/lesion/setBatchLesionAbnormalitiesDetection.m index 3b7efda7c..feea26789 100755 --- a/src/batches/lesion/setBatchLesionAbnormalitiesDetection.m +++ b/src/batches/lesion/setBatchLesionAbnormalitiesDetection.m @@ -16,8 +16,6 @@ % (C) Copyright 2021 bidspm developers - % TODO add test - printBatchName('Lesion abnormalities', opt); outliers_detection = opt.toolbox.ALI.outliers_detection; diff --git a/src/batches/stats/setBatchFactorialDesign.m b/src/batches/stats/setBatchFactorialDesign.m index cc24ca476..090d1cc3a 100644 --- a/src/batches/stats/setBatchFactorialDesign.m +++ b/src/batches/stats/setBatchFactorialDesign.m @@ -21,107 +21,209 @@ % (C) Copyright 2019 bidspm developers - % TODO implement other models than group average of contrast from lower levels - - % TODO implement Contrasts and not just dummy contrasts - % 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 = {}; - [status, groupBy] = checks(opt, nodeName); + [status, groupBy, glmType] = checks(opt, nodeName); if ~status 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 all(ismember(lower(groupBy), {'contrast'})) + switch glmType + case 'one_sample_t_test' - % collect all con images from all subjects - for iSub = 1:numel(opt.subjects) - subLabel = opt.subjects{iSub}; - conImages{iSub} = findSubjectConImage(opt, subLabel, contrasts); - end + 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 - % TODO further refactoring is possible? - for iCon = 1:numel(contrasts) + case {'two_sample_t_test', 'one_way_anova'} - contrastName = contrasts{iCon}; + if strcmp(glmType, 'two_sample_t_test') + label = '2samplesTTest'; + else + label = '1WayANOVA'; + end - contrastsList{end + 1, 1} = contrastName; - groups{end + 1, 1} = 'ALL'; + contrastsList = getContrastsListForFactorialDesign(opt, nodeName); - msg = sprintf(' Group contrast: "%s"', contrastName); - logger('INFO', msg, 'options', opt, 'filename', mfilename()); + % Sorting is important so that we know in which order + % the groups are entered in the design matrix. + % Otherwise it will be harder to properly design + % the contrast vectors later. + groupColumnHdr = opt.model.bm.getGroupColumnHdrFromDesignMatrix(nodeName); + availableGroups = getAvailableGroups(opt, groupColumnHdr); - icell = allocateSubjectsContrasts(opt, opt.subjects, conImages, iCon); + if strcmp(glmType, 'two_sample_t_test') && numel(availableGroups) > 2 + list = bids.internal.create_unordered_list(availableGroups); + msg = sprintf('Too many groups for 2 samples t-test: "%s"', list); + logger('ERROR', msg, 'options', opt, 'filename', filename, 'id', id); + end - matlabbatch = assignToBatch(matlabbatch, opt, nodeName, contrastName, icell); + for iCon = 1:numel(contrastsList) - end + groups{end + 1} = label; %#ok<*AGROW> - elseif all(ismember(lower(groupBy), {'contrast', 'group'})) + contrastName = contrastsList{iCon}; - groupColumnHdr = groupBy{ismember(lower(groupBy), {'group'})}; + rfxDir = getRFXdir(opt, nodeName, contrastName, label); + overwriteDir(rfxDir, opt); - checkColumnParticipantsTsv(BIDS, groupColumnHdr); + assert(~checkSpmMat(rfxDir, opt)); - availableGroups = unique(BIDS.raw.participants.content.(groupColumnHdr)); + if strcmp(glmType, 'two_sample_t_test') + factorialDesign = returnTwoSampleTTestBatch(rfxDir); + else + factorialDesign = returnOneWayAnovaBatch(rfxDir); + end - for iGroup = 1:numel(availableGroups) + for iGroup = 1:numel(availableGroups) - thisGroup = availableGroups{iGroup}; + 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); + subjectsLabel = returnSubjectLabelInGroup(opt, groupColumnHdr, thisGroup); - % collect all con images from all subjects - for iSub = 1:numel(subjectsLabel) - subLabel = subjectsLabel{iSub}; - conImages{iSub} = findSubjectConImage(opt, subLabel, contrasts); - end + % collect all con images from all subjects + for iSub = 1:numel(subjectsLabel) + subLabel = subjectsLabel{iSub}; + conImages{iSub} = findSubjectConImage(opt, subLabel, contrastName); + end + + msg = sprintf(' Group contrast "%s" for group "%s"', contrastName, thisGroup); + logger('INFO', msg, 'options', opt, 'filename', mfilename()); - for iCon = 1:numel(contrasts) + icell = allocateSubjectsConImages(opt, subjectsLabel, conImages, iCon); - contrastName = contrasts{iCon}; + if strcmp(glmType, 'two_sample_t_test') + if iGroup == 1 + factorialDesign.des.t2.scans1 = icell.scans; - contrastsList{end + 1, 1} = contrastName; - groups{end + 1, 1} = thisGroup; + elseif iGroup == 2 + factorialDesign.des.t2.scans2 = icell.scans; + end + elseif strcmp(glmType, 'one_way_anova') + factorialDesign.des.anova.icell(iGroup).scans = icell.scans; + end - msg = sprintf(' Group contrast "%s" for group "%s"', contrastName, thisGroup); - logger('INFO', msg, 'options', opt, 'filename', mfilename()); + end - icell = allocateSubjectsContrasts(opt, subjectsLabel, conImages, iCon); + mask = getInclusiveMask(opt, nodeName); + factorialDesign.masking.em = {mask}; - matlabbatch = assignToBatch(matlabbatch, opt, nodeName, contrastName, icell, thisGroup); + matlabbatch{end + 1}.spm.stats.factorial_design = factorialDesign; + + matlabbatch = setBatchPrintFigure(matlabbatch, opt, ... + fullfile(rfxDir, ... + designMatrixFigureName(opt, ... + 'before estimation'))); end + end +end + +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 = {}; + + contrasts = getContrastsListForFactorialDesign(opt, nodeName); + + % collect all con images from all subjects + for iSub = 1:numel(opt.subjects) + subLabel = opt.subjects{iSub}; + conImages{iSub} = findSubjectConImage(opt, subLabel, contrasts); + end + + % TODO further refactoring is possible? + for iCon = 1:numel(contrasts) + + contrastName = contrasts{iCon}; + + contrastsList{end + 1, 1} = contrastName; + groups{end + 1, 1} = 'ALL'; + + msg = sprintf(' Group contrast: "%s"', contrastName); + logger('INFO', msg, 'options', opt, 'filename', mfilename()); + + icell = allocateSubjectsConImages(opt, opt.subjects, conImages, iCon); + + matlabbatch = assignToBatch(matlabbatch, opt, nodeName, contrastName, icell); + + end +end + +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 = {}; + + contrasts = getContrastsListForFactorialDesign(opt, nodeName); + + participants = bids.util.tsvread(fullfile(opt.dir.raw, 'participants.tsv')); + groupColumnHdr = opt.model.bm.getGroupColumnHdrFromGroupBy(nodeName, participants); + availableGroups = getAvailableGroups(opt, groupColumnHdr); + + for iGroup = 1:numel(availableGroups) + + thisGroup = availableGroups{iGroup}; + subjectsLabel = returnSubjectLabelInGroup(opt, groupColumnHdr, thisGroup); + + % 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 = allocateSubjectsConImages(opt, subjectsLabel, conImages, iCon); + + matlabbatch = assignToBatch(matlabbatch, opt, nodeName, contrastName, icell, thisGroup); + end end + end -function icell = allocateSubjectsContrasts(opt, subjectsLabel, conImages, iCon) +function subjectsLabel = returnSubjectLabelInGroup(opt, groupColumnHdr, group) + % 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 + % + participants = bids.util.tsvread(fullfile(opt.dir.raw, 'participants.tsv')); + subjectsInGroup = strcmp(participants.(groupColumnHdr), group); + subjectsLabel = regexprep(participants.participant_id(subjectsInGroup), '^sub-', ''); + subjectsLabel = intersect(subjectsLabel, opt.subjects); +end + +function icell = allocateSubjectsConImages(opt, subjectsLabel, conImages, iCon) icell(1).scans = {}; @@ -157,13 +259,16 @@ rfxDir = getRFXdir(opt, nodeName, contrastName, thisGroup); overwriteDir(rfxDir, opt); - assert(exist(fullfile(rfxDir, 'SPM.mat'), 'file') == 0); + opt.verbosity = 0; + assert(~checkSpmMat(rfxDir, opt)); icell(1).levels = 1; assert(iscellstr(icell.scans)); - matlabbatch = returnFactorialDesignBatch(matlabbatch, rfxDir, icell, thisGroup); + factorialDesign = returnFactorialDesignBatch(rfxDir, icell, thisGroup); + + matlabbatch{end + 1}.spm.stats.factorial_design = factorialDesign; mask = getInclusiveMask(opt, nodeName); matlabbatch{end}.spm.stats.factorial_design.masking.em = {mask}; @@ -174,72 +279,103 @@ 'before estimation'))); end -function matlabbatch = returnFactorialDesignBatch(matlabbatch, directory, icell, thisGroup) - +function factorialDesign = returnFactorialDesignBatch(directory, icell, thisGroup) if isempty(thisGroup) thisGroup = 'GROUP'; end - factorialDesign.dir = {directory}; + factorialDesign = commonFactorialDesignBatch(directory); factorialDesign.des.fd.icell = icell; + factorialDesign.des.fd.fact = varianceStruct(); + % GROUP and the number of levels in the group. % If 2 groups, then number of levels = 2 factorialDesign.des.fd.fact.name = thisGroup; factorialDesign.des.fd.fact.levels = numel(icell); - factorialDesign.des.fd.fact.dept = 0; +end - % 1: Assumes that the variance is not the same across groups - % 0: There is no difference in the variance between groups - factorialDesign.des.fd.fact.variance = 1; - factorialDesign.des.fd.fact.gmsca = 0; - factorialDesign.des.fd.fact.ancova = 0; - % factorial_design.cov = []; +function factorialDesign = returnTwoSampleTTestBatch(directory) + factorialDesign = commonFactorialDesignBatch(directory); - factorialDesign = setBatchFactorialDesignImplicitMasking(factorialDesign); + factorialDesign.cov = struct('c', {}, 'cname', {}, 'iCFI', {}, 'iCC', {}); + factorialDesign.multi_cov = struct('files', {}, 'iCFI', {}, 'iCC', {}); - factorialDesign = setBatchFatorialDesignGlobalCalcAndNorm(factorialDesign); + factorialDesign.des.t2 = varianceStruct(); - matlabbatch{end + 1}.spm.stats.factorial_design = factorialDesign; + factorialDesign.des.t2.scans1 = {}; + factorialDesign.des.t2.scans2 = {}; +end + +function factorialDesign = returnOneWayAnovaBatch(directory) + factorialDesign = commonFactorialDesignBatch(directory); + + factorialDesign.cov = struct('c', {}, 'cname', {}, 'iCFI', {}, 'iCC', {}); + factorialDesign.multi_cov = struct('files', {}, 'iCFI', {}, 'iCC', {}); + + factorialDesign.des.anova = varianceStruct(); + factorialDesign.des.anova.icell(1).scans = {}; end -function [status, groupBy] = checks(opt, nodeName) +function factorialDesign = commonFactorialDesignBatch(directory) + factorialDesign.dir = {directory}; + factorialDesign = setBatchFactorialDesignImplicitMasking(factorialDesign); + factorialDesign = setBatchFactorialDesignGlobalCalcAndNorm(factorialDesign); +end - thisNode = opt.model.bm.get_nodes('Name', nodeName); +function value = varianceStruct() - commonMsg = sprintf('for the dataset level node: "%s"', nodeName); + value = struct(); - [status, groupBy] = checkGroupBy(thisNode); + % Assumes groups are independent + value(1).dept = 0; + % 1: Assumes that the variance is not the same across groups + % 0: There is no difference in the variance between groups + value.variance = 1; + value.gmsca = 0; + value.ancova = 0; - % only certain type of model supported for now - designMatrix = opt.model.bm.get_design_matrix('Name', nodeName); - if iscell(designMatrix) || (designMatrix ~= 1) +end - msg = sprintf('Models other than group average not implemented yet %s', commonMsg); - notImplemented(mfilename(), msg, opt); +function [status, groupBy, glmType] = checks(opt, nodeName) - status = false; + commonMsg = sprintf('for the dataset level node: "%s"', nodeName); - end + % TODO refactor + participants = bids.util.tsvread(fullfile(opt.dir.raw, 'participants.tsv')); - datasetLvlDummyContrasts = opt.model.bm.get_dummy_contrasts('Name', nodeName); - if isempty(datasetLvlDummyContrasts) || ~isfield(datasetLvlDummyContrasts, 'Test') + bm = opt.model.bm; - msg = sprintf('Only DummyContrasts are implemented %s', commonMsg); - notImplemented(mfilename(), msg, opt); + status = bm.validateGroupBy(nodeName, participants); - status = false; + [glmType, groupBy] = bm.groupLevelGlmType(nodeName, participants); + % only certain type of model supported for now + if ismember(glmType, {'unknown'}) + % TODO update message to with better info for 2 sample T-Test + msg = sprintf(['Models other than group average / comparisons ', ... + 'not implemented yet %s'], commonMsg); + notImplemented(mfilename(), msg, opt); + status = false; + return end - datasetLvlContrasts = opt.model.bm.get_contrasts('Name', nodeName); - if ~isempty(datasetLvlContrasts) + datasetLvlContrasts = bm.get_contrasts('Name', nodeName); + datasetLvlDummyContrasts = bm.get_dummy_contrasts('Name', nodeName); - msg = sprintf('Contrasts are not yet implemented %s', commonMsg); + if ismember(glmType, {'one_sample_t_test'}) && ... + (not(isempty(datasetLvlContrasts)) || isempty(datasetLvlDummyContrasts)) + msg = sprintf('For one-sample t-test only DummyContrasts are implemented %s', commonMsg); notImplemented(mfilename(), msg, opt); + end + if ismember(glmType, {'one_way_anova', 'two_sample_t_test'}) && ... + (isempty(datasetLvlContrasts) || not(isempty(datasetLvlDummyContrasts))) + msg = sprintf(['For one-way ANOVA or 2 samples t-test ', ... + 'only contrasts are implemented %s'], commonMsg); + notImplemented(mfilename(), msg, opt); end end diff --git a/src/batches/stats/setBatchFatorialDesignGlobalCalcAndNorm.m b/src/batches/stats/setBatchFactorialDesignGlobalCalcAndNorm.m similarity index 54% rename from src/batches/stats/setBatchFatorialDesignGlobalCalcAndNorm.m rename to src/batches/stats/setBatchFactorialDesignGlobalCalcAndNorm.m index b15235bee..886ae2578 100644 --- a/src/batches/stats/setBatchFatorialDesignGlobalCalcAndNorm.m +++ b/src/batches/stats/setBatchFactorialDesignGlobalCalcAndNorm.m @@ -1,8 +1,8 @@ -function factorialDesign = setBatchFatorialDesignGlobalCalcAndNorm(factorialDesign) +function factorialDesign = setBatchFactorialDesignGlobalCalcAndNorm(factorialDesign) % % USAGE:: % - % factorialDesign = setBatchFatorialDesignGlobalCalcAndNorm(factorialDesign) + % factorialDesign = setBatchFactorialDesignGlobalCalcAndNorm(factorialDesign) % % diff --git a/src/batches/stats/setBatchGroupLevelContrasts.m b/src/batches/stats/setBatchGroupLevelContrasts.m index 1db245f2a..825170f3f 100644 --- a/src/batches/stats/setBatchGroupLevelContrasts.m +++ b/src/batches/stats/setBatchGroupLevelContrasts.m @@ -26,15 +26,21 @@ printBatchName('group level contrast estimation', opt); - [groupGlmType, designMatrix, groupBy] = groupLevelGlmType(opt, nodeName); + bm = opt.model.bm; + participants = bids.util.tsvread(fullfile(opt.dir.raw, 'participants.tsv')); + + [groupGlmType, groupBy] = bm.groupLevelGlmType(nodeName, participants); switch groupGlmType case 'one_sample_t_test' + groupColumnHdr = bm.getGroupColumnHdrFromGroupBy(nodeName, participants); + availableGroups = getAvailableGroups(opt, groupColumnHdr); + contrastsList = getContrastsListForFactorialDesign(opt, nodeName); - if all(ismember(lower(groupBy), {'contrast'})) + if numel(groupBy) == 1 && ismember(lower(groupBy), {'contrast'}) for j = 1:numel(contrastsList) @@ -46,12 +52,7 @@ end - elseif all(ismember(lower(groupBy), {'contrast', 'group'})) - - participants = bids.util.tsvread(fullfile(opt.dir.raw, 'participants.tsv')); - - groupColumnHdr = groupBy{ismember(lower(groupBy), {'group'})}; - availableGroups = unique(participants.(groupColumnHdr)); + elseif numel(groupBy) == 2 && any(ismember(groupBy, fieldnames(participants))) for j = 1:numel(contrastsList) @@ -71,31 +72,53 @@ end - case 'two_sample_t_test' + case {'two_sample_t_test', 'one_way_anova' } + + % T test for ANOVA + % + % Loop over the subject level contrasts passed + % through the Edge filter. + % Then generate the between group contrasts. + + groupColumnHdr = bm.getGroupColumnHdrFromDesignMatrix(nodeName); + availableGroups = getAvailableGroups(opt, groupColumnHdr); + + edge = bm.get_edge('Destination', nodeName); + contrastsList = edge.Filter.contrast; - designMatrix = removeIntercept(designMatrix); + thisContrast = bm.get_contrasts('Name', nodeName); - if ismember(lower(designMatrix), {'group'}) - % TODO will this ignore the contrasts define at other levels and not - % passed through the filter ? - edge = opt.model.bm.get_edge('Destination', nodeName); - contrastsList = edge.Filter.contrast; + if strcmp(groupGlmType, 'two_sample_t_test') + label = '2samplesTTest'; + else + label = '1WayANOVA'; end for j = 1:numel(contrastsList) - thisContrast = opt.model.bm.get_contrasts('Name', nodeName); + spmMatFile = fullfile(getRFXdir(opt, ... + nodeName, ... + contrastsList{j}, ... + label), ... + 'SPM.mat'); - spmMatFile = fullfile(getRFXdir(opt, nodeName, contrastsList{j}), 'SPM.mat'); + for iCon = 1:numel(thisContrast) - if ~opt.dryRun - assert(exist(spmMatFile, 'file') == 2); - end + convec = generateConStruct(thisContrast{iCon}, ... + availableGroups, ... + groupColumnHdr); - for iCon = 1:numel(thisContrast) - consess{iCon}.tcon.name = thisContrast{iCon}.Name; - consess{iCon}.tcon.convec = thisContrast{iCon}.Weights; - consess{iCon}.tcon.sessrep = 'none'; + if strcmp(thisContrast{iCon}.Test, 't') + consess{iCon}.tcon.name = thisContrast{iCon}.Name; %#ok<*AGROW> + consess{iCon}.tcon.sessrep = 'none'; + consess{iCon}.tcon.convec = convec; + + elseif strcmp(thisContrast{iCon}.Test, 'F') + consess{iCon}.fcon.name = thisContrast{iCon}.Name; + consess{iCon}.fcon.sessrep = 'none'; + consess{iCon}.fcon.convec = convec; + + end end matlabbatch = setBatchContrasts(matlabbatch, opt, spmMatFile, consess); @@ -123,3 +146,40 @@ matlabbatch = setBatchContrasts(matlabbatch, opt, spmMatFile, consess); end + +function convec = generateConStruct(thisContrast, availableGroups, groupColumnHdr) + % Sort conditions and weights + ConditionList = strrep(thisContrast.ConditionList, ... + [groupColumnHdr, '.'], ... + ''); + [ConditionList, I] = sort(ConditionList); + + if strcmp(thisContrast.Test, 't') + + Weights = thisContrast.Weights(I); + + % Create contrast vectors by what was passed in the model + convec = zeros(size(availableGroups)); + for iGroup = 1:numel(availableGroups) + index = strcmp(availableGroups{iGroup}, ConditionList); + if any(index) + convec(iGroup) = Weights(index); + end + end + + elseif strcmp(thisContrast.Test, 'F') + + Weights = thisContrast.Weights(I, :); + + % Create contrast vectors by what was passed in the mode + convec = zeros(size(Weights, 1), numel(availableGroups)); + for iGroup = 1:numel(availableGroups) + index = strcmp(availableGroups{iGroup}, ConditionList); + if any(index) + convec(:, iGroup) = Weights(:, index); + end + end + + end + +end diff --git a/src/batches/stats/setBatchSubjectLevelContrasts.m b/src/batches/stats/setBatchSubjectLevelContrasts.m index 8e62a9dcf..707ea9e19 100644 --- a/src/batches/stats/setBatchSubjectLevelContrasts.m +++ b/src/batches/stats/setBatchSubjectLevelContrasts.m @@ -4,7 +4,7 @@ % % USAGE:: % - % matlabbatch = setBatchSubjectLevelContrasts(matlabbatch, opt, subLabel, funcFWHM) + % matlabbatch = setBatchSubjectLevelContrasts(matlabbatch, opt, subLabel, nodeName) % % :param matlabbatch: % :type matlabbatch: structure @@ -27,20 +27,21 @@ printBatchName('subject level contrasts specification', opt); - spmMatFile = fullfile(getFFXdir(subLabel, opt), 'SPM.mat'); - if noSPMmat(opt, subLabel, spmMatFile) + if ~checkSpmMat(getFFXdir(subLabel, opt), opt) return end + spmMatFile = fullfile(getFFXdir(subLabel, opt), 'SPM.mat'); load(spmMatFile, 'SPM'); - model = opt.model.bm; + bm = opt.model.bm; + bm.validateConstrasts(); % Create Contrasts if nargin < 4 || isempty(nodeName) - contrasts = specifyContrasts(model, SPM); + contrasts = specifyContrasts(bm, SPM); else - contrasts = specifyContrasts(model, SPM, nodeName); + contrasts = specifyContrasts(bm, SPM, nodeName); end validateContrasts(contrasts); diff --git a/src/batches/stats/setBatchSubjectLevelGLMSpec.m b/src/batches/stats/setBatchSubjectLevelGLMSpec.m index 8750bb98f..c6f42a62b 100644 --- a/src/batches/stats/setBatchSubjectLevelGLMSpec.m +++ b/src/batches/stats/setBatchSubjectLevelGLMSpec.m @@ -36,7 +36,7 @@ logger('ERROR', msg, 'filename', mfilename(), 'id', 'missingRawDir'); end - opt.model.bm.getModelType(); + bm = opt.model.bm; printBatchName('specify subject level fmri model', opt); @@ -87,11 +87,11 @@ fmri_spec.fact = struct('name', {}, 'levels', {}); - fmri_spec.mthresh = opt.model.bm.getInclusiveMaskThreshold(); + fmri_spec.mthresh = bm.getInclusiveMaskThreshold(); - fmri_spec.bases.hrf.derivs = opt.model.bm.getHRFderivatives(); + fmri_spec.bases.hrf.derivs = bm.getHRFderivatives(); - fmri_spec.cvi = opt.model.bm.getSerialCorrelationCorrection(); + fmri_spec.cvi = bm.getSerialCorrelationCorrection(); %% List scans, onsets, confounds for each task / session / run subLabel = regexify(subLabel); @@ -167,7 +167,7 @@ % multicondition selection fmri_spec.sess(iSpmSess).cond = struct('name', {}, 'onset', {}, 'duration', {}); - fmri_spec.sess(iSpmSess).hpf = opt.model.bm.getHighPassFilter(); + fmri_spec.sess(iSpmSess).hpf = bm.getHighPassFilter(); end @@ -185,7 +185,7 @@ matlabbatch{end + 1}.spm.stats.fmri_design = fmri_spec; else - node = opt.model.bm.get_root_node; + node = bm.get_root_node(); fmri_spec.mask = {getInclusiveMask(opt, node.Name, BIDS, subLabel)}; matlabbatch{end + 1}.spm.stats.fmri_spec = fmri_spec; diff --git a/src/batches/stats/setBatchTwoSampleTTest.m b/src/batches/stats/setBatchTwoSampleTTest.m deleted file mode 100644 index b3ece85ca..000000000 --- a/src/batches/stats/setBatchTwoSampleTTest.m +++ /dev/null @@ -1,210 +0,0 @@ -function [matlabbatch, contrastsList, groupList] = setBatchTwoSampleTTest(varargin) - % - % Sets up a group level GLM specification for a 2 sample T test - % - % USAGE:: - % - % matlabbatch = setBatchTwoSampleTTest(matlabbatch, opt, nodeName) - % - % - % :param matlabbatch: matlabbatch to append to. - % :type matlabbatch: cell - % - % :type opt: structure - % :param opt: Options chosen for the analysis. - % See :func:`checkOptions`. - % - % :param nodeName: - % :type nodeName: char - % - % :return: :matlabbatch: (cell) The matlabbatch ready to run the spm job - % - - % (C) Copyright 2022 bidspm developers - - % TODO so far this assumes contrasts are only passed through Edge.Filter - % this is too restrictive - % could want to do 2 sample t-tests on all contrast from lower levels - - args = inputParser; - - addRequired(args, 'matlabbatch', @iscell); - addRequired(args, 'opt', @isstruct); - addRequired(args, 'nodeName', @ischar); - - parse(args, varargin{:}); - - matlabbatch = args.Results.matlabbatch; - opt = args.Results.opt; - nodeName = args.Results.nodeName; - - %% - groupList = {}; - - %% - printBatchName('specify group level paired T-test fmri model', opt); - - [BIDS, opt] = getData(opt, opt.dir.preproc); - - node = opt.model.bm.get_nodes('Name', nodeName); - - % Assumes that group belonging was entered like this in the node contrast - % - % "ConditionList": [ - % "Group.blind", - % "Group.control" - % ], - group1 = regexp(node.Contrasts{1}.ConditionList{1}, '\.', 'split'); - group2 = regexp(node.Contrasts{1}.ConditionList{2}, '\.', 'split'); - - % for now we assume we can read the subject group belonging - % from the participant TSV in the raw dataset - % and from the same column - assert(strcmp(group1{1}, group2{1})); - groupField = group1{1}; - - group1 = group1{2}; - group2 = group2{2}; - - availableGroups = unique(BIDS.raw.participants.content.(groupField)); - - if any(~ismember({group1, group2}, availableGroups)) - error(['Some requested group is not present: %s.', ... - '\nAvailable groups in participants.tsv:\n%s'], ... - strjoin({group1, group2}, ', '), ... - bids.internal.create_unordered_list(availableGroups)); - end - - contrastsList = getContrastsListForDatasetLevel(opt, nodeName); - - % collect con images - for iSub = 1:numel(opt.subjects) - subLabel = opt.subjects{iSub}; - conImages{iSub} = findSubjectConImage(opt, subLabel, contrastsList); - end - - % set up the batch - for iCon = 1:numel(contrastsList) - - groupList{end + 1, 1} = 'ALL'; - - contrastName = contrastsList{iCon}; - - rfxDir = getRFXdir(opt, nodeName, contrastName, groupList{end}); - - overwriteDir(rfxDir, opt); - - factorialDesign.dir = {rfxDir}; - factorialDesign.des.t2.scans1 = {}; - factorialDesign.des.t2.scans2 = {}; - - for iSub = 1:numel(opt.subjects) - - subLabel = opt.subjects{iSub}; - - idx = strcmp(BIDS.raw.participants.content.participant_id, ['sub-' subLabel]); - participantGroup = BIDS.raw.participants.content.(groupField){idx}; - - if numel(contrastsList) == 1 - file = conImages{iSub}; - else - file = conImages{iSub}{iCon}; - end - if isempty(file) - continue - end - - if strcmp (participantGroup, group1) - factorialDesign.des.t2.scans1{end + 1, 1} = file; - - elseif strcmp (participantGroup, group2) - factorialDesign.des.t2.scans2{end + 1, 1} = file; - - else - error(['Unknown group: %s.', ... - '\nAvailable groups in participants.tsv:\n%s'], ... - participantGroup, ... - bids.internal.create_unordered_list(availableGroups)); - end - - printProcessingSubject(iSub, subLabel, opt); - msg = sprintf(' %s', file); - logger('INFO', msg, 'options', opt, 'filename', mfilename()); - - end - - mask = getInclusiveMask(opt, nodeName); - factorialDesign.masking.em = {mask}; - - factorialDesign.des.t2.dept = 0; - factorialDesign.des.t2.variance = 1; - factorialDesign.des.t2.gmsca = 0; - factorialDesign.des.t2.ancova = 0; - - factorialDesign.cov = struct('c', {}, 'cname', {}, 'iCFI', {}, 'iCC', {}); - factorialDesign.multi_cov = struct('files', {}, 'iCFI', {}, 'iCC', {}); - - factorialDesign = setBatchFactorialDesignImplicitMasking(factorialDesign); - - factorialDesign = setBatchFatorialDesignGlobalCalcAndNorm(factorialDesign); - - matlabbatch{end + 1}.spm.stats.factorial_design = factorialDesign; - - matlabbatch = setBatchPrintFigure(matlabbatch, opt, ... - fullfile(rfxDir, ... - designMatrixFigureName(opt, ... - 'before estimation'))); - - end - -end - -function contrastsList = getContrastsListForDatasetLevel(opt, nodeName) - - % TODO refactor with "getContrastsListForFactorialDesign" from - % "setBatchFactorialDesign" - - contrastsList = {}; - - % we try to grab the contrasts list from the Edge.Filter - % otherwise we dig in this in Node - % or the previous one to find the list of contrasts - - % If we assume that all the contrast we want to loop over - % are specified in the Filter of the edge of the BIDS stats model - % - % { - % "Source": "subject_level", - % "Destination": "between_groups", - % "Filter": { - % "contrast": [ - % "bar", "foo" - % ] - % } - % } - - edge = opt.model.bm.get_edge('Destination', nodeName); - - if isfield(edge, 'Filter') && ... - isfield(edge.Filter, 'contrast') && ... - ~isempty(edge.Filter.contrast) - - contrastsList = edge.Filter.contrast; - - else - - % TODO?? can't imagine a 2 sample t-test with dummy contrasts - node = opt.model.bm.get_nodes('Name', nodeName); - - % if no specific dummy contrasts mentioned also include all contrasts from previous levels - % or if contrasts are mentioned we grab them - if isfield(node, 'Contrasts') - tmp = getContrastsList(opt.model.bm, nodeName); - for i = 1:numel(tmp) - contrastsList{end + 1} = tmp{i}.Name; - end - end - - end - -end diff --git a/src/bids/checkColumnParticipantsTsv.m b/src/bids/checkColumnParticipantsTsv.m deleted file mode 100644 index a1f7b2950..000000000 --- a/src/bids/checkColumnParticipantsTsv.m +++ /dev/null @@ -1,21 +0,0 @@ -function checkColumnParticipantsTsv(BIDS, columnHdr) - % - % Check that a given column exists in participants.tsv - % - % USAGE:: - % - % checkColumnParticipantsTsv(BIDS, columnHdr) - % - % - - % (C) Copyright 2023 bidspm developers - - columns = fieldnames(BIDS.raw.participants.content); - if ~ismember(columnHdr, fieldnames(BIDS.raw.participants.content)) - msg = sprintf(['Could not find column "%s" in participants.tsv.\n', ... - 'Available columns are: %s'], ... - columnHdr, ... - bids.internal.create_unordered_list(columns)); - logger('ERROR', msg, 'filename', mfilename, 'id', 'missingColumn'); - end -end diff --git a/src/bids/getAvailableGroups.m b/src/bids/getAvailableGroups.m new file mode 100644 index 000000000..b7cb46153 --- /dev/null +++ b/src/bids/getAvailableGroups.m @@ -0,0 +1,25 @@ +function availableGroups = getAvailableGroups(opt, groupColumnHdr) + + % (C) Copyright 2024 bidspm developers + if isempty(groupColumnHdr) + availableGroups = {}; + return + end + + tsvFile = fullfile(opt.dir.raw, 'participants.tsv'); + + assert(exist(tsvFile, 'file') == 2); + + participants = bids.util.tsvread(tsvFile); + + columns = fieldnames(participants); + if ~ismember(groupColumnHdr, columns) + msg = sprintf(['Could not find column "%s" in participants.tsv.\n', ... + 'Available columns are: %s'], ... + columnHdr, ... + bids.internal.create_unordered_list(columns)); + logger('ERROR', msg, 'filename', mfilename, 'id', 'missingColumn'); + end + + availableGroups = sort(unique(participants.(groupColumnHdr))); +end diff --git a/src/bids_model/BidsModel.m b/src/bids_model/BidsModel.m index 3722a2e1c..0c49baeae 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); @@ -346,6 +347,82 @@ end + function [type, groupBy] = groupLevelGlmType(obj, nodeName, participants) + % + % Return type of GLM for a dataset level node. + % + % USAGE:: + % + % [type, groupBy] = bm.groupLevelGlmType(obj, nodeName, participants) + % + % :param nodeName: + % :type nodeName: char + % + % :param participants: content of participants.tsv + % :type participants: struct + % + if nargin < 3 + participants = struct(); + end + + node = obj.get_nodes('Name', nodeName); + + groupBy = node.GroupBy; + type = 'unknown'; + + if ~strcmpi(node.Level, 'Dataset') + return + end + + designMatrix = node.Model.X; + + if isnumeric(designMatrix) && designMatrix == 1 + type = 'one_sample_t_test'; + + elseif iscell(designMatrix) && numel(designMatrix) == 2 + + groupColumHdr = getGroupColumnHdrFromDesignMatrix(obj, nodeName); + + if isempty(groupColumHdr) || ~isfield(participants, groupColumHdr) + type = 'unknown'; + else + levels = participants.(groupColumHdr); + switch numel(unique(levels)) + case 1 + type = 'one_sample_t_test'; + case 2 + type = 'two_sample_t_test'; + otherwise + type = 'one_way_anova'; + end + end + + end + + end + + function groupColumnHdr = getGroupColumnHdrFromDesignMatrix(obj, nodeName) + node = obj.get_nodes('Name', nodeName); + designMatrix = node.Model.X; + if iscell(designMatrix) && numel(designMatrix) == 2 + designMatrix = removeIntercept(designMatrix); + groupColumnHdr = designMatrix{1}; + else + groupColumnHdr = ''; + end + end + + function groupColumnHdr = getGroupColumnHdrFromGroupBy(obj, nodeName, participants) + node = obj.get_nodes('Name', nodeName); + groupBy = node.GroupBy; + groupColumnHdr = intersect(groupBy, fieldnames(participants)); + if isempty(groupColumnHdr) + groupColumnHdr = ''; + else + groupColumnHdr = groupColumnHdr{1}; + end + end + function obj = addConfoundsToDesignMatrix(obj, varargin) % % Add some typical confounds to the design matrix of bids stat model. @@ -496,6 +573,108 @@ end + function validateRootNode(obj) + thisNode = obj.get_root_node(); + + if ismember(lower(thisNode.Level), {'session', 'subject'}) + + notImplemented(mfilename(), ... + '"session" and "subject" level Node not implemented yet'); + + elseif ismember(lower(thisNode.Level), {'dataset'}) + + msg = sprintf(['Your model seems to be having dataset Node at its root\n.', ... + 'Validate it: ', ... + 'https://bids-standard.github.io/stats-models/validator.html\n']); + obj.tolerant = false; + id = 'wrongLevel'; + bidsModelError(obj, id, msg); + + end + end + + function status = validateGroupBy(obj, nodeName, participants) + % + % USAGE:: + % + % bm = bm.validateGroupBy(nodeName, participants); + % + % Only certain type of GroupBy supported for now for each level + % + % This helps doing some defensive programming. + % + + % (C) Copyright 2022 bidspm developers + + status = false; + + node = obj.get_nodes('Name', nodeName); + if isempty(node) + return + end + + groupBy = sort(node.GroupBy); + + if nargin < 3 + participants = struct(); + end + + extraVar = fieldnames(participants); + + switch lower(node.Level) + + case 'run' + + % only certain type of GroupBy supported for now + supportedGroupBy = {'["run", "subject"]', '["run", "session", "subject"]'}; + if ismember('run', groupBy) && ... + all(ismember(groupBy, {'run', 'session', 'subject'})) + status = true; + end + + case 'session' + + % only certain type of GroupBy supported for now + supportedGroupBy = {'["contrast", "session", "subject"]'}; + if numel(groupBy) == 3 && ... + all(ismember(groupBy, {'contrast', 'session', 'subject'})) + status = true; + end + + case 'subject' + + supportedGroupBy = {'["contrast", "subject"]'}; + if numel(groupBy) == 2 && ... + all(ismember(groupBy, {'contrast', 'subject'})) + status = true; + end + + case 'dataset' + + % only certain type of GroupBy supported for now + supportedGroupBy = {'["contrast"]', ... + '["contrast", "x"] for "x" being a participant.tsv column name.'}; + + if numel(groupBy) == 1 && ismember(lower(groupBy), {'contrast'}) + status = true; + elseif numel(groupBy) == 2 && any(ismember(lower(groupBy), {'contrast'})) && ... + iscellstr(extraVar) && numel(extraVar) > 0 && any(ismember(groupBy, extraVar)) + status = true; + end + + end + + if status + return + end + + template = 'only "GroupBy": %s supported %s node level'; + msg = sprintf(template, ... + bids.internal.create_unordered_list(supportedGroupBy), ... + node.Level); + notImplemented(mfilename(), msg); + end + function validateConstrasts(obj) % validate all contrasts spec in the model @@ -643,3 +822,12 @@ function bidsModelError(obj, id, msg) end end + +function designMatrix = removeIntercept(designMatrix) + % + % remove intercept because SPM includes it anyway + % + + isIntercept = cellfun(@(x) (numel(x) == 1) && (x == 1), designMatrix, 'UniformOutput', true); + designMatrix(isIntercept) = []; +end diff --git a/src/bids_model/checkContrast.m b/src/bids_model/checkContrast.m index 76b52c003..d40c514f7 100644 --- a/src/bids_model/checkContrast.m +++ b/src/bids_model/checkContrast.m @@ -15,8 +15,7 @@ model.validate_constrasts(node); - if ~ismember(lower(node.Level), {'run', 'session', 'subject'}) && ... - ~isTtest(node.Contrasts(iCon)) + if ismember(lower(node.Level), {'session'}) && ~isTtest(node.Contrasts(iCon)) notImplemented(mfilename(), 'Only t test implemented for Contrasts'); contrast = []; return diff --git a/src/bids_model/checkGroupBy.m b/src/bids_model/checkGroupBy.m deleted file mode 100644 index 93a0bcac5..000000000 --- a/src/bids_model/checkGroupBy.m +++ /dev/null @@ -1,71 +0,0 @@ -function [status, groupBy] = checkGroupBy(node) - % - % Only certain type of GroupBy supported for now for each level - % - % This helps doing some defensive programming - % - - % (C) Copyright 2022 bidspm developers - - status = true; - groupBy = sort(node.GroupBy); - - switch lower(node.Level) - - case 'run' - - % only certain type of GroupBy supported for now - if ~ismember('run', groupBy) || ... - ~all(ismember(groupBy, {'run', 'session', 'subject'})) - - status = false; - - supportedGroupBy = {'["run", "subject"]', ... - '["run", "session", "subject"]'}; - - end - - case 'session' - - if ~(numel(groupBy) == 3) || ... - ~all(ismember(groupBy, {'contrast', 'session', 'subject'})) - - status = false; - - supportedGroupBy = {'["contrast", "session", "subject"]'}; - - end - - case 'subject' - - if ~(numel(groupBy) == 2) || ... - not(all([ismember('contrast', groupBy) ismember('subject', groupBy)])) - - status = false; - - supportedGroupBy = {'["contrast", "subject"]'}; - - end - - case 'dataset' - - % only certain type of GroupBy supported for now - if numel(groupBy) > 2 || ~all(ismember(lower(groupBy), {'contrast', 'group'})) - - status = false; - - supportedGroupBy = {'["contrast"]', '["contrast", "group"]'}; - - end - - end - - if ~status - template = 'only "GroupBy": %s supported %s node level'; - msg = sprintf(template, ... - bids.internal.create_unordered_list(supportedGroupBy), ... - node.Level); - notImplemented(mfilename(), msg); - end - -end diff --git a/src/bids_model/getContrastsList.m b/src/bids_model/getContrastsList.m index 1655f9cb5..331b3125e 100644 --- a/src/bids_model/getContrastsList.m +++ b/src/bids_model/getContrastsList.m @@ -1,4 +1,4 @@ -function contrastsList = getContrastsList(model, node) +function contrastsList = getContrastsList(model, node, participants) % % Get list of names of Contrast from this Node % or gets its from the previous Nodes @@ -18,6 +18,10 @@ % (C) Copyright 2022 bidspm developers + if nargin < 3 + participants = struct(); + end + contrastsList = {}; if ischar(node) @@ -44,7 +48,7 @@ % TODO relax those assumptions % assumptions - assert(checkGroupBy(node)); + assert(model.validateGroupBy(node.Name, participants)); assert(node.Model.X == 1); contrastsList = getContrastsFromParentNode(model, node); diff --git a/src/bids_model/getContrastsListForFactorialDesign.m b/src/bids_model/getContrastsListForFactorialDesign.m index 48310c792..80d468778 100644 --- a/src/bids_model/getContrastsListForFactorialDesign.m +++ b/src/bids_model/getContrastsListForFactorialDesign.m @@ -14,10 +14,16 @@ % (C) Copyright 2022 bidspm developers - % assuming we want to only average at the group level - if opt.model.bm.get_design_matrix('Name', nodeName) == 1 + % assuming we want to only average / comparisons at the group level + participants = bids.util.tsvread(fullfile(opt.dir.raw, 'participants.tsv')); - edge = opt.model.bm.get_edge('Destination', nodeName); + bm = opt.model.bm; + + groupGlmType = bm.groupLevelGlmType(nodeName, participants); + + if ismember(groupGlmType, {'one_sample_t_test', 'two_sample_t_test', 'one_way_anova'}) + + edge = bm.get_edge('Destination', nodeName); if isfield(edge, 'Filter') && ... isfield(edge.Filter, 'contrast') && ... @@ -28,14 +34,16 @@ else % this assumes DummyContrasts exist - contrastsList = getDummyContrastsList(opt.model.bm, nodeName); + contrastsList = getDummyContrastsList(bm, nodeName, participants); - node = opt.model.bm.get_nodes('Name', nodeName); + node = bm.get_nodes('Name', nodeName); % if no specific dummy contrasts mentioned also include all contrasts from previous levels % or if contrasts are mentioned we grab them - if ~isfield(node.DummyContrasts, 'Contrasts') || isfield(node, 'Contrasts') - tmp = getContrastsList(opt.model.bm, nodeName); + if ~isfield(node, 'DummyContrasts') || ... + ~isfield(node.DummyContrasts, 'Contrasts') || ... + isfield(node, 'Contrasts') + tmp = getContrastsList(bm, nodeName, participants); for i = 1:numel(tmp) contrastsList{end + 1} = tmp{i}.Name; end @@ -46,9 +54,12 @@ else commonMsg = sprintf('for the dataset level node: "%s"', nodeName); - msg = sprintf('Models other than group average not implemented yet %s', commonMsg); + msg = sprintf(['Models other than group average / comparisons ', ... + 'not implemented yet %s'], commonMsg); notImplemented(mfilename(), msg, opt); + contrastsList = {}; + end end diff --git a/src/bids_model/getDummyContrastsList.m b/src/bids_model/getDummyContrastsList.m index 3c03e6bb2..aefe88137 100644 --- a/src/bids_model/getDummyContrastsList.m +++ b/src/bids_model/getDummyContrastsList.m @@ -1,4 +1,4 @@ -function dummyContrastsList = getDummyContrastsList(model, node) +function dummyContrastsList = getDummyContrastsList(model, node, participants) % % Get list of names of DummyContrast from this Node or gets its from the % previous Nodes @@ -18,6 +18,10 @@ % (C) Copyright 2022 bidspm developers + if nargin < 3 + participants = struct(); + end + dummyContrastsList = {}; if ischar(node) @@ -27,19 +31,19 @@ end end - if isfield(node.DummyContrasts, 'Contrasts') + if isfield(node, 'DummyContrasts') && isfield(node.DummyContrasts, 'Contrasts') dummyContrastsList = node.DummyContrasts.Contrasts; else - assert(checkGroupBy(node)); + assert(model.validateGroupBy(node.Name, participants)); switch lower(node.Level) case 'run' - dummyContrastsList = node.Model.X; + dummyContrastsList = node.Model.HRF.Variables; case {'session', 'subject'} diff --git a/src/bids_model/getInclusiveMask.m b/src/bids_model/getInclusiveMask.m index 8f9f15232..cc735a337 100644 --- a/src/bids_model/getInclusiveMask.m +++ b/src/bids_model/getInclusiveMask.m @@ -13,13 +13,15 @@ % (C) Copyright 2022 bidspm developers + bm = opt.model.bm; + if nargin < 2 - [mask, nodeName] = opt.model.bm.getModelMask(); + [mask, nodeName] = bm.getModelMask(); else - [mask, nodeName] = opt.model.bm.getModelMask('Name', nodeName); + [mask, nodeName] = bm.getModelMask('Name', nodeName); end - node = opt.model.bm.get_nodes('Name', nodeName); + node = bm.get_nodes('Name', nodeName); % TODO refactor with bidsResults part for checking background for montage if isstruct(mask) diff --git a/src/messages/noSPMmat.m b/src/messages/noSPMmat.m deleted file mode 100644 index 01351c860..000000000 --- a/src/messages/noSPMmat.m +++ /dev/null @@ -1,45 +0,0 @@ -function status = noSPMmat(varargin) - % - % USAGE:: - % - % status = noSPMmat(opt, subLabel, spmMatFile) - % - % :type opt: structure - % :param opt: Options chosen for the analysis. - % See :func:`checkOptions`. - % :param subLabel: - % :type subLabel: char - % :param spmMatFile: - % :type spmMatFile: path - % - % :return: status - % :rtype: logical - % - - % (C) Copyright 2022 bidspm developers - - args = inputParser; - - addRequired(args, 'opt', @isstruct); - addRequired(args, 'subLabel', @ischar); - addRequired(args, 'spmMatFile', @ischar); - - parse(args, varargin{:}); - - status = false; - - % could be refactored with an equivalent function of group level result - if ~(exist(args.Results.spmMatFile, 'file') == 2) - - status = true; - - msg = sprintf(['No SPM.mat found for subject %s in folder:\n%s\n', ... - 'Run bidsFFX(''specify'', opt)\n'], ... - args.Results.subLabel, ... - bids.internal.format_path(spm_fileparts(args.Results.spmMatFile))); - id = 'noSpecifiedModel'; - logger('WARNING', msg, 'id', id, 'filename', mfilename(), 'options', args.Results.opt); - - end - -end diff --git a/src/stats/group_level/getRFXdir.m b/src/stats/group_level/getRFXdir.m index f3da84d13..a64e3d152 100644 --- a/src/stats/group_level/getRFXdir.m +++ b/src/stats/group_level/getRFXdir.m @@ -66,11 +66,8 @@ end sub = 'ALL'; - if ~isempty(contrastName) - thisNode = opt.model.bm.get_nodes('Name', nodeName); - if all(ismember(lower(thisNode.GroupBy), {'contrast', 'group'})) && ~isempty(thisGroup) - sub = thisGroup; - end + if ~isempty(contrastName) && ~isempty(thisGroup) + sub = thisGroup; end glmDirName = ['sub-', sub, '_', glmDirName]; diff --git a/src/stats/group_level/groupLevelGlmType.m b/src/stats/group_level/groupLevelGlmType.m deleted file mode 100644 index e7f28feef..000000000 --- a/src/stats/group_level/groupLevelGlmType.m +++ /dev/null @@ -1,25 +0,0 @@ -function [type, designMatrix, groupBy] = groupLevelGlmType(opt, nodeName) - % - - % (C) Copyright 2022 bidspm developers - - designMatrix = opt.model.bm.getBidsDesignMatrix('Name', nodeName); - - node = opt.model.bm.get_nodes('Name', nodeName); - - groupBy = node.GroupBy; - - if isnumeric(designMatrix) && designMatrix == 1 - - type = 'one_sample_t_test'; - - elseif iscell(designMatrix) && numel(designMatrix) == 2 - - type = 'two_sample_t_test'; - - else - type = 'unknown'; - - end - -end diff --git a/src/stats/subject_level/convertOnsetTsvToMat.m b/src/stats/subject_level/convertOnsetTsvToMat.m index 78b1f73fd..d0f6bc9b6 100644 --- a/src/stats/subject_level/convertOnsetTsvToMat.m +++ b/src/stats/subject_level/convertOnsetTsvToMat.m @@ -113,8 +113,10 @@ opt.model.bm = BidsModel('file', opt.model.file); end - varToConvolve = opt.model.bm.getVariablesToConvolve(); - designMatrix = opt.model.bm.getBidsDesignMatrix(); + bm = opt.model.bm; + + varToConvolve = bm.getVariablesToConvolve(); + designMatrix = bm.getBidsDesignMatrix(); designMatrix = removeIntercept(designMatrix); % conditions to be filled according to the conditions present in each run @@ -125,7 +127,7 @@ condToModel.idx = 1; % TODO get / apply transformers from a specific node - transformers = opt.model.bm.getBidsTransformers(); + transformers = bm.getBidsTransformers(); tsv.content = bids.transformers(transformers, tsv.content); for iVar = 1:numel(varToConvolve) diff --git a/src/stats/subject_level/createAndReturnCounfoundMatFile.m b/src/stats/subject_level/createAndReturnCounfoundMatFile.m index 40790cb3d..bc2ff9470 100644 --- a/src/stats/subject_level/createAndReturnCounfoundMatFile.m +++ b/src/stats/subject_level/createAndReturnCounfoundMatFile.m @@ -42,9 +42,11 @@ return end - designMatrix = opt.model.bm.getBidsDesignMatrix(); + bm = opt.model.bm; - transformers = opt.model.bm.getBidsTransformers(); + designMatrix = bm.getBidsDesignMatrix(); + + transformers = bm.getBidsTransformers(); content = bids.transformers(transformers, content); [names, R] = createConfounds(content, designMatrix, opt.glm.maxNbVols); %#ok<*ASGLU> diff --git a/src/stats/subject_level/specifyContrasts.m b/src/stats/subject_level/specifyContrasts.m index e729af96e..3e13cf3e9 100644 --- a/src/stats/subject_level/specifyContrasts.m +++ b/src/stats/subject_level/specifyContrasts.m @@ -57,6 +57,10 @@ node = node{1}; end + if ismember(lower(node.Level), {'session', 'subject'}) && ~model.validateGroupBy(node.Name) + continue + end + [contrasts, counter] = specifyDummyContrasts(model, node, contrasts, counter); switch lower(node.Level) @@ -65,15 +69,9 @@ [contrasts, counter] = specifyRunLvlContrasts(model, node, contrasts, counter); case 'session' - if ~checkGroupBy(node) - continue - end [contrasts, counter] = specifySessionLvlContrasts(model, node, contrasts, counter); case 'subject' - if ~checkGroupBy(node) - continue - end [contrasts, counter] = specifySubLvlContrasts(model, node, contrasts, counter); end @@ -91,7 +89,7 @@ end -function nodeList = getNodeList(model) +function nodeList = getNodeList(bm) % % Return all the nodes present in the bids stats model % - if there is only one, @@ -99,33 +97,33 @@ % nodeList = {}; - if isempty(model.Nodes) + if isempty(bm.Nodes) return end - if numel(model.Nodes) == 1 - if iscell(model.Nodes) - nodeList = model.Nodes{1}; + if numel(bm.Nodes) == 1 + if iscell(bm.Nodes) + nodeList = bm.Nodes{1}; else % should not be necessary % but in case nodes were not coerced to cell % during test set up - nodeList = model.Nodes(1); + nodeList = bm.Nodes(1); end return end - if isempty(model.Edges) - model = model.get_edges_from_nodes(); + if isempty(bm.Edges) + bm = bm.get_edges_from_nodes(); end - for i = 1:numel(model.Edges) - nodeList{end + 1} = model.Edges{i}.Source; %#ok<*AGROW> - nodeList{end + 1} = model.Edges{i}.Destination; + for i = 1:numel(bm.Edges) + nodeList{end + 1} = bm.Edges{i}.Source; %#ok<*AGROW> + nodeList{end + 1} = bm.Edges{i}.Destination; end nodeList = unique(nodeList); for iNode = 1:length(nodeList) - nodeList{iNode} = model.get_nodes('Name', nodeList{iNode}); + nodeList{iNode} = bm.get_nodes('Name', nodeList{iNode}); end end diff --git a/src/stats/subject_level/specifyDummyContrasts.m b/src/stats/subject_level/specifyDummyContrasts.m index e1cb742e4..217c8bd93 100644 --- a/src/stats/subject_level/specifyDummyContrasts.m +++ b/src/stats/subject_level/specifyDummyContrasts.m @@ -26,13 +26,10 @@ return end - if ismember(node.Level, {'session'}) && ~checkGroupBy(node) + if ismember(node.Level, {'Subject', 'Session'}) && ~model.validateGroupBy(node.Name) return end - if strcmp(node.Level, 'subject') && ~checkGroupBy(node) - return - end - if ismember(node.Level, {'dataset'}) + if ismember(node.Level, {'Dataset'}) % see setBatchGroupLevelContrasts return end diff --git a/src/stats/utils/checkSpmMat.m b/src/stats/utils/checkSpmMat.m new file mode 100644 index 000000000..d205e82f6 --- /dev/null +++ b/src/stats/utils/checkSpmMat.m @@ -0,0 +1,19 @@ +function status = checkSpmMat(dir, opt, strict) + % (C) Copyright 2019 bidspm developers + status = exist(fullfile(dir, 'SPM.mat'), 'file'); + if nargin < 2 + opt = struct('verbosity', 2); + end + if nargin < 3 + strict = false; + end + if ~status + msg = sprintf('\nCould not find a SPM.mat file in directory:\n%s\n', dir); + id = 'noSpmMatFile'; + if strict + logger('ERROR', msg, 'id', id, 'options', opt, 'filename', mfilename); + else + logger('WARNING', msg, 'id', id, 'options', opt, 'filename', mfilename); + end + end +end diff --git a/src/stats/utils/getContrastNb.m b/src/stats/utils/getContrastNb.m index 55ae510f1..773740ebb 100644 --- a/src/stats/utils/getContrastNb.m +++ b/src/stats/utils/getContrastNb.m @@ -58,6 +58,8 @@ id = 'noMatchingContrastName'; logger('WARNING', msg, 'id', id, 'filename', mfilename()); + contrastNb = []; + return end diff --git a/src/stats/subject_level/removeIntercept.m b/src/stats/utils/removeIntercept.m similarity index 100% rename from src/stats/subject_level/removeIntercept.m rename to src/stats/utils/removeIntercept.m diff --git a/src/workflows/bidsChangeSuffix.m b/src/workflows/bidsChangeSuffix.m index debf1b2e2..2ce47532f 100644 --- a/src/workflows/bidsChangeSuffix.m +++ b/src/workflows/bidsChangeSuffix.m @@ -8,7 +8,7 @@ function bidsChangeSuffix(varargin) % See :func:`checkOptions`. % :type opt: structure % - % :param newSuffix: TODO: add checks on newSuffix to make sure it only contains [a-zA-Z0-9] + % :param newSuffix: % :type newSuffix: char % % :param filter: structure to decide which files to include. Default: @@ -36,6 +36,8 @@ function bidsChangeSuffix(varargin) % (C) Copyright 2022 bidspm developers + % TODO: add checks on newSuffix to make sure it only contains [a-zA-Z0-9] + args = inputParser; addRequired(args, 'opt', @isstruct); diff --git a/src/workflows/bidsRename.m b/src/workflows/bidsRename.m index aaeac18ea..8f187926a 100644 --- a/src/workflows/bidsRename.m +++ b/src/workflows/bidsRename.m @@ -80,7 +80,6 @@ function renameFileAndUpdateMetadata(opt, data, newFilename, json, createdFiles) return end - % TODO write test for this if ~opt.rename.overwrite && exist(outputFile, 'file') || ... ismember(newFilename, createdFiles) diff --git a/src/workflows/lesion/bidsLesionOverlapMap.m b/src/workflows/lesion/bidsLesionOverlapMap.m index 27c57e52b..a3d852933 100755 --- a/src/workflows/lesion/bidsLesionOverlapMap.m +++ b/src/workflows/lesion/bidsLesionOverlapMap.m @@ -21,8 +21,6 @@ function bidsLesionOverlapMap(opt) % (C) Copyright 2021 bidspm developers - % TODO add test - if checkToolbox('ALI', 'verbose', opt.verbosity > 0) opt = setFields(opt, ALI_my_defaults()); end diff --git a/src/workflows/preproc/bidsWholeBrainFuncMask.m b/src/workflows/preproc/bidsWholeBrainFuncMask.m index 81f7cfa9a..11aa5d6c5 100644 --- a/src/workflows/preproc/bidsWholeBrainFuncMask.m +++ b/src/workflows/preproc/bidsWholeBrainFuncMask.m @@ -9,8 +9,6 @@ % (C) Copyright 2020 bidspm developers - % TODO add test - opt.pipeline.type = 'preproc'; opt.dir.input = opt.dir.preproc; diff --git a/src/workflows/roi/bidsRoiBasedGLM.m b/src/workflows/roi/bidsRoiBasedGLM.m index ddf39d820..8f558746a 100644 --- a/src/workflows/roi/bidsRoiBasedGLM.m +++ b/src/workflows/roi/bidsRoiBasedGLM.m @@ -76,16 +76,14 @@ logger('INFO', msg, 'options', opt, 'filename', mfilename()); outputDir = getFFXdir(subLabel, opt); - - spmFile = fullfile(outputDir, 'SPM.mat'); - - if noSPMmat(opt, subLabel, spmFile) + if ~checkSpmMat(outputDir, opt) continue end - load(spmFile, 'SPM'); + spmMatFile = fullfile(outputDir, 'SPM.mat'); + load(spmMatFile, 'SPM'); model = mardo(SPM); - eventSpec = getEventSpecificationRoiGlm(spmFile, opt.model.file); + eventSpec = getEventSpecificationRoiGlm(spmMatFile, opt.model.file); subTempDir = fullfile(tempDir, ['sub-' subLabel]); spm_mkdir(subTempDir); diff --git a/src/workflows/stats/bidsConcatBetaTmaps.m b/src/workflows/stats/bidsConcatBetaTmaps.m index fc0708b79..d84069da7 100644 --- a/src/workflows/stats/bidsConcatBetaTmaps.m +++ b/src/workflows/stats/bidsConcatBetaTmaps.m @@ -119,7 +119,7 @@ function bidsConcatBetaTmaps(opt, deleteTmaps) spm_save(fullfile(ffxDir, bf.filename), tsvContent); - % TODO in the dev branch make those output filenames "BIDS derivatives" compliant + % TODO make those output filenames "BIDS derivatives" compliant % beta maps bf = bids.File(nameStructure); bf.extension = '.nii'; diff --git a/src/workflows/stats/bidsFFX.m b/src/workflows/stats/bidsFFX.m index d3b2c5a6b..4747d7e72 100644 --- a/src/workflows/stats/bidsFFX.m +++ b/src/workflows/stats/bidsFFX.m @@ -103,7 +103,7 @@ case 'estimate' - if noSPMmat(opt, subLabel, fullfile(outputDir, 'SPM.mat')) + if ~checkSpmMat(outputDir, opt) continue end matlabbatch = setAction(action, matlabbatch, BIDS, opt, subLabel); @@ -115,11 +115,9 @@ case 'contrasts' - if isempty(nodeName) - matlabbatch = setBatchSubjectLevelContrasts(matlabbatch, opt, subLabel); - else - matlabbatch = setBatchSubjectLevelContrasts(matlabbatch, opt, subLabel, nodeName); - end + opt.model.bm.validateConstrasts(); + + matlabbatch = setBatchSubjectLevelContrasts(matlabbatch, opt, subLabel, nodeName); end @@ -161,23 +159,11 @@ function checkRootNode(opt) % This only concerns 'specify' and 'specifyAndEstimate' % - thisNode = opt.model.bm.get_root_node; - - if ismember(lower(thisNode.Level), {'session', 'subject'}) - - notImplemented(mfilename(), ... - '"session" and "subject" level Node not implemented yet'); - - elseif ismember(lower(thisNode.Level), {'dataset'}) - - msg = sprintf(['Your model seems to be having dataset Node at its root\n.', ... - 'Validate it: https://bids-standard.github.io/stats-models/validator.html\n']); - id = 'wrongLevel'; - logger('ERROR', msg, 'id', id, 'filename', mfilename()); - - end + bm = opt.model.bm; + bm.validateRootNode(); - checkGroupBy(thisNode); + thisNode = bm.get_root_node; + bm.validateGroupBy(thisNode.Name); % should not be necessary % mostly in case users did not validate the model inputs diff --git a/src/workflows/stats/bidsModelSelection.m b/src/workflows/stats/bidsModelSelection.m index 24bbc5c8e..e2b6f1bc2 100644 --- a/src/workflows/stats/bidsModelSelection.m +++ b/src/workflows/stats/bidsModelSelection.m @@ -190,15 +190,10 @@ end ffxDir = getFFXdir(subLabel, opt); - - spmMatFile = spm_select('FPList', ffxDir, 'SPM.mat'); - - if isempty(spmMatFile) - msg = sprintf('no SPM.mat found in:\n%s\n\n', ffxDir); - id = 'noSPMmat'; - logger('WARNING', msg, 'id', id, 'filename', mfilename()); + if not(checkSpmMat(ffxDir, opt)) continue end + spmMatFile = spm_select('FPList', ffxDir, 'SPM.mat'); msg = struct('Subject', subLabel); msg.Model = names{iModel}; diff --git a/src/workflows/stats/bidsRFX.m b/src/workflows/stats/bidsRFX.m index 659c5cde2..2bdbf786d 100644 --- a/src/workflows/stats/bidsRFX.m +++ b/src/workflows/stats/bidsRFX.m @@ -59,19 +59,20 @@ % TODO refactor % - extract function for anat and mask computation - % - merge rfx and ffx into a single "stats" workflow if ismember(lower(action), {'meananatandmask', 'rfx', 'contrasts'}) opt.dir.output = fullfile(opt.dir.stats, 'derivatives', 'bidspm-groupStats'); opt.dir.jobs = fullfile(opt.dir.output, 'jobs', strjoin(opt.taskName, '')); end + bm = opt.model.bm; + if ismember(lower(action), {'rfx', 'contrasts'}) % TODO add possibility to pass several nodeNames at once if ~isempty(nodeName) - datasetNodes = opt.model.bm.get_nodes('Name', nodeName); + datasetNodes = bm.get_nodes('Name', nodeName); else - datasetNodes = opt.model.bm.get_nodes('Level', 'Dataset'); + datasetNodes = bm.get_nodes('Level', 'Dataset'); end if isstruct(datasetNodes) datasetNodes = {datasetNodes}; @@ -89,25 +90,25 @@ case 'rfx' - for i = 1:numel(datasetNodes) + [~, opt] = getData(opt, opt.dir.preproc); - msg = sprintf('\n PROCESSING NODE: %s\n', nodeName); - logger('INFO', msg, 'options', opt, 'filename', mfilename()); + participants = bids.util.tsvread(fullfile(opt.dir.raw, 'participants.tsv')); + + for i = 1:numel(datasetNodes) matlabbatch = {}; nodeName = datasetNodes{i}.Name; - switch groupLevelGlmType(opt, nodeName) + msg = sprintf('\n PROCESSING NODE: %s\n', nodeName); + logger('INFO', msg, 'options', opt, 'filename', mfilename()); - case 'one_sample_t_test' + switch bm.groupLevelGlmType(nodeName, participants) + + case {'one_sample_t_test', 'two_sample_t_test', 'one_way_anova'} [matlabbatch, contrastsList, groups] = setBatchFactorialDesign(matlabbatch, ... opt, ... datasetNodes{i}.Name); - case 'two_sample_t_test' - [matlabbatch, contrastsList, groups] = setBatchTwoSampleTTest(matlabbatch, ... - opt, ... - datasetNodes{i}.Name); otherwise msg = sprintf('Node %s has has model type I cannot handle.\n', nodeName); notImplemented(mfilename(), msg); @@ -132,6 +133,8 @@ case 'contrasts' + bm.validateConstrasts(); + for i = 1:numel(datasetNodes) matlabbatch = setBatchGroupLevelContrasts(matlabbatch, opt, datasetNodes{i}.Name); saveAndRunWorkflow(matlabbatch, 'contrasts_rfx', opt); diff --git a/src/workflows/stats/bidsResults.m b/src/workflows/stats/bidsResults.m index 8de9e3562..d675a1159 100644 --- a/src/workflows/stats/bidsResults.m +++ b/src/workflows/stats/bidsResults.m @@ -159,7 +159,9 @@ opt.dir.output = opt.dir.stats; - modelResults = opt.model.bm.getResults(); + bm = opt.model.bm; + + modelResults = bm.getResults(); if ~isempty(modelResults) opt.results = modelResults; end @@ -183,7 +185,7 @@ % loop through the steps to compute for each contrast mentioned for each node for iRes = 1:length(opt.results) - node = opt.model.bm.get_nodes('Name', opt.results(iRes).nodeName); + node = bm.get_nodes('Name', opt.results(iRes).nodeName); if isempty(node) @@ -393,9 +395,7 @@ tmp.dir = getFFXdir(subLabel, opt); - status = checkSpmMat(tmp.dir, opt); - - if ~status + if ~checkSpmMat(tmp.dir, opt) return end @@ -443,10 +443,14 @@ matlabbatch = {}; - node = opt.model.bm.get_nodes('Name', opt.results(iRes).nodeName); + bm = opt.model.bm; + + node = bm.get_nodes('Name', opt.results(iRes).nodeName); opt = checkMontage(opt, iRes, node); + participants = bids.util.tsvread(fullfile(opt.dir.raw, 'participants.tsv')); + for i = 1:length(opt.results(iRes).name) result = opt.results(iRes); @@ -459,11 +463,11 @@ logger('WARNING', msg, 'id', id, 'options', opt, 'filename', mfilename()); end - switch groupLevelGlmType(opt, result.nodeName) + [glmType, groupBy] = bm.groupLevelGlmType(result.nodeName, participants); - case 'one_sample_t_test' + switch glmType - [~, ~, groupBy] = groupLevelGlmType(opt, result.nodeName); + case 'one_sample_t_test' if all(ismember(lower(groupBy), {'contrast'})) @@ -473,12 +477,12 @@ matlabbatch = appendToBatch(matlabbatch, opt, result); + % TODO make more general than just with group elseif all(ismember(lower(groupBy), {'contrast', 'group'})) - participants = bids.util.tsvread(fullfile(opt.dir.raw, 'participants.tsv')); - + % TODO make more general than just with group groupColumnHdr = groupBy{ismember(lower(groupBy), {'group'})}; - availableGroups = unique(participants.(groupColumnHdr)); + availableGroups = getAvailableGroups(opt, groupColumnHdr); for iGroup = 1:numel(availableGroups) @@ -493,16 +497,31 @@ end - case 'two_sample_t_test' + case {'two_sample_t_test', 'one_way_anova'} + + edge = bm.get_edge('Destination', result.nodeName); + contrastsList = edge.Filter.contrast; - thisContrast = opt.model.bm.get_contrasts('Name', result.nodeName); + for iCon = 1:numel(contrastsList) + + if strcmp(glmType, 'two_sample_t_test') + label = '2samplesTTest'; + else + label = '1WayANOVA'; + end - result.dir = getRFXdir(opt, result.nodeName, name); + result.dir = getRFXdir(opt, result.nodeName, contrastsList{iCon}, label); + + load(fullfile(result.dir, 'SPM.mat'), 'SPM'); + + result.name = name; + contrastNb = getContrastNb(result, opt, SPM); + result.contrastNb = contrastNb; + + result.name = [bids.internal.camel_case(contrastsList{iCon}) ' - ' name]; - for iCon = 1:numel(thisContrast) - result.name = [thisContrast{iCon}.Name ' - ' name]; - result.contrastNb = iCon; matlabbatch = appendToBatch(matlabbatch, opt, result); + end otherwise @@ -515,18 +534,6 @@ end -function status = checkSpmMat(dir, opt) - status = exist(fullfile(dir, 'SPM.mat'), 'file'); - if ~status - if nargin < 2 - opt = struct('verbosity', 2); - end - msg = sprintf('\nCould not find a SPM.mat file in directory:\n%s\n', dir); - id = 'noSpmMatFile'; - logger('WARNING', msg, 'id', id, 'options', opt, 'filename', mfilename); - end -end - function matlabbatch = appendToBatch(matlabbatch, opt, result) if ~checkSpmMat(result.dir, opt) diff --git a/tests/.gitignore b/tests/.gitignore index f41fa0222..3780f205f 100644 --- a/tests/.gitignore +++ b/tests/.gitignore @@ -10,6 +10,9 @@ tmp */*.json +data/3_groups/derivatives +data/3_groups/bidspm* + data/derivatives/bidspm-*/jobs data/derivatives/bidspm-preproc/sub* diff --git a/tests/Makefile b/tests/Makefile index f52fa6034..9c898e747 100644 --- a/tests/Makefile +++ b/tests/Makefile @@ -11,6 +11,7 @@ clean_local: create_dummy_dataset: clean_local python create_dummy_dataset.py + python create_3_groups_dataset.py bids_examples: rm -rf bids-examples/ diff --git a/tests/create_3_groups_dataset.py b/tests/create_3_groups_dataset.py new file mode 100644 index 000000000..c42274851 --- /dev/null +++ b/tests/create_3_groups_dataset.py @@ -0,0 +1,89 @@ +from __future__ import annotations + +import shutil +from pathlib import Path + +import pandas as pd + +from utils import ( + create_preproc_func_vismotion, + create_raw_func_vismotion, + touch, + write_ds_desc, +) + +ROOT_DIR = Path(__file__).parent.parent +START_DIR = Path(__file__).parent +SUBJECTS = ["ctrl01", "blind01", "relative01", "ctrl02", "blind02", "relative02"] +# SUBJECTS = ["ctrl01", "blind01", "relative01"] +SESSIONS = ["01", "02"] +SESSIONS = ["01"] + + +def main(start_dir, subject_list, session_list): + # Create BIDS directory structure + raw_dir = start_dir / "data" / "3_groups" / "bidspm-raw" + write_ds_desc(raw_dir) + derivative_dir = start_dir / "data" / "3_groups" / "derivatives" + preproc_dir = derivative_dir / "bidspm-preproc" + write_ds_desc(preproc_dir) + stats_dir = derivative_dir / "bidspm-stats" + write_ds_desc(stats_dir) + + create_raw_dataset(raw_dir, subject_list, session_list) + create_preproc_dataset(preproc_dir, subject_list, session_list) + create_stats_dataset(stats_dir, subject_list) + + +def create_raw_dataset(target_dir, subject_list, session_list): + # Create the raw BIDS dataset + create_participants_tsv(target_dir, subject_list) + for sub in subject_list: + for ses in session_list: + create_raw_func_vismotion(target_dir, sub, ses) + + +def create_preproc_dataset(target_dir, subject_list, session_list): + for sub in subject_list: + for ses in session_list: + create_preproc_func_vismotion(target_dir, sub, ses) + + +def create_participants_tsv(bids_dir, subject_list): + bids_dir.parent.mkdir(parents=True, exist_ok=True) + content = { + "participant_id": subject_list, + "diagnostic": [x[:-2] for x in subject_list], + } + df = pd.DataFrame(content) + df.to_csv(bids_dir / "participants.tsv", sep="\t", index=False) + + +def create_stats_dataset(stats_dir, subject_list): + task_list = ["vismotion"] + + for sub in subject_list: + for task in task_list: + this_dir = stats_dir / f"sub-{sub}" / f"task-{task}_space-IXI549Space_FWHM-6" + + touch(this_dir / "mask.nii") + + for i in range(1, 10): + touch(this_dir / f"beta_000{i}.nii") + + for i in range(1, 5): + touch(this_dir / f"spmT_000{i}.nii") + touch(this_dir / f"con_000{i}.nii") + + shutil.copy( + ROOT_DIR / "tests" / "data" / "mat_files" / "SPM.mat", + this_dir / "SPM.mat", + ) + + +if __name__ == "__main__": + main( + start_dir=START_DIR, + subject_list=SUBJECTS, + session_list=SESSIONS, + ) diff --git a/tests/create_dummy_dataset.py b/tests/create_dummy_dataset.py index b5fce77dd..d912d1b3c 100644 --- a/tests/create_dummy_dataset.py +++ b/tests/create_dummy_dataset.py @@ -4,9 +4,15 @@ import shutil from pathlib import Path -import pandas as pd +from utils import ( + ROOT_DIR, + create_events_tsv, + create_preproc_func_vismotion, + create_raw_func_vismotion, + touch, + write_ds_desc, +) -ROOT_DIR = Path(__file__).parent.parent START_DIR = Path(__file__).parent SUBJECTS = ["ctrl01", "blind01", "01"] SESSIONS = ["01", "02"] @@ -14,15 +20,10 @@ HEMISPHERES = ["L", "R"] -def touch(path: Path): - path.parent.mkdir(parents=True, exist_ok=True) - with open(path, "a"): - pass - - def main(start_dir, subject_list, session_list, roi_list, hemispheres): # Create BIDS directory structure raw_dir = start_dir / "data" / "bidspm-raw" + write_ds_desc(start_dir / "data" / "bidspm-raw") derivative_dir = start_dir / "data" / "derivatives" preproc_dir = derivative_dir / "bidspm-preproc" stats_dir = derivative_dir / "bidspm-stats" @@ -47,46 +48,6 @@ def create_raw_dataset(target_dir, subject_list, session_list): create_raw_anat(target_dir, sub) -def create_raw_func_vismotion(target_dir, sub, ses): - suffix = "bold" - task = "vismotion" - this_dir = target_dir / f"sub-{sub}" / f"ses-{ses}" / "func" - this_dir.mkdir(parents=True, exist_ok=True) - - basename = f"sub-{sub}_ses-{ses}_task-{task}" - - create_events_tsv( - filename=this_dir / f"{basename}_run-1_events.tsv", - onset=[2, 4], - duration=[2, 2], - trial_type=["VisMot", "VisStat"], - ) - - for run in range(1, 3): - filename = this_dir / f"{basename}_run-{run}_part-mag_{suffix}.nii" - touch(filename) - filename = this_dir / f"{basename}_run-{run}_part-phase_{suffix}.nii" - touch(filename) - - create_events_tsv( - filename=this_dir / f"{basename}_run-2_events.tsv", - onset=[3, 6], - duration=[2, 2], - trial_type=["VisStat", "VisMot"], - ) - - touch(this_dir / f"{basename}_acq-1p60mm_run-1_{suffix}.nii") - - for run in [1, 2]: - create_events_tsv( - filename=this_dir / f"{basename}_acq-1p60mm_run-{run}_events.tsv", - onset=[4, 8], - duration=[2, 2], - trial_type=["VisMot", "VisStat"], - ) - touch(this_dir / f"{basename}_acq-1p60mm_dir-PA_run-{run}_{suffix}.nii") - - def create_raw_func_vislocalizer(target_dir, sub, ses): suffix = "bold" task = "vislocalizer" @@ -104,16 +65,6 @@ def create_raw_func_vislocalizer(target_dir, sub, ses): ) -def create_events_tsv(filename, onset, duration, trial_type): - events = { - "onset": onset, - "duration": duration, - "trial_type": trial_type, - } - df = pd.DataFrame(events) - df.to_csv(filename, sep="\t", index=False) - - def create_raw_func_rest(target_dir, sub, ses): suffix = "bold" task = "rest" @@ -192,40 +143,6 @@ def create_preproc_dataset(target_dir, subject_list, session_list): create_preproc_anat(target_dir, sub) -def create_preproc_func_vismotion(target_dir, sub, ses): - suffix = "bold" - task = "vismotion" - this_dir = target_dir / f"sub-{sub}" / f"ses-{ses}" / "func" - this_dir.mkdir(parents=True, exist_ok=True) - - for acq_entity in ["", "_acq-1p60mm"]: - basename = f"sub-{sub}_ses-{ses}_task-{task}{acq_entity}_part-mag" - - if ses == "01": - touch(this_dir / f"{basename}_space-individual_desc-mean_{suffix}.nii") - touch(this_dir / f"{basename}_space-IXI549Space_desc-mean_{suffix}.nii") - touch(this_dir / f"mean_{basename}_{suffix}.nii") - - for run in range(1, 3): - basename = f"sub-{sub}_ses-{ses}_task-{task}{acq_entity}_run-{run}_part-mag" - - filename = this_dir / f"rp_{basename}_{suffix}.txt" - shutil.copy(ROOT_DIR / "tests" / "data" / "tsv_files" / "rp.txt", filename) - - shutil.copy( - ROOT_DIR / "tests" / "data" / "tsv_files" / "rp.tsv", - this_dir / f"{basename}_desc-confounds_regressors.tsv", - ) - - desc_label_list = ["preproc", "mean", "smth6"] - for desc in desc_label_list: - touch(this_dir / f"{basename}_space-individual_desc-{desc}_{suffix}.nii") - touch(this_dir / f"{basename}_space-IXI549Space_desc-{desc}_{suffix}.nii") - - touch(this_dir / f"{basename}_space-individual_desc-stc_{suffix}.nii") - touch(this_dir / f"{basename}_space-IXI549Space_desc-brain_mask.nii") - - def create_preproc_func_vislocalizer(target_dir, sub, ses): suffix = "bold" task = "vislocalizer" @@ -277,7 +194,6 @@ def create_preproc_anat(target_dir, sub): suffix = "T1w" this_dir = target_dir / f"sub-{sub}" / f"ses-{ses}" / "anat" - this_dir.mkdir(parents=True, exist_ok=True) basename = f"sub-{sub}_ses-{ses}" @@ -312,15 +228,14 @@ def create_stats_dataset(stats_dir, subject_list): for sub in subject_list: for task in task_list: this_dir = stats_dir / f"sub-{sub}" / f"task-{task}_space-IXI549Space_FWHM-6" - this_dir.mkdir(parents=True, exist_ok=True) + + touch(this_dir / "mask.nii") shutil.copy( ROOT_DIR / "tests" / "data" / "mat_files" / "SPM.mat", this_dir / "SPM.mat", ) - touch(this_dir / "mask.nii") - for i in range(1, 10): touch(this_dir / f"beta_000{i}.nii") @@ -336,15 +251,14 @@ def create_stats_dataset(stats_dir, subject_list): / f"sub-{sub}" / f"task-{task}_space-IXI549Space_FWHM-6_node-globalSignal" ) - this_dir.mkdir(parents=True, exist_ok=True) + + touch(this_dir / "mask.nii") shutil.copy( ROOT_DIR / "tests" / "data" / "mat_files" / "SPM.mat", this_dir / "SPM.mat", ) - touch(this_dir / "mask.nii") - for i in range(1, 11): touch(this_dir / f"beta_000{i}.nii") @@ -356,7 +270,6 @@ def create_stats_dataset(stats_dir, subject_list): def create_roi_directories(roi_dir, roi_list, subject_list, hemispheres): # Create ROI directories roi_group_dir = roi_dir / "group" - roi_group_dir.mkdir(parents=True, exist_ok=True) for roi in roi_list: touch(roi_group_dir / f"{roi}_mask.nii") diff --git a/tests/data/bidspm-raw/dataset_description.json b/tests/data/bidspm-raw/dataset_description.json index d1673b1fe..3b58ac5b3 100755 --- a/tests/data/bidspm-raw/dataset_description.json +++ b/tests/data/bidspm-raw/dataset_description.json @@ -1,5 +1,5 @@ { - "BIDSVersion": "1.6.0", + "BIDSVersion": "1.9.0", "DatasetType": "raw", "Name": "bidspm raw" } 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..e947c36e2 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": { @@ -114,13 +114,14 @@ }, { "Name": "ctrl_gt_blind", + "Description": "Check that non sorted conditions are handled properly.", "ConditionList": [ - "Group.blind", - "Group.ctrl" + "Group.ctrl", + "Group.blind" ], "Weights": [ - -1, - 1 + 1, + -1 ], "Test": "t" } diff --git a/tests/data/models/model-vislocalizerSeveralDatasetLevels_smdl.json b/tests/data/models/model-vislocalizerSeveralDatasetLevels_smdl.json index 546c80439..bc419f8cb 100644 --- a/tests/data/models/model-vislocalizerSeveralDatasetLevels_smdl.json +++ b/tests/data/models/model-vislocalizerSeveralDatasetLevels_smdl.json @@ -171,5 +171,19 @@ ] } } + ], + "Edges": [ + { + "Source": "run_level", + "Destination": "subject_level" + }, + { + "Source": "subject_level", + "Destination": "simple contrast" + }, + { + "Source": "subject_level", + "Destination": "complex contrast" + } ] } diff --git a/tests/data/models/model-vismotionNoOverWrite_smdl.json b/tests/data/models/model-vismotionNoOverWrite_smdl.json index 1aa9a9ff6..0c66cdf62 100644 --- a/tests/data/models/model-vismotionNoOverWrite_smdl.json +++ b/tests/data/models/model-vismotionNoOverWrite_smdl.json @@ -191,7 +191,7 @@ { "Level": "Dataset", "Name": "between_groups", - "Description": "2 sample t-test of the all_olf condition", + "Description": "2 sample t-test", "GroupBy": [ "contrast" ], @@ -199,7 +199,7 @@ "Type": "glm", "X": [ 1, - "group" + "Group" ], "Options": { "Mask": { @@ -227,5 +227,35 @@ } ] } + ], + "Edges": [ + { + "Source": "run_level", + "Destination": "subject_level" + }, + { + "Source": "subject_level", + "Destination": "simple contrast" + }, + { + "Source": "subject_level", + "Destination": "complex contrast" + }, + { + "Source": "subject_level", + "Destination": "within_group" + }, + { + "Source": "subject_level", + "Destination": "between_groups", + "Filter": { + "contrast": [ + "VisMot", + "VisStat", + "VisMot_gt_VisStat", + "VisStat_gt_VisMot" + ] + } + } ] } diff --git a/tests/data/models/model-vismotionOneWayANOVA_smdl.json b/tests/data/models/model-vismotionOneWayANOVA_smdl.json new file mode 100644 index 000000000..584e2b43b --- /dev/null +++ b/tests/data/models/model-vismotionOneWayANOVA_smdl.json @@ -0,0 +1,202 @@ +{ + "Name": "vismotion", + "BIDSModelVersion": "1.0.0", + "Description": "contrasts for the motion dataset", + "Input": { + "task": [ + "vismotion" + ], + "space": [ + "IXI549Space" + ] + }, + "Nodes": [ + { + "Level": "Run", + "Name": "run_level", + "GroupBy": [ + "run", + "subject" + ], + "Model": { + "Type": "glm", + "X": [ + "trial_type.VisMot", + "trial_type.VisStat", + "trans_?", + "rot_?" + ], + "HRF": { + "Variables": [ + "trial_type.VisMot", + "trial_type.VisStat" + ], + "Model": "spm" + } + }, + "Contrasts": [ + { + "Name": "VisMot_gt_VisStat", + "ConditionList": [ + "trial_type.VisMot", + "trial_type.VisStat" + ], + "Weights": [ + 1, + -1 + ], + "Test": "t" + }, + { + "Name": "VisStat_gt_VisMot", + "ConditionList": [ + "trial_type.VisMot", + "trial_type.VisStat" + ], + "Weights": [ + -1, + 1 + ], + "Test": "t" + } + ] + }, + { + "Level": "Subject", + "Name": "subject_level", + "GroupBy": [ + "contrast", + "subject" + ], + "Model": { + "Type": "glm", + "X": [ + 1 + ] + }, + "DummyContrasts": { + "Test": "t" + } + }, + { + "Level": "Dataset", + "Name": "dataset_level", + "GroupBy": [ + "contrast", + "diagnostic" + ], + "Model": { + "Type": "glm", + "X": [ + 1 + ] + }, + "DummyContrasts": { + "Test": "t" + } + }, + { + "Level": "Dataset", + "Name": "between_groups", + "Description": "one way anova", + "GroupBy": [ + "contrast" + ], + "Model": { + "Type": "glm", + "X": [ + 1, + "diagnostic" + ], + "Software": { + "bidspm": { + "Results": [ + { + "name": [ + "relative_gt_ctrl", + "blind_gt_relative" + ], + "p": 0.05, + "MC": "none", + "png": true, + "binary": false, + "nidm": false, + "montage": { + "do": false + } + } + ] + } + } + }, + "Contrasts": [ + { + "Name": "relative_gt_ctrl", + "Description": "3 groups but not sorted alphabetically whereas the groups are entered alphabetically in the design matrix.", + "ConditionList": [ + "diagnostic.relative", + "diagnostic.ctrl", + "diagnostic.blind" + ], + "Weights": [ + 1, + -1, + 0 + ], + "Test": "t" + }, + { + "Name": "blind_gt_relative", + "Description": "Only 2 groups should be OK but needs to be handled properly: should give a vector [-1, 0, 1]", + "ConditionList": [ + "diagnostic.blind", + "diagnostic.relative" + ], + "Weights": [ + -1, + 1 + ], + "Test": "t" + }, + { + "Name": "F_test", + "ConditionList": [ + "diagnostic.blind", + "diagnostic.relative" + ], + "Weights": [ + [ + 1, + 0 + ], + [ + 0, + 1 + ] + ], + "Test": "F" + } + ] + } + ], + "Edges": [ + { + "Source": "run_level", + "Destination": "subject_level" + }, + { + "Source": "subject_level", + "Destination": "dataset_level" + }, + { + "Source": "subject_level", + "Destination": "between_groups", + "Filter": { + "contrast": [ + "VisMot_gt_VisStat", + "VisStat_gt_VisMot" + ] + } + } + ] +} diff --git a/tests/data/models/model-vismotionSeveralDatasetLevel_smdl.json b/tests/data/models/model-vismotionSeveralDatasetLevel_smdl.json index e28a198bb..ea24e9d5a 100644 --- a/tests/data/models/model-vismotionSeveralDatasetLevel_smdl.json +++ b/tests/data/models/model-vismotionSeveralDatasetLevel_smdl.json @@ -160,5 +160,19 @@ ] } } + ], + "Edges": [ + { + "Source": "run_level", + "Destination": "subject_level" + }, + { + "Source": "subject_level", + "Destination": "simple contrast" + }, + { + "Source": "subject_level", + "Destination": "complex contrast" + } ] } diff --git a/tests/tests_batches/preproc/test_setBatchGenerateT1map.m b/tests/tests_batches/preproc/test_setBatchGenerateT1map.m index 7a2911252..161d9ebb4 100644 --- a/tests/tests_batches/preproc/test_setBatchGenerateT1map.m +++ b/tests/tests_batches/preproc/test_setBatchGenerateT1map.m @@ -48,7 +48,7 @@ function test_setBatchGenerateT1map_warning() BIDS = getLayout(opt); matlabbatch = {}; - opt.verbosity = 2; + opt.verbosity = 1; assertWarning(@() setBatchGenerateT1map(matlabbatch, BIDS, opt, subLabel), ... 'setBatchGenerateT1map:missingNonBIDSMetadata'); diff --git a/tests/tests_batches/stats/test_setBatchGroupLevelContrasts.m b/tests/tests_batches/stats/test_setBatchGroupLevelContrasts.m index d91dd517b..7439081ce 100644 --- a/tests/tests_batches/stats/test_setBatchGroupLevelContrasts.m +++ b/tests/tests_batches/stats/test_setBatchGroupLevelContrasts.m @@ -8,7 +8,7 @@ initTestSuite; end -function test_setBatchGroupLevelContrasts_between_groups() +function test_setBatchGroupLevelContrasts_two_sample_ttest() opt = setOptions('vislocalizer', '', 'pipelineType', 'stats'); @@ -22,7 +22,9 @@ function test_setBatchGroupLevelContrasts_between_groups() assertEqual(numel(matlabbatch{1}.spm.stats.con.consess), 2); assertEqual(matlabbatch{1}.spm.stats.con.consess{1}.tcon.name, 'blind_gt_ctrl'); + assertEqual(matlabbatch{1}.spm.stats.con.consess{1}.tcon.convec, [1; -1]); assertEqual(matlabbatch{1}.spm.stats.con.consess{2}.tcon.name, 'ctrl_gt_blind'); + assertEqual(matlabbatch{1}.spm.stats.con.consess{2}.tcon.convec, [-1; 1]); end diff --git a/tests/tests_bids/test_fileFilterForBold.m b/tests/tests_bids/test_fileFilterForBold.m index 50efff7f5..3de226fbb 100644 --- a/tests/tests_bids/test_fileFilterForBold.m +++ b/tests/tests_bids/test_fileFilterForBold.m @@ -13,7 +13,7 @@ function test_fileFilterForBold_events() opt.bidsFilterFile.bold = struct('modality', 'func', 'suffix', 'bold'); - opt.verbosity = 2; + opt.verbosity = 0; opt.taskName = 'foo'; opt.fwhm.func = 6; opt.space = {'IXI549Space'}; @@ -38,7 +38,7 @@ function test_fileFilterForBold_events_aroma() opt.bidsFilterFile.bold = struct('modality', 'func', ... 'suffix', 'bold', ... 'desc', {'smoothAROMAnonaggr'}); - opt.verbosity = 2; + opt.verbosity = 0; opt.taskName = 'foo'; opt.fwhm.func = 6; opt.space = {'MNI152NLin6Asym'}; @@ -61,7 +61,7 @@ function test_fileFilterForBold_events_aroma() function test_fileFilterForBold_confounds() opt.bidsFilterFile.bold = struct('modality', 'func', 'suffix', 'bold'); - opt.verbosity = 2; + opt.verbosity = 0; opt.taskName = 'foo'; opt.fwhm.func = 6; opt.space = {'IXI549Space'}; @@ -84,7 +84,7 @@ function test_fileFilterForBold_confounds() function test_fileFilterForBold_basic() opt.bidsFilterFile.bold = struct('modality', 'func', 'suffix', 'bold'); - opt.verbosity = 2; + opt.verbosity = 0; opt.taskName = 'foo'; opt.fwhm.func = 6; opt.space = {'IXI549Space'}; @@ -111,7 +111,7 @@ function test_fileFilterForBold_desc() opt.bidsFilterFile.bold = struct('modality', 'func', ... 'suffix', 'bold', ... 'desc', {'smoothAROMAnonaggr'}); - opt.verbosity = 2; + opt.verbosity = 0; opt.taskName = 'foo'; opt.fwhm.func = 6; opt.space = {'MNI152NLin6Asym'}; @@ -138,7 +138,7 @@ function test_fileFilterForBold_desc() function test_fileFilterForBold_no_smooth() opt.bidsFilterFile.bold = struct('modality', 'func', 'suffix', 'bold'); - opt.verbosity = 2; + opt.verbosity = 0; opt.taskName = 'foo'; opt.fwhm.func = 0; opt.space = {'IXI549Space'}; @@ -163,7 +163,7 @@ function test_fileFilterForBold_no_smooth() function test_fileFilterForBold_stc() opt.bidsFilterFile.bold = struct('modality', 'func', 'suffix', 'bold'); - opt.verbosity = 2; + opt.verbosity = 0; opt.taskName = 'foo'; opt.fwhm.func = 0; opt.space = {'IXI549Space'}; diff --git a/tests/tests_bids/unit_tests/test_checkColumnParticipantsTsv.m b/tests/tests_bids/unit_tests/test_checkColumnParticipantsTsv.m deleted file mode 100644 index fe2b2b7af..000000000 --- a/tests/tests_bids/unit_tests/test_checkColumnParticipantsTsv.m +++ /dev/null @@ -1,23 +0,0 @@ -function test_suite = test_checkColumnParticipantsTsv %#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_checkColumnParticipantsTsv_basic() - - BIDS.raw.participants.content = struct('foo', 1); - checkColumnParticipantsTsv(BIDS, 'foo'); - -end - -function test_checkColumnParticipantsTsv_fail() - columnHdr = 'foo'; - BIDS.raw.participants.content = struct('bar', 1); - assertExceptionThrown(@() checkColumnParticipantsTsv(BIDS, columnHdr), ... - 'checkColumnParticipantsTsv:missingColumn'); - -end diff --git a/tests/tests_bids_model/test_checkGroupBy.m b/tests/tests_bids_model/test_checkGroupBy.m deleted file mode 100644 index 0385975f8..000000000 --- a/tests/tests_bids_model/test_checkGroupBy.m +++ /dev/null @@ -1,57 +0,0 @@ -% (C) Copyright 2022 bidspm developers - -function test_suite = test_checkGroupBy %#ok<*STOUT> - 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_checkGroupBy_run() - - opt = setOptions('vismotion', {'01' 'ctrl01'}, 'pipelineType', 'stats'); - - opt.model.bm.Nodes{1}.GroupBy = {'subject'}; - - assertWarning(@()checkGroupBy(opt.model.bm.Nodes{1}), ... - 'checkGroupBy:notImplemented'); - - opt.model.bm.Nodes{1}.GroupBy = {'run', 'dataset'}; - - assertWarning(@()checkGroupBy(opt.model.bm.Nodes{1}), ... - 'checkGroupBy:notImplemented'); - -end - -function test_checkGroupBy_subject() - - opt = setOptions('vismotion', {'01' 'ctrl01'}, 'pipelineType', 'stats'); - - opt.model.bm.Nodes{2}.GroupBy = {'subject'}; - - assertWarning(@()checkGroupBy(opt.model.bm.Nodes{2}), ... - 'checkGroupBy:notImplemented'); - - opt.model.bm.Nodes{2}.GroupBy = {'session', 'subject'}; - - assertWarning(@()checkGroupBy(opt.model.bm.Nodes{2}), ... - 'checkGroupBy:notImplemented'); - -end - -function test_checkGroupBy_dataset() - - opt = setOptions('vismotion', {'01' 'ctrl01'}, 'pipelineType', 'stats'); - - opt.model.bm.Nodes{3}.GroupBy = {'subject'}; - - assertWarning(@()checkGroupBy(opt.model.bm.Nodes{3}), ... - 'checkGroupBy:notImplemented'); - - opt.model.bm.Nodes{3}.GroupBy = {'session', 'subject'}; - - assertWarning(@()checkGroupBy(opt.model.bm.Nodes{3}), ... - 'checkGroupBy:notImplemented'); - -end diff --git a/tests/tests_bids_model/test_getOptionsFromModel.m b/tests/tests_bids_model/test_getOptionsFromModel.m index fa71291e0..665a7ce1d 100644 --- a/tests/tests_bids_model/test_getOptionsFromModel.m +++ b/tests/tests_bids_model/test_getOptionsFromModel.m @@ -27,7 +27,7 @@ function test_getOptionsFromModel_no_model() opt.pipeline.type = 'stats'; opt.model.file = ''; opt.tolerant = true; - opt.verbosity = 3; + opt.verbosity = 1; opt.pipeline.isBms = false; assertExceptionThrown(@() getOptionsFromModel(opt), 'getOptionsFromModel:modelFileMissing'); @@ -79,7 +79,7 @@ function test_getOptionsFromModel_task() %% opt.taskName = {'foo'}; - opt.verbosity = 2; + opt.verbosity = 1; opt.tolerant = true; assertWarning(@() getOptionsFromModel(opt), 'getOptionsFromModel:modelOverridesOptions'); @@ -104,7 +104,7 @@ function test_getOptionsFromModel_subject() %% opt.subjects = '02'; - opt.verbosity = 2; + opt.verbosity = 1; assertWarning(@() getOptionsFromModel(opt), 'getOptionsFromModel:modelOverridesOptions'); @@ -147,7 +147,7 @@ function test_getOptionsFromModel_query() %% opt.query.acq = 'foo'; - opt.verbosity = 2; + opt.verbosity = 1; assertWarning(@() getOptionsFromModel(opt), 'getOptionsFromModel:modelOverridesOptions'); diff --git a/tests/tests_bids_model/test_groupLevelGlmType.m b/tests/tests_bids_model/test_groupLevelGlmType.m new file mode 100644 index 000000000..e26dfb2ba --- /dev/null +++ b/tests/tests_bids_model/test_groupLevelGlmType.m @@ -0,0 +1,47 @@ +function test_suite = test_groupLevelGlmType %#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_groupLevelGlmType_basic() + + opt = setOptions('vismotion', {'01' 'ctrl01'}, 'pipelineType', 'stats'); + + model = opt.model.bm; + model.Nodes{3}.GroupBy = {'contrast', 'group'}; + type = model.groupLevelGlmType('dataset_level'); + assertEqual(type, 'one_sample_t_test'); + + model.Nodes{3}.Model.X = {1, 'group'}; + type = model.groupLevelGlmType('dataset_level'); + assertEqual(type, 'unknown'); + + model.Nodes{3}.Model.X = {1, 'diagnostic'}; + participants = struct('participant_id', {{'01', '02'}}, ... + 'diagnostic', {{'ctrl', 'ctrl'}}); + type = model.groupLevelGlmType('dataset_level', participants); + assertEqual(type, 'one_sample_t_test'); + + model.Nodes{3}.Model.X = {1, 'diagnostic'}; + participants = struct('participant_id', {{'01', '02'}}, ... + 'diagnostic', {{'ctrl', 'patient'}}); + type = model.groupLevelGlmType('dataset_level', participants); + assertEqual(type, 'two_sample_t_test'); + + model.Nodes{3}.Model.X = {1, 'group'}; + participants = struct('participant_id', {{'01', '02'}}, ... + 'diagnostic', {{'ctrl', 'patient'}}); + type = model.groupLevelGlmType('dataset_level', participants); + assertEqual(type, 'unknown'); + + model.Nodes{3}.Model.X = {1, 'group'}; + participants = struct('participant_id', {{'01', '02'}}, ... + 'group', {{'ctrl', 'patient', 'relative'}}); + type = model.groupLevelGlmType('dataset_level', participants); + assertEqual(type, 'one_way_anova'); + +end diff --git a/tests/tests_bids_model/test_validateGroupBy.m b/tests/tests_bids_model/test_validateGroupBy.m new file mode 100644 index 000000000..3ca85b5e1 --- /dev/null +++ b/tests/tests_bids_model/test_validateGroupBy.m @@ -0,0 +1,76 @@ +% (C) Copyright 2022 bidspm developers + +function test_suite = test_validateGroupBy %#ok<*STOUT> + 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_validateGroupBy_run() + + opt = setOptions('vismotion', {'01' 'ctrl01'}, 'pipelineType', 'stats'); + + bm = BidsModel('file', opt.model.file); + bm.verbose = false; + nodeName = bm.Nodes{1}.Name; + + bm.Nodes{1}.GroupBy = {'subject'}; + assertWarning(@()bm.validateGroupBy(nodeName), 'BidsModel:notImplemented'); + + bm.Nodes{1}.GroupBy = {'run', 'dataset'}; + assertWarning(@()bm.validateGroupBy(nodeName), 'BidsModel:notImplemented'); + +end + +function test_validateGroupBy_subject() + + opt = setOptions('vismotion', {'01' 'ctrl01'}, 'pipelineType', 'stats'); + bm = BidsModel('file', opt.model.file); + bm.verbose = false; + nodeName = bm.Nodes{2}.Name; + + bm.Nodes{2}.GroupBy = {'subject'}; + assertWarning(@()bm.validateGroupBy(nodeName), 'BidsModel:notImplemented'); + + bm.Nodes{2}.GroupBy = {'session', 'subject'}; + assertWarning(@()bm.validateGroupBy(nodeName), 'BidsModel:notImplemented'); +end + +function test_validateGroupBy_dataset() + + opt = setOptions('vismotion', {'01' 'ctrl01'}, 'pipelineType', 'stats'); + + bm = BidsModel('file', opt.model.file); + bm.verbose = false; + nodeName = bm.Nodes{3}.Name; + + % should be fine + bm.Nodes{3}.GroupBy = {'contrast'}; + status = bm.validateGroupBy(nodeName); + assertEqual(status, true); + + bm.Nodes{3}.GroupBy = {'subject'}; + assertWarning(@()bm.validateGroupBy(nodeName), 'BidsModel:notImplemented'); + + bm.Nodes{3}.GroupBy = {'session', 'subject'}; + assertWarning(@()bm.validateGroupBy(nodeName), 'BidsModel:notImplemented'); + + bm.Nodes{3}.GroupBy = {'session', 'subject', 'foo'}; + assertWarning(@()bm.validateGroupBy(nodeName), 'BidsModel:notImplemented'); + +end + +function test_validateGroupBy_dataset_group_from_participant() + + opt = setOptions('vismotion', {'01' 'ctrl01'}, 'pipelineType', 'stats'); + + bm = BidsModel('file', opt.model.file); + bm.verbose = false; + nodeName = bm.Nodes{3}.Name; + + bm.Nodes{3}.GroupBy = {'contrast', 'diagnostic'}; + bm.validateGroupBy(nodeName, struct('diagnostic', {{'foo', 'bar'}})); + +end diff --git a/tests/tests_slow/tests_batches/stats/test_setBatchFactorialDesign.m b/tests/tests_slow/tests_batches/stats/test_setBatchFactorialDesign.m index 28fbce855..e73fe54b2 100644 --- a/tests/tests_slow/tests_batches/stats/test_setBatchFactorialDesign.m +++ b/tests/tests_slow/tests_batches/stats/test_setBatchFactorialDesign.m @@ -14,6 +14,8 @@ function test_setBatchFactorialDesign_within_group() opt = setOptions('vislocalizer', '', 'pipelineType', 'stats'); + [~, opt] = getData(opt, opt.dir.preproc); + opt.model.file = spm_file(opt.model.file, ... 'basename', ... 'model-vislocalizerWithinGroup_smdl'); @@ -37,6 +39,8 @@ function test_setBatchFactorialDesign_complex() opt = setOptions('vismotion', {'^01'}, 'pipelineType', 'stats'); + [~, opt] = getData(opt, opt.dir.preproc); + opt.model.file = spm_file(opt.model.file, ... 'basename', ... 'model-vismotionSeveralDatasetLevel_smdl'); @@ -83,6 +87,8 @@ function test_setBatchFactorialDesign_basic() opt = setOptions('vismotion', {'01' 'ctrl01'}, 'pipelineType', 'stats'); + [~, opt] = getData(opt, opt.dir.preproc); + datasetNode = opt.model.bm.get_nodes('Level', 'dataset'); opt.model.bm.Nodes{3}.Model.Options = rmfield(opt.model.bm.Nodes{3}.Model.Options, ... 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 new file mode 100644 index 000000000..5d814c2c2 --- /dev/null +++ b/tests/tests_slow/tests_batches/stats/test_setBatchFactorialDesign_3_groups.m @@ -0,0 +1,23 @@ +% (C) Copyright 2020 bidspm developers + +function test_suite = test_setBatchFactorialDesign_3_groups %#ok<*STOUT> + 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_setBatchFactorialDesign_between_3_groups() + + opt = setOptions('3_groups', '', 'pipelineType', 'stats'); + opt.verbosity = 0; + + [~, opt] = getData(opt, opt.dir.preproc); + + matlabbatch = {}; + [matlabbatch, ~, ~] = setBatchFactorialDesign(matlabbatch, ... + opt, ... + 'between_groups'); + +end diff --git a/tests/tests_slow/tests_workflows/stats/test_bidsRFX.m b/tests/tests_slow/tests_workflows/stats/test_bidsRFX.m index 0e78db53c..5ab6bbb38 100644 --- a/tests/tests_slow/tests_workflows/stats/test_bidsRFX.m +++ b/tests/tests_slow/tests_workflows/stats/test_bidsRFX.m @@ -45,104 +45,6 @@ function test_bidsRFX_rfx() end -function test_bidsRFX_no_overwrite_smoke_test() - - markTestAs('slow'); - - opt = setOptions('vislocalizer', '', 'pipelineType', 'stats'); - opt.model.file = spm_file(opt.model.file, ... - 'basename', ... - 'model-vismotionNoOverWrite_smdl'); - opt.model.bm = BidsModel('file', opt.model.file); - - matlabbatch = bidsRFX('RFX', opt); - - expectedNbBatch = 7; - if bids.internal.is_octave() - expectedNbBatch = 5; - end - assertEqual(numel(matlabbatch), expectedNbBatch); - -end - -function test_bidsRFX_within_group_ttest() - - markTestAs('slow'); - - opt = setOptions('vislocalizer', '', 'pipelineType', 'stats'); - opt.model.file = spm_file(opt.model.file, ... - 'basename', ... - 'model-vislocalizerWithinGroup_smdl'); - opt.model.bm = BidsModel('file', opt.model.file); - - matlabbatch = bidsRFX('RFX', opt); - - if bids.internal.is_octave() - assertEqual(numel(matlabbatch), 10); - else - assertEqual(numel(matlabbatch), 14); - end - - % creates 1 batch for - % - specify - % - figure - % - estimate - % - MACS: model space - % - MACS: information criteria - % - review - % - figure - % for each group - - batchOrder = {'stats', 'factorial_design'; ... - 'util', 'print'; ... - 'stats', 'factorial_design'; ... - 'util', 'print'}; - batchOrder = extendBatchOrder(batchOrder); - batchOrder = extendBatchOrder(batchOrder); - summary = batchSummary(matlabbatch); - assertEqual(summary, batchOrder); - - assertEqual(matlabbatch{1}.spm.stats.factorial_design.dir{1}, ... - fileparts(matlabbatch{5}.spm.stats.fmri_est.spmmat{1})); - if bids.internal.is_octave() - assertEqual(matlabbatch{3}.spm.stats.factorial_design.dir{1}, ... - fileparts(matlabbatch{8}.spm.stats.fmri_est.spmmat{1})); - else - assertEqual(matlabbatch{3}.spm.stats.factorial_design.dir{1}, ... - fileparts(matlabbatch{10}.spm.stats.fmri_est.spmmat{1})); - end - -end - -function test_bidsRFX_two_sample_ttest() - - markTestAs('slow'); - - opt = setOptions('vislocalizer', '', 'pipelineType', 'stats'); - opt.model.file = spm_file(opt.model.file, ... - 'basename', ... - 'model-vislocalizer2sampleTTest_smdl'); - opt.model.bm = BidsModel('file', opt.model.file); - - matlabbatch = bidsRFX('RFX', opt); - - summary = batchSummary(matlabbatch); - - % creates 1 batch for (specify, figure, estimate, review, figure) - batchOrder = {'stats', 'factorial_design'; ... - 'util', 'print'}; - batchOrder = extendBatchOrder(batchOrder); - assertEqual(summary, batchOrder); - - assertEqual(matlabbatch{1}.spm.stats.factorial_design.dir{1}, ... - fileparts(matlabbatch{3}.spm.stats.fmri_est.spmmat{1})); - - % 2 blind and 1 ctrl - assertEqual(numel(matlabbatch{1}.spm.stats.factorial_design.des.t2.scans1), 2); - assertEqual(numel(matlabbatch{1}.spm.stats.factorial_design.des.t2.scans2), 1); - -end - function test_bidsRFX_select_datasets_level_to_run() markTestAs('slow'); @@ -168,36 +70,6 @@ function test_bidsRFX_select_datasets_level_to_run() end -function test_bidsRFX_several_datasets_level() - - markTestAs('slow'); - - opt = setOptions('vislocalizer', '', 'pipelineType', 'stats'); - opt.model.file = spm_file(opt.model.file, ... - 'basename', ... - 'model-vislocalizerSeveralDatasetLevels_smdl'); - opt.model.bm = BidsModel('file', opt.model.file); - - matlabbatch = bidsRFX('RFX', opt); - - summary = batchSummary(matlabbatch); - - % only the batches from the last node is returned - % creates 1 batch for (specify, figure, estimate, review, figure) - batchOrder = {'stats', 'factorial_design'; ... - 'util', 'print'}; - batchOrder = extendBatchOrder(batchOrder); - assertEqual(summary, batchOrder); - - nbGroupLevelModelsReturned = 1; - nbBatchPerModel = 7; - if bids.internal.is_octave() - nbBatchPerModel = 5; - end - assertEqual(numel(matlabbatch), nbGroupLevelModelsReturned * nbBatchPerModel); - -end - function test_bidsRFX_rfx_on_empty_dir() markTestAs('slow'); diff --git a/tests/tests_slow/tests_workflows/stats/test_bidsRFX_2_groups.m b/tests/tests_slow/tests_workflows/stats/test_bidsRFX_2_groups.m new file mode 100644 index 000000000..f6415c385 --- /dev/null +++ b/tests/tests_slow/tests_workflows/stats/test_bidsRFX_2_groups.m @@ -0,0 +1,121 @@ +% (C) Copyright 2021 bidspm developers + +function test_suite = test_bidsRFX_2_groups %#ok<*STOUT> + 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_bidsRFX_within_group_ttest() + + markTestAs('slow'); + + opt = setOptions('vislocalizer', '', 'pipelineType', 'stats'); + opt.model.file = spm_file(opt.model.file, ... + 'basename', ... + 'model-vislocalizerWithinGroup_smdl'); + opt.model.bm = BidsModel('file', opt.model.file); + + matlabbatch = bidsRFX('RFX', opt); + + if bids.internal.is_octave() + assertEqual(numel(matlabbatch), 10); + else + assertEqual(numel(matlabbatch), 14); + end + + % creates 1 batch for + % - specify + % - figure + % - estimate + % - MACS: model space + % - MACS: information criteria + % - review + % - figure + % for each group + + batchOrder = {'stats', 'factorial_design'; ... + 'util', 'print'; ... + 'stats', 'factorial_design'; ... + 'util', 'print'}; + batchOrder = extendBatchOrder(batchOrder); + batchOrder = extendBatchOrder(batchOrder); + summary = batchSummary(matlabbatch); + assertEqual(summary, batchOrder); + + [~, folder] = fileparts(matlabbatch{1}.spm.stats.factorial_design.dir{1}); + assert(bids.internal.starts_with(folder, 'sub-blind')); + + [~, folder] = fileparts(matlabbatch{3}.spm.stats.factorial_design.dir{1}); + assert(bids.internal.starts_with(folder, 'sub-ctrl')); + + assertEqual(matlabbatch{1}.spm.stats.factorial_design.dir{1}, ... + fileparts(matlabbatch{5}.spm.stats.fmri_est.spmmat{1})); + if bids.internal.is_octave() + assertEqual(matlabbatch{3}.spm.stats.factorial_design.dir{1}, ... + fileparts(matlabbatch{8}.spm.stats.fmri_est.spmmat{1})); + else + assertEqual(matlabbatch{3}.spm.stats.factorial_design.dir{1}, ... + fileparts(matlabbatch{10}.spm.stats.fmri_est.spmmat{1})); + end + +end + +function test_bidsRFX_two_sample_ttest() + + markTestAs('slow'); + + opt = setOptions('vislocalizer', '', 'pipelineType', 'stats'); + opt.model.file = spm_file(opt.model.file, ... + 'basename', ... + 'model-vislocalizer2sampleTTest_smdl'); + opt.model.bm = BidsModel('file', opt.model.file); + + matlabbatch = bidsRFX('RFX', opt); + + summary = batchSummary(matlabbatch); + + % creates 1 batch for (specify, figure, estimate, review, figure) + batchOrder = {'stats', 'factorial_design'; ... + 'util', 'print'}; + batchOrder = extendBatchOrder(batchOrder); + assertEqual(summary, batchOrder); + + [~, folder] = fileparts(matlabbatch{1}.spm.stats.factorial_design.dir{1}); + assert(bids.internal.starts_with(folder, 'sub-2samplesTTest')); + + assertEqual(matlabbatch{1}.spm.stats.factorial_design.dir{1}, ... + fileparts(matlabbatch{3}.spm.stats.fmri_est.spmmat{1})); + + % 2 blind and 1 ctrl + assertEqual(numel(matlabbatch{1}.spm.stats.factorial_design.des.t2.scans1), 2); + assertEqual(numel(matlabbatch{1}.spm.stats.factorial_design.des.t2.scans2), 1); + +end + +function batchOrder = extendBatchOrder(batchOrder) + if nargin < 1 + batchOrder = {}; + end + extension = {'stats', 'fmri_est'; ... + 'tools', 'MACS'; ... + 'tools', 'MACS'; ... + 'stats', 'review'; ... + 'util', 'print'}; + if bids.internal.is_octave() + extension(2:3, :) = []; + end + batchOrder = cat(1, batchOrder, extension); +end + +function value = batchSummary(matlabbatch) + value = {}; + for i = 1:numel(matlabbatch) + field = fieldnames(matlabbatch{i}.spm); + value{i, 1} = field{1}; %#ok<*AGROW> + field = fieldnames(matlabbatch{i}.spm.(value{i, 1})); + value{i, 2} = field{1}; + end +end diff --git a/tests/tests_slow/tests_workflows/stats/test_bidsRFX_3_groups.m b/tests/tests_slow/tests_workflows/stats/test_bidsRFX_3_groups.m new file mode 100644 index 000000000..aa8dc7e94 --- /dev/null +++ b/tests/tests_slow/tests_workflows/stats/test_bidsRFX_3_groups.m @@ -0,0 +1,93 @@ +% (C) Copyright 2024 bidspm developers + +function test_suite = test_bidsRFX_3_groups %#ok<*STOUT> + 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_bidsRFX_one_way_anova() + + markTestAs('slow'); + + opt = setOptions('3_groups', '', 'pipelineType', 'stats'); + + opt.model.bm = BidsModel('file', opt.model.file); + opt.verbosity = 0; + + matlabbatch = bidsRFX('RFX', opt, 'nodeName', 'between_groups'); + + assertEqual(batchSummary(matlabbatch), expectedBatchOrder()); + + nodeName = 'between_groups'; + contrastName = 'VisMot_gt_VisStat'; + + [~, folder] = fileparts(matlabbatch{1}.spm.stats.factorial_design.dir{1}); + assert(bids.internal.starts_with(folder, 'sub-1WayANOVA')); + + rfxDir = getRFXdir(opt, nodeName, contrastName, '1WayANOVA'); + assertEqual(matlabbatch{1}.spm.stats.factorial_design.dir{1}, rfxDir); + + assertEqual(numel(matlabbatch{1}.spm.stats.factorial_design.des.anova.icell), 3); + assertEqual(numel(matlabbatch{1}.spm.stats.factorial_design.des.anova.icell(1).scans), 2); + assertEqual(fileparts(matlabbatch{5}.spm.stats.fmri_est.spmmat{1}), rfxDir); + if ~bids.internal.is_octave + assertEqual(matlabbatch{6}.spm.tools.MACS.MA_model_space.dir{1}, rfxDir); + end + + nodeName = 'between_groups'; + contrastName = 'VisStat_gt_VisMot'; + rfxDir = getRFXdir(opt, nodeName, contrastName, '1WayANOVA'); + assertEqual(matlabbatch{3}.spm.stats.factorial_design.dir{1}, rfxDir); + assertEqual(numel(matlabbatch{3}.spm.stats.factorial_design.des.anova.icell), 3); + assertEqual(numel(matlabbatch{3}.spm.stats.factorial_design.des.anova.icell(1).scans), 2); + if ~bids.internal.is_octave + assertEqual(fileparts(matlabbatch{10}.spm.stats.fmri_est.spmmat{1}), rfxDir); + assertEqual(matlabbatch{11}.spm.tools.MACS.MA_model_space.dir{1}, rfxDir); + end + +end + +function value = batchSummary(matlabbatch) + value = {}; + for i = 1:numel(matlabbatch) + field = fieldnames(matlabbatch{i}.spm); + value{i, 1} = field{1}; %#ok<*AGROW> + field = fieldnames(matlabbatch{i}.spm.(value{i, 1})); + value{i, 2} = field{1}; + end +end + +function value = expectedBatchOrder() + + value = {'stats', 'factorial_design'; ... + 'util', 'print'; ... + 'stats', 'factorial_design'; ... + 'util', 'print'; ... + 'stats', 'fmri_est'; ... + 'tools', 'MACS'; ... + 'tools', 'MACS'; ... + 'stats', 'review'; ... + 'util', 'print'; ... + 'stats', 'fmri_est'; ... + 'tools', 'MACS'; ... + 'tools', 'MACS'; ... + 'stats', 'review'; ... + 'util', 'print' }; + + if bids.internal.is_octave + value = {'stats', 'factorial_design'; ... + 'util', 'print'; ... + 'stats', 'factorial_design'; ... + 'util', 'print'; ... + 'stats', 'fmri_est'; ... + 'stats', 'review'; ... + 'util', 'print'; ... + 'stats', 'fmri_est'; ... + 'stats', 'review'; ... + 'util', 'print' }; + end + +end diff --git a/tests/tests_slow/tests_workflows/stats/test_bidsRFX_chain_several_nodes.m b/tests/tests_slow/tests_workflows/stats/test_bidsRFX_chain_several_nodes.m new file mode 100644 index 000000000..39bfe93f1 --- /dev/null +++ b/tests/tests_slow/tests_workflows/stats/test_bidsRFX_chain_several_nodes.m @@ -0,0 +1,77 @@ +% (C) Copyright 2021 bidspm developers + +function test_suite = test_bidsRFX_chain_several_nodes %#ok<*STOUT> + 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_bidsRFX_no_overwrite() + + markTestAs('slow'); + + opt = setOptions('vislocalizer', '', 'pipelineType', 'stats'); + opt.model.file = spm_file(opt.model.file, ... + 'basename', ... + 'model-vismotionNoOverWrite_smdl'); + opt.model.bm = BidsModel('file', opt.model.file); + + opt.ignore = {'qa'}; + + matlabbatch = bidsRFX('RFX', opt); + + % 2 simple dummy contrasts + % 1 complex + % 8 within group: 4 contrast from run level * 2 groups + % 4 between group: 4 contrast from run level * 1 group comparison + % but only the last 4 are returned + expected_nb_dsigns = 2 + 1 + 8 + 4; + summary = batchSummary(matlabbatch); + nb_designs = sum(sum(cellfun(@(x) strcmp(x, 'factorial_design'), summary))); + + assertEqual(nb_designs, 4); + + % folders = {}; + % for i = 1:numel(matlabbatch) + % if isfield(matlabbatch{i}.spm, 'stats') && ... + % isfield(matlabbatch{i}.spm.stats, 'factorial_design') + % [~, tmp] = fileparts(matlabbatch{i}.spm.stats.factorial_design.dir{1}); + % folders{end+1, 1} = tmp; + % + % end + % end + % folders + +end + +function test_bidsRFX_several_datasets_level() + + markTestAs('slow'); + + opt = setOptions('vislocalizer', '', 'pipelineType', 'stats'); + opt.model.file = spm_file(opt.model.file, ... + 'basename', ... + 'model-vislocalizerSeveralDatasetLevels_smdl'); + opt.model.bm = BidsModel('file', opt.model.file); + + matlabbatch = bidsRFX('RFX', opt); + + summary = batchSummary(matlabbatch); + nb_designs = sum(sum(cellfun(@(x) strcmp(x, 'factorial_design'), summary))); + + % 2 simple dummy contrasts and one complex + % but only the last one is returned + assertEqual(nb_designs, 1); +end + +function value = batchSummary(matlabbatch) + value = {}; + for i = 1:numel(matlabbatch) + field = fieldnames(matlabbatch{i}.spm); + value{i, 1} = field{1}; %#ok<*AGROW> + field = fieldnames(matlabbatch{i}.spm.(value{i, 1})); + value{i, 2} = field{1}; + end +end diff --git a/tests/tests_slow/tests_workflows/stats/test_bidsResults.m b/tests/tests_slow/tests_workflows/stats/test_bidsResults.m index b9a390034..2a7dfdc42 100644 --- a/tests/tests_slow/tests_workflows/stats/test_bidsResults.m +++ b/tests/tests_slow/tests_workflows/stats/test_bidsResults.m @@ -8,21 +8,6 @@ initTestSuite; end -function test_bidsResults_no_results() - - markTestAs('slow'); - - skipIfOctave('mixed-string-concat warning thrown'); - - opt = setOptions('vismotion', '', 'pipelineType', 'stats'); - opt.verbosity = 2; - - opt = rmfield(opt, 'results'); - - assertWarning(@() bidsResults(opt), 'bidsResults:noResultsAsked'); - -end - function test_bidsResults_subject_lvl() markTestAs('slow'); @@ -143,20 +128,6 @@ function test_bidsResults_filter_by_nodeName() end -function test_bidsResults_filter_by_nodeName_empty() - - markTestAs('slow'); - - opt = setOptions('vislocalizer', '', 'pipelineType', 'stats'); - - opt.results = defaultResultsStructure(); - - opt.results.nodeName = 'subject_level'; - - bidsResults(opt, 'nodeName', 'foo'); - -end - function test_bidsResults_too_many_backgrounds() markTestAs('slow'); @@ -282,21 +253,6 @@ function test_bidsResults_dataset_lvl() end -function test_bidsResults_error_missing_node() - - markTestAs('slow'); - - opt = setOptions('vismotion', '', 'pipelineType', 'stats'); - - opt.results = defaultResultsStructure(); - - opt.results.nodeName = 'egg'; - opt.results.name = {'spam'}; - - assertWarning(@()bidsResults(opt), 'Model:missingNode'); - -end - function opt = rmResultsFromModel(opt) for i = 1:numel(opt.model.bm.Nodes) opt.model.bm.Nodes{i}.Model.Software = rmfield(opt.model.bm.Nodes{i}.Model.Software, ... diff --git a/tests/tests_slow/tests_workflows/test_renameSegmentParameter.m b/tests/tests_slow/tests_workflows/test_renameSegmentParameter.m index 933a74b15..a107a9f88 100644 --- a/tests/tests_slow/tests_workflows/test_renameSegmentParameter.m +++ b/tests/tests_slow/tests_workflows/test_renameSegmentParameter.m @@ -30,7 +30,7 @@ function test_renameSegmentParameter_basic() 'index_dependencies', false); opt.dryRun = false; - opt.verbosity = 2; + opt.verbosity = 0; segParamFiles = spm_select('FPListRec', tmpDir, '^.*_T1w_seg8.mat$'); assertEqual(size(segParamFiles, 1), 1); diff --git a/tests/tests_slow/tests_workflows/test_renameUnwarpParameter.m b/tests/tests_slow/tests_workflows/test_renameUnwarpParameter.m index 998d9ee1a..02df29050 100644 --- a/tests/tests_slow/tests_workflows/test_renameUnwarpParameter.m +++ b/tests/tests_slow/tests_workflows/test_renameUnwarpParameter.m @@ -25,7 +25,7 @@ function test_renameUnwarpParameter_basic() 'index_dependencies', false); opt.dryRun = false; - opt.verbosity = 2; + opt.verbosity = 0; unwarpParamFiles = spm_select('FPListRec', tmpDir, '^.*_bold_uw.mat$'); assertEqual(size(unwarpParamFiles, 1), 1); diff --git a/tests/tests_stats/subject_level/test_specifyContrasts.m b/tests/tests_stats/subject_level/test_specifyContrasts.m index e6b2e4be8..72856d6fc 100644 --- a/tests/tests_stats/subject_level/test_specifyContrasts.m +++ b/tests/tests_stats/subject_level/test_specifyContrasts.m @@ -9,6 +9,9 @@ end function test_specifyContrasts_bug_854() + % no error when no contrast to build + + skipIfOctave('mixed-string-concat warning thrown'); % GIVEN subLabel = '01'; @@ -29,8 +32,8 @@ function test_specifyContrasts_bug_854() spmMatFile = cellstr(fullfile(ffxDir, 'SPM.mat')); load(spmMatFile{1}, 'SPM'); - % WHEN - contrasts = specifyContrasts(opt.model.bm, SPM); + assertWarning(@()specifyContrasts(opt.model.bm, SPM), ... + 'specifyContrasts:noContrast'); end @@ -62,6 +65,80 @@ function test_specifyContrasts_bug_815() end +function test_specifyContrasts_subject_level() + + taskName = 'motion'; + + % GIVEN + DummyContrasts{1} = 'motion'; + DummyContrasts{2} = 'static'; + + Contrasts.Name = 'motion_gt_static'; + Contrasts.ConditionList = {'motion', 'static'}; + Contrasts.Weights = [1, -1]; + Contrasts.Test = 't'; + + model = BidsModel('init', true); + + model.Input.task = taskName; + + model.Nodes{1, 1}.DummyContrasts.Contrasts = DummyContrasts; + model.Nodes{1, 1}.Contrasts{1} = Contrasts; + model.Nodes{1, 1}.GroupBy = {'run', 'subject'}; + model.Nodes{1, 1}.Model.HRF.Variables = {'motion', 'static'}; + + model.Nodes{2, 1} = BidsModel.empty_node('subject'); + model.Nodes{2, 1}.GroupBy = {'subject', 'contrast'}; + model.Nodes{2, 1}.DummyContrasts = struct('Test', 't'); + model.Nodes{2, 1} = rmfield(model.Nodes{2}, 'Contrasts'); + model.Nodes{2, 1}.Model.X = 1; + + SPM.Sess(1).col = [1, 2]; + % skip Sess 2 to make sure contrast naming is based on the Sess number + SPM.Sess(3).col = [3, 4]; + SPM.Sess(4).col = [5, 6]; + + SPM.xX.name = { ... + ' motion*bf(1)' + ' static*bf(1)' + ' motion*bf(1)' + ' static*bf(1)' + ' motion*bf(1)' + ' static*bf(1)' + }; + + SPM.xX.X = ones(1, numel(SPM.xX.name)); + + % WHEN + contrasts = specifyContrasts(model, SPM); + + % THEN + names_contrast = { ... + 'motion_1', [1 0 0 0 0 0] + 'motion_3', [0 0 1 0 0 0] + 'motion_4', [0 0 0 0 1 0] + 'static_1', [0 1 0 0 0 0] + 'static_3', [0 0 0 1 0 0] + 'static_4', [0 0 0 0 0 1] + 'motion_gt_static_1', [1 -1 0 0 0 0] + 'motion_gt_static_3', [0 0 1 -1 0 0] + 'motion_gt_static_4', [0 0 0 0 1 -1] + 'motion', [1 0 1 0 1 0] + 'static', [0 1 0 1 0 1] + 'motion_gt_static', [1 -1 1 -1 1 -1] + }; + + assertEqual(numel(contrasts), size(names_contrast, 1)); + + for i = 1:size(names_contrast, 1) + expected(i).name = names_contrast{i, 1}; + expected(i).C = names_contrast{i, 2}; + expected(i).type = 't'; + assertEqual(contrasts(i), expected(i)); + end + +end + function test_specifyContrasts_subject_level_select_node() taskName = 'motion'; @@ -132,6 +209,7 @@ function test_specifyContrasts_run_level_dummy_contrast_from_X() model.Input.task = taskName; model.Nodes = []; model.Nodes{1}.Model.X = {'motion', 'static'}; + model.Nodes{1}.Model.HRF.Variables = {'motion', 'static'}; model.Nodes{1}.Name = 'run_level'; model.Nodes{1}.Level = 'Run'; model.Nodes{1}.DummyContrasts = struct('Test', 't'); @@ -188,8 +266,12 @@ function test_specifyContrasts_missing_condition_for_dummy_contrasts() model.Input.task = taskName; model.Nodes = []; model.Nodes{1}.Level = 'Run'; + model.Nodes{1}.Name = 'Run'; + model.Nodes{1}.GroupBy = {'run', 'subject'}; model.Nodes{1}.DummyContrasts.Contrasts = DummyContrasts; model.Nodes{1}.DummyContrasts.Test = 't'; + model.Nodes{1}.Model = struct('X', {{'foo', 'bar'}}, ... + 'HRF', struct('Variables', {{'foo', 'bar'}})); SPM.Sess(1).col = [1, 2]; @@ -247,75 +329,6 @@ function test_specifyContrasts_missing_condition() end -function test_specifyContrasts_subject_level() - - taskName = 'motion'; - - % GIVEN - DummyContrasts{1} = 'motion'; - DummyContrasts{2} = 'static'; - - Contrasts.Name = 'motion_gt_static'; - Contrasts.ConditionList = {'motion', 'static'}; - Contrasts.Weights = [1, -1]; - Contrasts.Test = 't'; - - model = BidsModel('init', true); - model.Nodes{2, 1} = BidsModel.empty_node('subject'); - model.Input.task = taskName; - model.Nodes{1, 1}.DummyContrasts.Contrasts = DummyContrasts; - model.Nodes{1, 1}.Contrasts{1} = Contrasts; - model.Nodes{2, 1}.GroupBy = {'subject', 'contrast'}; - model.Nodes{2, 1}.DummyContrasts = struct('Test', 't'); - model.Nodes{2, 1} = rmfield(model.Nodes{2}, 'Contrasts'); - model.Nodes{2, 1}.Model.X = 1; - - SPM.Sess(1).col = [1, 2]; - % skip Sess 2 to make sure contrast naming is based on the Sess number - SPM.Sess(3).col = [3, 4]; - SPM.Sess(4).col = [5, 6]; - - SPM.xX.name = { ... - ' motion*bf(1)' - ' static*bf(1)' - ' motion*bf(1)' - ' static*bf(1)' - ' motion*bf(1)' - ' static*bf(1)' - }; - - SPM.xX.X = ones(1, numel(SPM.xX.name)); - - % WHEN - contrasts = specifyContrasts(model, SPM); - - % THEN - names_contrast = { ... - 'motion_1', [1 0 0 0 0 0] - 'motion_3', [0 0 1 0 0 0] - 'motion_4', [0 0 0 0 1 0] - 'static_1', [0 1 0 0 0 0] - 'static_3', [0 0 0 1 0 0] - 'static_4', [0 0 0 0 0 1] - 'motion_gt_static_1', [1 -1 0 0 0 0] - 'motion_gt_static_3', [0 0 1 -1 0 0] - 'motion_gt_static_4', [0 0 0 0 1 -1] - 'motion', [1 0 1 0 1 0] - 'static', [0 1 0 1 0 1] - 'motion_gt_static', [1 -1 1 -1 1 -1] - }; - - assertEqual(numel(contrasts), size(names_contrast, 1)); - - for i = 1:size(names_contrast, 1) - expected(i).name = names_contrast{i, 1}; - expected(i).C = names_contrast{i, 2}; - expected(i).type = 't'; - assertEqual(contrasts(i), expected(i)); - end - -end - function test_specifyContrasts_complex() % % to test the generation of contrasts when there are several runs @@ -338,6 +351,10 @@ function test_specifyContrasts_complex() model.Nodes{1}.DummyContrasts.Contrasts = DummyContrasts; model.Nodes{1}.Contrasts{1} = Contrasts; model.Nodes{1}.Level = 'Run'; + model.Nodes{1}.Name = 'Run'; + model.Nodes{1}.Model = struct('X', {{'motion', 'static'}}, ... + 'HRF', struct('Variables', {{'motion', 'static'}})); + model.Nodes{1}.GroupBy = {'run', 'subject'}; SPM.Sess(1).col = [1, 2]; % skip Sess 2 to make sure contrast naming is based on the Sess number diff --git a/tests/tests_stats/subject_level/test_specifyDummyContrasts_Session.m b/tests/tests_stats/subject_level/test_specifyDummyContrasts_Session.m index a21717b5b..3d622003d 100644 --- a/tests/tests_stats/subject_level/test_specifyDummyContrasts_Session.m +++ b/tests/tests_stats/subject_level/test_specifyDummyContrasts_Session.m @@ -58,6 +58,9 @@ function test_specifyDummyContrasts_session() model.SPM = SPM; model.Nodes{1}.DummyContrasts.Contrasts{1} = 'sign_Stim1F'; + model.Nodes{1}.GroupBy = {'run', 'subject'}; + model.Nodes{1}.Model.X = {'sign_Stim1F', 'sign_Stim2F', 'no_resp'}; + model.Nodes{1}.Model.HRF.Variables = {'sign_Stim1F', 'sign_Stim2F', 'no_resp'}; model.Nodes{2, 1} = model.Nodes{1}; model.Nodes{end}.Name = 'Session'; diff --git a/tests/tests_stats/unit_tests/test_getInclusiveMask.m b/tests/tests_stats/unit_tests/test_getInclusiveMask.m index 263d23a70..d59dc909f 100644 --- a/tests/tests_stats/unit_tests/test_getInclusiveMask.m +++ b/tests/tests_stats/unit_tests/test_getInclusiveMask.m @@ -15,7 +15,7 @@ function test_getInclusiveMask_too_many() opt = setOptions('vismotion', '01', 'useTempDir', true); - opt.verbosity = 2; + opt.verbosity = 1; opt.model.bm = BidsModel('file', opt.model.file); @@ -63,7 +63,7 @@ function test_getInclusiveMask_no_image() opt = setOptions('vismotion', '01'); - opt.verbosity = 2; + opt.verbosity = 1; opt.model.bm = BidsModel('file', opt.model.file); diff --git a/tests/tests_unit/test_warnings.m b/tests/tests_unit/test_warnings.m index 22773c784..5acaa5a77 100644 --- a/tests/tests_unit/test_warnings.m +++ b/tests/tests_unit/test_warnings.m @@ -10,31 +10,6 @@ initTestSuite; end -function test_noSPMmat() - - opt.verbosity = 0; - subLabel = '01'; - - % GIVEN - spmMatFile = fullfile(getTestDataDir(), 'mat_files', 'SPM.mat'); - % WHEN - status = noSPMmat(opt, subLabel, spmMatFile); - % THEN - assertEqual(status, false); - - % GIVEN - spmMatFile = fullfile(pwd, 'sub-01', 'stats', 'foo', 'SPM.mat'); - % WHEN - status = noSPMmat(opt, subLabel, spmMatFile); - % THEN - assertEqual(status, true); - - skipIfOctave('mixed-string-concat warning thrown'); - opt.verbosity = 1; - assertWarning(@()noSPMmat(opt, subLabel, spmMatFile), 'noSPMmat:noSpecifiedModel'); - -end - function test_noRoiFound() % GIVEN diff --git a/tests/tests_workflows/stats/test_bidsRFX_3_groups.m b/tests/tests_workflows/stats/test_bidsRFX_3_groups.m new file mode 100644 index 000000000..613f966d5 --- /dev/null +++ b/tests/tests_workflows/stats/test_bidsRFX_3_groups.m @@ -0,0 +1,64 @@ +% (C) Copyright 2024 bidspm developers + +function test_suite = test_bidsRFX_3_groups %#ok<*STOUT> + 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_bidsRFX_contrast() + + opt = setOptions('vislocalizer', '', 'pipelineType', 'stats'); + + matlabbatch = bidsRFX('contrasts', opt); + + assertEqual(numel(matlabbatch), 4); + + assertEqual(matlabbatch{1}.spm.stats.con.consess{1}.tcon.name, 'VisMot'); + assertEqual(matlabbatch{2}.spm.stats.con.consess{1}.tcon.name, 'VisStat'); + assertEqual(matlabbatch{3}.spm.stats.con.consess{1}.tcon.name, 'VisMot_&_VisStat'); + assertEqual(matlabbatch{4}.spm.stats.con.consess{1}.tcon.name, 'VisMot_&_VisStat_lt_baseline'); + +end + +function test_bidsRFX_one_way_anova_contrast() + + opt = setOptions('3_groups', '', 'pipelineType', 'stats'); + + [~, opt] = getData(opt, opt.dir.preproc); + + opt.model.bm = BidsModel('file', opt.model.file); + opt.verbosity = 0; + + nodeName = 'between_groups'; + + matlabbatch = bidsRFX('contrasts', opt, 'nodeName', nodeName); + + contrastName = 'VisMot_gt_VisStat'; + rfxDir = getRFXdir(opt, nodeName, contrastName, '1WayANOVA'); + assertEqual(fileparts(matlabbatch{1}.spm.stats.con.spmmat{1}), rfxDir); + + assertEqual(numel(matlabbatch{1}.spm.stats.con.consess), 3); + + assertEqual(matlabbatch{1}.spm.stats.con.consess{1}.tcon.name, 'relative_gt_ctrl'); + assertEqual(matlabbatch{1}.spm.stats.con.consess{1}.tcon.convec, [0; -1; 1]); + assertEqual(matlabbatch{1}.spm.stats.con.consess{2}.tcon.name, 'blind_gt_relative'); + assertEqual(matlabbatch{1}.spm.stats.con.consess{2}.tcon.convec, [-1; 0; 1]); + assertEqual(matlabbatch{1}.spm.stats.con.consess{3}.fcon.name, 'F_test'); + assertEqual(matlabbatch{1}.spm.stats.con.consess{3}.fcon.convec, [1, 0, 0; 0, 0, 1]); + + contrastName = 'VisStat_gt_VisMot'; + rfxDir = getRFXdir(opt, nodeName, contrastName, '1WayANOVA'); + assertEqual(fileparts(matlabbatch{2}.spm.stats.con.spmmat{1}), rfxDir); + + assertEqual(numel(matlabbatch{2}.spm.stats.con.consess), 3); + assertEqual(matlabbatch{2}.spm.stats.con.consess{1}.tcon.name, 'relative_gt_ctrl'); + assertEqual(matlabbatch{2}.spm.stats.con.consess{1}.tcon.convec, [0; -1; 1]); + assertEqual(matlabbatch{2}.spm.stats.con.consess{2}.tcon.name, 'blind_gt_relative'); + assertEqual(matlabbatch{2}.spm.stats.con.consess{2}.tcon.convec, [-1; 0; 1]); + assertEqual(matlabbatch{1}.spm.stats.con.consess{3}.fcon.name, 'F_test'); + assertEqual(matlabbatch{1}.spm.stats.con.consess{3}.fcon.convec, [1, 0, 0; 0, 0, 1]); + +end diff --git a/tests/tests_workflows/stats/test_bidsResults.m b/tests/tests_workflows/stats/test_bidsResults.m new file mode 100644 index 000000000..d9edae9c9 --- /dev/null +++ b/tests/tests_workflows/stats/test_bidsResults.m @@ -0,0 +1,88 @@ +% (C) Copyright 2024 bidspm developers + +function test_suite = test_bidsResults %#ok<*STOUT> + 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_bidsResults_filter_by_nodeName_empty() + + opt = setOptions('vislocalizer', '', 'pipelineType', 'stats'); + + opt.results = defaultResultsStructure(); + + opt.results.nodeName = 'subject_level'; + + bidsResults(opt, 'nodeName', 'foo'); + +end + +function test_bidsResults_error_missing_node() + + opt = setOptions('vismotion', '', 'pipelineType', 'stats'); + + opt.results = defaultResultsStructure(); + + opt.results.nodeName = 'egg'; + opt.results.name = {'spam'}; + + assertWarning(@()bidsResults(opt), 'Model:missingNode'); + +end + +function test_bidsResults_no_results() + + skipIfOctave('mixed-string-concat warning thrown'); + + opt = setOptions('vismotion', '', 'pipelineType', 'stats'); + opt.verbosity = 1; + + opt = rmfield(opt, 'results'); + + assertWarning(@() bidsResults(opt), 'bidsResults:noResultsAsked'); + +end + +function test_bidsResults_one_way_anova_results() + + opt = setOptions('3_groups', '', 'pipelineType', 'stats'); + + opt.model.bm = BidsModel('file', opt.model.file); + opt.verbosity = 0; + + nodeName = 'between_groups'; + + contrastName = 'VisMot_gt_VisStat'; + rfxDir = getRFXdir(opt, nodeName, contrastName, '1WayANOVA'); + spm_mkdir(rfxDir); + xCon(1) = struct('name', 'relative_gt_ctrl'); + xCon(2) = struct('name', 'blind_gt_relative'); + SPM = struct('nscan', 10, 'xCon', xCon); + save(fullfile(rfxDir, 'SPM.mat'), 'SPM'); + + contrastName = 'VisStat_gt_VisMot'; + rfxDir = getRFXdir(opt, nodeName, contrastName, '1WayANOVA'); + spm_mkdir(rfxDir); + % here we switch the order of the contrast in the SPM.mat + % to make sure the batch setting picks up the correct ones. + xCon(2) = struct('name', 'relative_gt_ctrl'); + xCon(1) = struct('name', 'blind_gt_relative'); + SPM = struct('nscan', 10, 'xCon', xCon); + save(fullfile(rfxDir, 'SPM.mat'), 'SPM'); + + matlabbatch = bidsResults(opt, 'nodeName', nodeName); + + assertEqual(length(matlabbatch), 4); + assertEqual(matlabbatch{1}.result.name, 'VisMotGtVisStat - relative_gt_ctrl'); + assertEqual(matlabbatch{1}.spm.stats.results.conspec.contrasts, 1); + assertEqual(matlabbatch{2}.result.name, 'VisStatGtVisMot - relative_gt_ctrl'); + assertEqual(matlabbatch{2}.spm.stats.results.conspec.contrasts, 2); + assertEqual(matlabbatch{3}.result.name, 'VisMotGtVisStat - blind_gt_relative'); + assertEqual(matlabbatch{3}.spm.stats.results.conspec.contrasts, 2); + assertEqual(matlabbatch{4}.result.name, 'VisStatGtVisMot - blind_gt_relative'); + assertEqual(matlabbatch{4}.spm.stats.results.conspec.contrasts, 1); + +end diff --git a/tests/utils.py b/tests/utils.py new file mode 100644 index 000000000..db6b502f5 --- /dev/null +++ b/tests/utils.py @@ -0,0 +1,107 @@ +import json +import shutil +from pathlib import Path + +import pandas as pd + +ROOT_DIR = Path(__file__).parent.parent + + +def touch(path: Path): + path.parent.mkdir(parents=True, exist_ok=True) + with open(path, "a"): + pass + + +def write_ds_desc(bids_dir: Path, DatasetType="raw"): + content = { + "BIDSVersion": "1.9.0", + "DatasetType": DatasetType, + "Name": f"bidspm {DatasetType}", + } + bids_dir.mkdir(parents=True, exist_ok=True) + with open(bids_dir / "dataset_description.json", "w") as f: + json.dump(content, f, indent=4) + + +def create_events_tsv(filename, onset, duration, trial_type): + filename.parent.mkdir(parents=True, exist_ok=True) + events = { + "onset": onset, + "duration": duration, + "trial_type": trial_type, + } + df = pd.DataFrame(events) + df.to_csv(filename, sep="\t", index=False) + + +def create_raw_func_vismotion(target_dir, sub, ses): + suffix = "bold" + task = "vismotion" + this_dir = target_dir / f"sub-{sub}" / f"ses-{ses}" / "func" + + basename = f"sub-{sub}_ses-{ses}_task-{task}" + + create_events_tsv( + filename=this_dir / f"{basename}_run-1_events.tsv", + onset=[2, 4], + duration=[2, 2], + trial_type=["VisMot", "VisStat"], + ) + + for run in range(1, 3): + filename = this_dir / f"{basename}_run-{run}_part-mag_{suffix}.nii" + touch(filename) + filename = this_dir / f"{basename}_run-{run}_part-phase_{suffix}.nii" + touch(filename) + + create_events_tsv( + filename=this_dir / f"{basename}_run-2_events.tsv", + onset=[3, 6], + duration=[2, 2], + trial_type=["VisStat", "VisMot"], + ) + + touch(this_dir / f"{basename}_acq-1p60mm_run-1_{suffix}.nii") + + for run in [1, 2]: + create_events_tsv( + filename=this_dir / f"{basename}_acq-1p60mm_run-{run}_events.tsv", + onset=[4, 8], + duration=[2, 2], + trial_type=["VisMot", "VisStat"], + ) + touch(this_dir / f"{basename}_acq-1p60mm_dir-PA_run-{run}_{suffix}.nii") + + +def create_preproc_func_vismotion(target_dir, sub, ses): + suffix = "bold" + task = "vismotion" + this_dir = target_dir / f"sub-{sub}" / f"ses-{ses}" / "func" + + for acq_entity in ["", "_acq-1p60mm"]: + basename = f"sub-{sub}_ses-{ses}_task-{task}{acq_entity}_part-mag" + + if ses == "01": + touch(this_dir / f"{basename}_space-individual_desc-mean_{suffix}.nii") + touch(this_dir / f"{basename}_space-IXI549Space_desc-mean_{suffix}.nii") + touch(this_dir / f"mean_{basename}_{suffix}.nii") + + for run in range(1, 3): + basename = f"sub-{sub}_ses-{ses}_task-{task}{acq_entity}_run-{run}_part-mag" + + desc_label_list = ["preproc", "mean", "smth6"] + for desc in desc_label_list: + touch(this_dir / f"{basename}_space-individual_desc-{desc}_{suffix}.nii") + touch(this_dir / f"{basename}_space-IXI549Space_desc-{desc}_{suffix}.nii") + + touch(this_dir / f"{basename}_space-individual_desc-stc_{suffix}.nii") + touch(this_dir / f"{basename}_space-IXI549Space_desc-brain_mask.nii") + + filename = this_dir / f"rp_{basename}_{suffix}.txt" + shutil.copy(ROOT_DIR / "tests" / "data" / "tsv_files" / "rp.txt", filename) + + shutil.copy( + ROOT_DIR / "tests" / "data" / "tsv_files" / "rp.tsv", + this_dir / f"{basename}_desc-confounds_regressors.tsv", + ) diff --git a/tests/utils/createDummyData.m b/tests/utils/createDummyData.m index 0091d50a9..01d6d14df 100644 --- a/tests/utils/createDummyData.m +++ b/tests/utils/createDummyData.m @@ -16,6 +16,7 @@ function createDummyData() [status, result] = system('rm -fr data/derivatives/bidspm-*/jobs'); system('python create_dummy_dataset.py'); + system('python create_3_groups_dataset.py'); cd(startDir); diff --git a/tests/utils/setOptions.m b/tests/utils/setOptions.m index ca6c02fbe..543f98a01 100644 --- a/tests/utils/setOptions.m +++ b/tests/utils/setOptions.m @@ -110,6 +110,18 @@ opt.model.file = fullfile(getTestDataDir(), 'models', ... ['model-' strjoin(task, '') '_smdl.json']); + elseif strcmp(task, '3_groups') + + task = {'vismotion'}; + + dataDir = fullfile(getTestDir(), 'data', '3_groups'); + + opt.dir.raw = fullfile(dataDir, 'bidspm-raw'); + opt.dir.preproc = fullfile(dataDir, 'derivatives', 'bidspm-preproc'); + opt.dir.stats = fullfile(dataDir, 'derivatives', 'bidspm-stats'); + + opt.model.file = fullfile(getTestDataDir(), 'models', 'model-vismotionOneWayANOVA_smdl.json'); + else opt.dir.raw = getTestDataDir('raw');