Merge pull request #191 from adelavega/enh/fixed-effects

Add fixed effects - FEMA contrasts

Author: adelavega

examples/models/ds000003/models/model-001_smdl.json

@@ -55,6 +55,12 @@
         }
       ]
     },
+    {
+      "Level": "subject",
+      "DummyContrasts": {
+        "Type": "FEMA"
+      }
+    },
     {
       "Level": "dataset",
       "DummyContrasts": {

fitlins/interfaces/nistats.py

@@ -135,7 +135,7 @@ def _run_interface(self, runtime):
         out_ents = self.inputs.contrast_info[0]['entities']
         fname_fmt = os.path.join(runtime.cwd, '{}_{}.nii.gz').format
         for name, weights, contrast_type in prepare_contrasts(
-                self.inputs.contrast_info, mat.columns.tolist()):
+              self.inputs.contrast_info, mat.columns):
             contrast_metadata.append(
                 {'contrast': name,
                  'stat': contrast_type,
@@ -178,12 +178,12 @@ def _match(query, metadata):
 class SecondLevelModel(NistatsBaseInterface, SecondLevelEstimatorInterface, SimpleInterface):
     def _run_interface(self, runtime):
         from nistats import second_level_model as level2
+        from nistats.contrasts import compute_fixed_effects
+
         smoothing_fwhm = self.inputs.smoothing_fwhm
         if not isdefined(smoothing_fwhm):
             smoothing_fwhm = None
 
-        model = level2.SecondLevelModel(smoothing_fwhm=smoothing_fwhm)
-
         effect_maps = []
         variance_maps = []
         stat_maps = []
@@ -196,46 +196,76 @@ def _run_interface(self, runtime):
         # Only keep files which match all entities for contrast
         stat_metadata = _flatten(self.inputs.stat_metadata)
         input_effects = _flatten(self.inputs.effect_maps)
+        input_variances = _flatten(self.inputs.variance_maps)
 
         filtered_effects = []
+        filtered_variances = []
         names = []
-        for m, eff in zip(stat_metadata, input_effects):
+        for m, eff, var in zip(stat_metadata, input_effects, input_variances):
             if _match(out_ents, m):
                 filtered_effects.append(eff)
+                filtered_variances.append(var)
                 names.append(m['contrast'])
 
-        # Dummy code contrast of input effects
-        design_matrix = pd.get_dummies(names)
+        mat = pd.get_dummies(names)
+        contrasts = prepare_contrasts(self.inputs.contrast_info, mat.columns)
 
-        # Fit single model for all inputs
-        model.fit(filtered_effects, design_matrix=design_matrix)
+        # Only fit model if any non-FEMA contrasts at this level
+        if any(c[2] != 'FEMA' for c in contrasts):
+            if len(filtered_effects) < 2:
+                raise RuntimeError(
+                    "At least two inputs are required for a 't' for 'F' "
+                    "second level contrast")
+            model = level2.SecondLevelModel(smoothing_fwhm=smoothing_fwhm)
+            model.fit(filtered_effects, design_matrix=mat)
 
-        for name, weights, contrast_type in prepare_contrasts(
-          self.inputs.contrast_info, design_matrix.columns.to_list()):
+        for name, weights, contrast_type in contrasts:
             contrast_metadata.append(
                 {'contrast': name,
                  'stat': contrast_type,
                  **out_ents})
 
-            maps = model.compute_contrast(
-                second_level_contrast=weights,
-                second_level_stat_type=contrast_type,
-                output_type='all')
+            # Pass-through happens automatically as it can handle 1 input
+            if contrast_type == 'FEMA':
+                # Filter effects and variances based on weights
+                ix = weights[0].astype(bool)
+
+                ffx_res = compute_fixed_effects(
+                    np.array(filtered_effects)[ix],
+                    np.array(filtered_variances)[ix]
+                    )
+
+                maps = {
+                    'effect_size': ffx_res[0],
+                    'effect_variance': ffx_res[1],
+                    'stat': ffx_res[2]
+                    }
+            else:
+                maps = model.compute_contrast(
+                    second_level_contrast=weights,
+                    second_level_stat_type=contrast_type,
+                    output_type='all'
+                    )
 
             for map_type, map_list in (('effect_size', effect_maps),
                                        ('effect_variance', variance_maps),
                                        ('z_score', zscore_maps),
                                        ('p_value', pvalue_maps),
                                        ('stat', stat_maps)):
-                fname = fname_fmt(name, map_type)
-                maps[map_type].to_filename(fname)
-                map_list.append(fname)
+                if map_type in maps:
+                    fname = fname_fmt(name, map_type)
+                    maps[map_type].to_filename(fname)
+                    map_list.append(fname)
 
         self._results['effect_maps'] = effect_maps
         self._results['variance_maps'] = variance_maps
         self._results['stat_maps'] = stat_maps
-        self._results['zscore_maps'] = zscore_maps
-        self._results['pvalue_maps'] = pvalue_maps
         self._results['contrast_metadata'] = contrast_metadata
 
+        # These are "optional" as fixed effects do not support these
+        if zscore_maps:
+            self._results['zscore_maps'] = zscore_maps
+        if pvalue_maps:
+            self._results['pvalue_maps'] = pvalue_maps
+
         return runtime

fitlins/interfaces/utils.py

@@ -32,9 +32,11 @@ def _list_outputs(self):
         outputs = self._outputs().get()
         for key in self._fields:
             val = getattr(self.inputs, key)
-            if self._check_lengths is True:
-                self._calculate_length(val)
-            outputs[key] = [elem for sublist in val for elem in sublist]
+            # Allow for empty inputs
+            if isdefined(val):
+                if self._check_lengths is True:
+                    self._calculate_length(val)
+                outputs[key] = [elem for sublist in val for elem in sublist]
         self._lengths = None
 
         return outputs
@@ -72,12 +74,14 @@ def _run_interface(self, runtime):
         self._results.update({'metadata': [], 'out': []})
         for key in self._fields:
             val = getattr(self.inputs, key)
-            if len(val) != n:
-                raise ValueError(f"List lengths must match metadata. Failing list: {key}")
-            for md, obj in zip(orig_metadata, val):
-                metadata = md.copy()
-                metadata.update(md_map.get(key, {}))
-                self._results['metadata'].append(metadata)
-                self._results['out'].append(obj)
+            # Allow for missing values
+            if isdefined(val):
+                if len(val) != n:
+                    raise ValueError(f"List lengths must match metadata. Failing list: {key}")
+                for md, obj in zip(orig_metadata, val):
+                    metadata = md.copy()
+                    metadata.update(md_map.get(key, {}))
+                    self._results['metadata'].append(metadata)
+                    self._results['out'].append(obj)
 
         return runtime

requirements.txt

@@ -1 +1 @@
-nistats @ git+https://github.com/nistats/nistats.git@b69e127020df05aa8ba540f9a8d799e1dad3d602
+nistats @ git+https://github.com/nistats/nistats.git@e88a5f

setup.cfg

@@ -28,8 +28,8 @@ install_requires =
     nilearn>=0.4
     pandas>=0.19
     tables>=3.2.1
-    nistats>=0.0.1b0
-    pybids @ git+https://github.com/bids-standard/pybids.git@24d2032020fd407e5bcbbbf20943726f48717574
+    nistats@git+https://github.com/nistats/nistats.git@e88a5f
+    pybids>=0.10.0
     jinja2
 
 [options.extras_require]