fix contrasts

cpp-lln-lab · Jul 25, 2024 · 8efb9ad · 8efb9ad
1 parent e79deeb
commit 8efb9ad
Show file tree

Hide file tree

Showing 8 changed files with 216 additions and 41 deletions.
diff --git a/demos/openneuro/ds003397_run.m b/demos/openneuro/ds003397_run.m
@@ -35,10 +35,6 @@
                       'models', ...
                       'model-ds003397_smdl.json');
 
-opt.results(1) = defaultResultsStructure();
-opt.results(1).nodeName = 'subject_level';
-opt.results(1).name = 'flashing checkerboard';
-
 % bidspm(bids_dir, output_dir, 'subject', ...
 %        'participant_label', participant_label, ...
 %        'action', 'results', ...
@@ -48,8 +44,7 @@
 %        'space', space, ...
 %        'fwhm', 0, ...
 %        'skip_validation', true, ...
-%        'verbosity', 3, ...
-%        'opt', opt);
+%        'verbosity', 3);
 
 bidspm(bids_dir, output_dir, 'dataset', ...
        'participant_label', participant_label, ...
@@ -60,5 +55,4 @@
        'space', space, ...
        'fwhm', 0, ...
        'skip_validation', true, ...
-       'verbosity', 3, ...
-       'opt', opt);
+       'verbosity', 3);
diff --git a/demos/openneuro/models/model-ds003397_smdl.json b/demos/openneuro/models/model-ds003397_smdl.json
@@ -55,7 +55,38 @@
         "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"
@@ -73,22 +104,55 @@
         "X": [
           1,
           "group"
-        ]
+        ],
+        "Software": {
+          "bidspm": {
+            "Results": [
+              {
+                "name": [
+                  "B > I",
+                  "average across groups"
+                ],
+                "p": 0.05,
+                "MC": "none",
+                "png": true,
+                "binary": false,
+                "nidm": false,
+                "montage": {
+                  "do": false
+                }
+              }
+            ]
+          }
+        }
       },
       "Contrasts": [
         {
-          "Name": "flashing checkerboard",
+          "Name": "B > I",
           "ConditionList": [
             "group.B",
-            "group.I",
-            "group.BI"
+            "group.I"
           ],
           "Weights": [
             1,
             -1,
             0
           ],
           "Test": "t"
+        },
+        {
+          "Name": "average across groups",
+          "ConditionList": [
+            "group.B",
+            "group.I",
+            "group.BI"
+          ],
+          "Weights": [
+            1,
+            1,
+            1
+          ],
+          "Test": "t"
         }
       ]
     }

diff --git a/src/batches/stats/setBatchFactorialDesign.m b/src/batches/stats/setBatchFactorialDesign.m
@@ -61,10 +61,10 @@
       contrastsList = {};
       groups = {};
 
+      % TODO all contrasts should have the same name
       contrasts = getContrastsListForFactorialDesign(opt, nodeName);
 
       % TODO refactor
-
       designMatrix = opt.model.bm.get_design_matrix('Name', nodeName);
       designMatrix = cellfun(@(x) num2str(x), designMatrix, 'uniformoutput', false);
 
@@ -73,7 +73,11 @@
 
       checkColumnParticipantsTsv(BIDS, groupColumnHdr);
 
-      availableGroups = unique(BIDS.raw.participants.content.(groupColumnHdr));
+      % 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.
+      availableGroups = sort(unique(BIDS.raw.participants.content.(groupColumnHdr)));
 
       label = '1WayANOVA';
 

diff --git a/src/batches/stats/setBatchGroupLevelContrasts.m b/src/batches/stats/setBatchGroupLevelContrasts.m
@@ -48,8 +48,9 @@
 
         end
 
+        % TODO make more general than just with group
       elseif all(ismember(lower(groupBy), {'contrast', 'group'}))
-
+        % TODO make more general than just with group
         groupColumnHdr = groupBy{ismember(lower(groupBy), {'group'})};
         availableGroups = unique(participants.(groupColumnHdr));
 
@@ -71,13 +72,66 @@
 
       end
 
