[FIX] fix group level results after contrasts smoothing (#1021)

* fix warning

* fix typo

* update test

* add cumulative smoothing to report

Author: Remi-Gau

demos/openneuro/ds002799_run.m

@@ -14,6 +14,7 @@
 
 % The directory where the data are located
 root_dir = fullfile(fileparts(mfilename('fullpath')));
+root_dir = pwd;
 
 %% Parameters
 
@@ -75,6 +76,9 @@
 % update the dataset level node
 run_lvl_idx = 3;
 bm.Nodes{run_lvl_idx}.GroupBy = {'contrast'};
+bm.Nodes{run_lvl_idx}.Model.Software.bidspm.Results(1) = defaultResultsStructure();
+bm.Nodes{run_lvl_idx}.Model.Software.bidspm.Results(1).nodeName = 'dataset';
+bm.Nodes{run_lvl_idx}.Model.Software.bidspm.Results(1).name = {'es'};
 
 % write
 bm.write(model_file);

src/reports/boilerplate.m

@@ -89,17 +89,13 @@
       opt.unwarp = true;
     end
 
-    opt.smoothing = true;
-    if opt.fwhm.func == 0
-      opt.smoothing = false;
-    end
+    opt.smoothing = opt.fwhm.func > 0;
 
   elseif strcmp(pipelineType, 'stats')
 
-    opt.smoothing = true;
-
-    if opt.fwhm.contrast == 0
-      opt.smoothing = false;
+    opt.smoothing = opt.fwhm.contrast > 0;
+    if opt.smoothing
+      opt.fwhm.cumulative = computeCumulativeFwhm(opt);
     end
 
     bm = BidsModel('file', opt.model.file);

src/reports/boilerplate_stats.mustache

@@ -43,7 +43,7 @@ to censor any time points flagged as outlier,
 {{/confounds}}.
 
 {{#smoothing}}
-Contrast images were spatially smoothed using a 3D gaussian kernel (FWHM = {{fwhm.contrast}} mm).
+Contrast images were spatially smoothed using a 3D gaussian kernel (FWHM = {{fwhm.contrast}} mm resulting in a final smoothing equivalent to {{fwhm.cumulative}} mm).
 {{/smoothing}}
 {{! TODO Table of conditions with duration of each event}}
 

src/stats/group_level/computeCumulativeFwhm.m

@@ -0,0 +1,16 @@
+function fwhm = computeCumulativeFwhm(opt)
+  %
+  % Compute resulting fwhm when smoothing both time series and contrasts.
+  %
+  % USAGE::
+  %
+  %  fwhm = computeCumulativeFwhm(opt)
+  %
+
+  % (C) Copyright 2023 Remi Gau
+
+  if opt.fwhm.contrast > 0
+    fwhm = (opt.fwhm.func^2 + opt.fwhm.contrast^2)^0.5;
+  end
+
+end

src/workflows/stats/bidsResults.m

@@ -156,11 +156,9 @@
 
   opt.dir.output = opt.dir.stats;
 
-  % TODO
-  % bids stats model should override options
-  if ~isfield(opt, 'results') || isempty(opt.results) || ...
-          strcmp(opt.results(1).name{1}, '')
-    opt.results = opt.model.bm.getResults();
+  modelResults = opt.model.bm.getResults();
+  if ~isempty(modelResults)
+    opt.results = modelResults;
   end
 
   status = checks(opt);
@@ -365,7 +363,7 @@
 
     tmp.dir = getFFXdir(subLabel, opt);
 
-    status = checkSpmMat(tmp.dir);
+    status = checkSpmMat(tmp.dir, opt);
 
     if ~status
       return
@@ -486,18 +484,21 @@
 
 end
 
-function status = checkSpmMat(dir)
+function status = checkSpmMat(dir, opt)
   status = exist(fullfile(dir, 'SPM.mat'), 'file');
   if ~status
-    msg = sprintf('\nCould not find a SPM.mat file in directory %s\n', dir);
+    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);
+    logger('WARNING', msg, 'id', id, 'options', opt, 'filename', mfilename);
   end
 end
 
 function [matlabbatch, results] = appendToBatch(matlabbatch, opt, results, result)
 
-  if ~checkSpmMat(result.dir)
+  if ~checkSpmMat(result.dir, opt)
     return
   end
 

tests/tests_workflows/stats/test_bidsResults.m

@@ -22,6 +22,8 @@ function test_bidsResults_subject_lvl_regex()
 
   opt.subjects = {'01'};
 
+  opt = rmResultsFromModel(opt);
+
   %% WHEN
   matlabbatch = bidsResults(opt);
 
@@ -64,6 +66,8 @@ function test_bidsResults_filter_by_nodeName()
   opt.results(2).nodeName = 'dataset_level';
   opt.results(2).name = {'VisMot_gt_VisStat'};
 
+  opt = rmResultsFromModel(opt);
+
   %% WHEN
   matlabbatch = bidsResults(opt, 'nodeName', 'subject_level');
 
@@ -103,6 +107,8 @@ function test_bidsResults_too_many_backgrounds()
   opt.results.montage.do = true;
   opt.results.montage.background = struct('suffix', 'probseg');
 
+  opt = rmResultsFromModel(opt);
+
   if bids.internal.is_octave()
     return
   end
@@ -114,7 +120,8 @@ function test_bidsResults_too_many_backgrounds()
 function test_bidsResults_background_for_subject()
 
   %% GIVEN
-  opt = setOptions('vislocalizer', {'01', 'ctrl01', 'blind01'}, 'pipelineType', 'stats');
+  opt = setOptions('vislocalizer', {'01', 'ctrl01', 'blind01'}, ...
+                   'pipelineType', 'stats');
 
   % Specify what output we want
   opt.results = defaultResultsStructure();
@@ -132,14 +139,20 @@ function test_bidsResults_background_for_subject()
                                           'suffix', 'T1w');
 
   opt.subjects = {'^01'};
+
+  opt = rmResultsFromModel(opt);
+
   matlabbatch = bidsResults(opt);
+
   bf = bids.File(matlabbatch{1}.spm.stats.results.export{4}.montage.background{1});
   assertEqual(bf.entities.desc, 'preproc');
   assertEqual(bf.entities.sub, '01');
   assertEqual(bf.suffix, 'T1w');
 
   opt.subjects = {'ctrl01'};
+
   matlabbatch = bidsResults(opt);
+
   bf = bids.File(matlabbatch{1}.spm.stats.results.export{4}.montage.background{1});
   assertEqual(bf.entities.desc, 'preproc');
   assertEqual(bf.entities.sub, 'ctrl01');
@@ -165,6 +178,8 @@ function test_bidsResults_no_background_for_montage()
   opt.results.montage.do = true;
   opt.results.montage.background = 'aFileThatDoesNotExist.nii';
 
+  opt = rmResultsFromModel(opt);
+
   if bids.internal.is_octave()
     return
   end
@@ -190,7 +205,7 @@ function test_bidsResults_dataset_lvl()
 
   %% WHEN
 
-  % matlabbatch = bidsResults(opt);
+  matlabbatch = bidsResults(opt);
 
 end
 
@@ -222,17 +237,14 @@ function test_bidsResults_subject_lvl()
   opt.results.MC =  'FWE';
   opt.results.p = 0.05;
   opt.results.k = 5;
-
   opt.results.png = true();
-
   opt.results.csv = true();
-
   opt.results.threshSpm = true();
-
   opt.results.binary = true();
-
   opt.results.nidm = true();
 
+  opt = rmResultsFromModel(opt);
+
   %% WHEN
   matlabbatch = bidsResults(opt);
 
@@ -272,3 +284,10 @@ function test_bidsResults_subject_lvl()
               expected{4}.binary.basename);
 
 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, ...
+                                                   'bidspm');
+  end
+end