+    case 'one_way_anova'
+
+      % T test for ANOVA
+      %
+      % Loop over the subject level contrasts passed
+      % through the Edge filter.
+      % Then generate the between group contrasts.
+
+      designMatrix = removeIntercept(designMatrix);
+
+      groups = unique(participants.(designMatrix{1}));
+
+      edge = opt.model.bm.get_edge('Destination', nodeName);
+      contrastsList = edge.Filter.contrast;
+
+      thisContrast = opt.model.bm.get_contrasts('Name', nodeName);
+
+      for j = 1:numel(contrastsList)
+
+        spmMatFile = fullfile(getRFXdir(opt, ...
+                                        nodeName, ...
+                                        contrastsList{j}, ...
+                                        '1WayANOVA'), ...
+                              'SPM.mat');
+
+        for iCon = 1:numel(thisContrast)
+
+          % Sort conditions and weights
+          [ConditionList, I] = sort(thisContrast{iCon}.ConditionList);
+          for iCdt = 1:numel(ConditionList)
+            ConditionList{iCdt} =  strrep(ConditionList{iCdt}, [designMatrix{1}, '.'], '');
+          end
+          Weights = thisContrast{iCon}.Weights(I);
+
+          % Create contrast vectors by what was passed in the model
+          convec = zeros(size(groups));
+          for iGroup = 1:numel(groups)
+            index = strcmp(groups{iGroup}, ConditionList);
+            if any(index)
+              convec(iGroup) = Weights(index);
+            end
+          end
+
+          consess{iCon}.tcon.name = thisContrast{iCon}.Name;
+          consess{iCon}.tcon.convec = convec;
+          consess{iCon}.tcon.sessrep = 'none';
+        end
+
+        matlabbatch = setBatchContrasts(matlabbatch, opt, spmMatFile, consess);
+
+      end
+
     case 'two_sample_t_test'
 
       designMatrix = removeIntercept(designMatrix);
 
+      % TODO make more general than just with group
       if ismember(lower(designMatrix), {'group'})
-        % TODO will this ignore the contrasts define at other levels and not
-        % passed through the filter ?
+        % 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;
       end

diff --git a/src/stats/group_level/getRFXdir.m b/src/stats/group_level/getRFXdir.m
@@ -68,6 +68,8 @@
   sub = 'ALL';
   if ~isempty(contrastName)
     thisNode = opt.model.bm.get_nodes('Name', nodeName);
+
+    % TODO make more general than just with group
     if all(ismember(lower(thisNode.GroupBy), {'contrast', 'group'})) && ~isempty(thisGroup)
       sub = thisGroup;
     end

diff --git a/src/workflows/stats/bidsResults.m b/src/workflows/stats/bidsResults.m
@@ -475,8 +475,11 @@
 
           matlabbatch = appendToBatch(matlabbatch, opt, result);
 
+          % TODO make more general than just with group
         elseif all(ismember(lower(groupBy), {'contrast', 'group'}))
 
+          % TODO make more general than just with group
+
           groupColumnHdr = groupBy{ismember(lower(groupBy), {'group'})};
           availableGroups = unique(participants.(groupColumnHdr));
 

diff --git a/tests/data/models/model-vismotionOneWayANOVA_smdl.json b/tests/data/models/model-vismotionOneWayANOVA_smdl.json
@@ -46,6 +46,18 @@
             -1
           ],
           "Test": "t"
+        },
+        {
+          "Name": "VisStat_gt_VisMot",
+          "ConditionList": [
+            "trial_type.VisMot",
+            "trial_type.VisStat"
+          ],
+          "Weights": [
+            -1,
+            1
+          ],
+          "Test": "t"
         }
       ]
     },
@@ -100,17 +112,31 @@
       "Contrasts": [
         {
           "Name": "VisMot_gt_VisStat",
+          "Description": "3 groups but not sorted alphabetically whereas the groups are entered alphabetically in the design matrix.",
           "ConditionList": [
+            "diagnostic.relative",
             "diagnostic.ctrl",
-            "diagnostic.blind",
-            "diagnostic.relative"
+            "diagnostic.blind"
           ],
           "Weights": [
             1,
             -1,
             0
           ],
           "Test": "t"
+        },
+        {
+          "Name": "VisMot_gt_VisStat",
+          "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"
         }
       ]
     }
@@ -129,7 +155,8 @@
       "Destination": "between_groups",
       "Filter": {
         "contrast": [
-          "VisMot_gt_VisStat"
+          "VisMot_gt_VisStat",
+          "VisStat_gt_VisMot"
         ]
       }
     }