Meta fix (addresses issue #290)

fitlins/interfaces/nistats.py

@@ -3,6 +3,8 @@
 import pandas as pd
 from functools import partial
 
+
+
 from nipype.interfaces.base import LibraryBaseInterface, SimpleInterface, isdefined
 
 from .abstract import (
@@ -322,7 +324,7 @@ def _run_interface(self, runtime):
             compute_fixed_effects,
             _compute_fixed_effects_params,
         )
-
+        from ..interfaces.pymare_extension import pymare_model
         spec = self.inputs.spec
         smoothing_fwhm = self.inputs.smoothing_fwhm
         smoothing_type = self.inputs.smoothing_type
@@ -332,7 +334,6 @@ def _run_interface(self, runtime):
             raise NotImplementedError(
                 "Only the iso smoothing type is available for the nistats estimator."
             )
-
         effect_maps = []
         variance_maps = []
         stat_maps = []
@@ -342,12 +343,11 @@ def _run_interface(self, runtime):
         spec_metadata = spec['metadata'].to_dict('records')
         out_ents = spec['entities'].copy()  # Same for all
         out_ents.setdefault("contrast", spec['contrasts'][0]['name'])
-
+    
         # 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 = []
@@ -361,17 +361,15 @@ def _run_interface(self, runtime):
                     names.append(m['contrast'])
 
         contrasts = prepare_contrasts(spec['contrasts'], spec['X'].columns)
-
         is_cifti = filtered_effects[0].endswith('dscalar.nii')
         if is_cifti:
             fname_fmt = os.path.join(runtime.cwd, '{}_{}.dscalar.nii').format
         else:
             fname_fmt = os.path.join(runtime.cwd, '{}_{}.nii.gz').format
 
         model_type = spec["model"].get("type", "")
-
-        # Do not fit model for meta-analyses
-        if model_type != 'Meta':
+
+        if model_type.lower() != "meta":
             if len(filtered_effects) < 2:
                 raise RuntimeError(
                     "At least two inputs are required for a 't' for 'F' " "second level contrast"
@@ -386,6 +384,13 @@ def _run_interface(self, runtime):
             else:
                 model = level2.SecondLevelModel(smoothing_fwhm=smoothing_fwhm)
                 model.fit(filtered_effects, design_matrix=spec['X'])
+        else:
+            if is_cifti:
+                model = pymare_model(is_cifti=True)
+                model.fit(filtered_effects, filtered_variances, spec['X'])
+            else:
+                model = pymare_model(is_cifti=False)
+                model.fit(filtered_effects, filtered_variances, spec['X'])
 
         for name, weights, cont_ents, contrast_test in contrasts:
             contrast_metadata.append(
@@ -397,40 +402,7 @@ def _run_interface(self, runtime):
                 }
             )
 
-            # Pass-through happens automatically as it can handle 1 input
-            if model_type == 'Meta':
-                # Index design identity matrix on non-zero contrasts weights
-                con_ix = weights[0].astype(bool)
-                # Index of all input files "involved" with that contrast
-                dm_ix = spec['X'].iloc[:, con_ix].any(axis=1)
-
-                contrast_imgs = np.array(filtered_effects)[dm_ix]
-                variance_imgs = np.array(filtered_variances)[dm_ix]
-                if is_cifti:
-                    ffx_cont, ffx_var, ffx_t = _compute_fixed_effects_params(
-                        np.squeeze(
-                            [nb.load(fname).get_fdata(dtype='f4') for fname in contrast_imgs]
-                        ),
-                        np.squeeze(
-                            [nb.load(fname).get_fdata(dtype='f4') for fname in variance_imgs]
-                        ),
-                        precision_weighted=False,
-                    )
-                    img = nb.load(filtered_effects[0])
-                    maps = {
-                        'effect_size': dscalar_from_cifti(img, ffx_cont, "effect_size"),
-                        'effect_variance': dscalar_from_cifti(img, ffx_var, "effect_variance"),
-                        'stat': dscalar_from_cifti(img, ffx_t, "stat"),
-                    }
-
-                else:
-                    ffx_res = compute_fixed_effects(contrast_imgs, variance_imgs)
-                    maps = {
-                        'effect_size': ffx_res[0],
-                        'effect_variance': ffx_res[1],
-                        'stat': ffx_res[2],
-                    }
-            else:
+            if model_type.lower() != "meta":
                 if is_cifti:
                     contrast = compute_contrast(
                         labels, estimates, weights, contrast_type=contrast_test
@@ -452,6 +424,16 @@ def _run_interface(self, runtime):
                         second_level_stat_type=contrast_test,
                         output_type='all',
                     )
+            else:
+                if is_cifti:
+                    dict_maps = model.compute_contrast(weights)
+                    img = nb.load(filtered_effects[0])
+                    maps = {
+                        map_type: dscalar_from_cifti(img, dict_maps[map_type], map_type)
+                        for map_type in dict_maps.keys()
+                    } 
+                else:
+                    maps = model.compute_contrast(weights)
 
             for map_type, map_list in (
                 ('effect_size', effect_maps),

fitlins/interfaces/pymare_extension.py

@@ -0,0 +1,80 @@
+import nibabel as nb
+import numpy as np
+from nilearn.input_data import NiftiMasker
+from pymare import estimators
+from nilearn.image import mean_img
+
+
+def remove_zero_var_voxels_(effect_data, variance_data):
+    nonzero_var_mask = np.sum(variance_data == 0, 0) == 0
+    variance_data_masked = variance_data[:, nonzero_var_mask]
+    effect_data_masked = effect_data[:, nonzero_var_mask]
+    return effect_data_masked, variance_data_masked, nonzero_var_mask
+
+
+
+class pymare_model:
+    def __init__(
+        self,
+        smoothing_fwhm=None,
+        is_cifti=False
+    ):
+        self.smoothing_fwhm = smoothing_fwhm
+        self.is_cifti = is_cifti
+
+    def fit(
+            self,
+            filtered_effects=None,
+            filtered_variances=None,
+            design_matrix=None
+    ):
+        self.design_matrix_ = design_matrix.to_numpy()
+        if self.is_cifti:
+            effect_data = np.squeeze(
+                [nb.load(effect).get_fdata(dtype='f4') for effect in filtered_effects]
+            )
+            variance_data = np.squeeze(
+                [nb.load(variance).get_fdata(dtype='f4') for variance in filtered_variances]
+            )
+            self.effect_data_, self.variance_data_, self.nonzero_var_mask_ = \
+                remove_zero_var_voxels_(effect_data, variance_data)
+
+        else:
+            self.masker_ = NiftiMasker(
+                smoothing_fwhm=self.smoothing_fwhm,
+            )
+            sample_map = mean_img(filtered_effects)
+            self.masker_.fit(sample_map)
+            effect_data = self.masker_.transform(filtered_effects)
+            variance_data = self.masker_.transform(filtered_variances)
+            self.effect_data_, self.variance_data_, self.nonzero_var_mask_ = \
+                remove_zero_var_voxels_(effect_data, variance_data)
+
+        self.wls_ = estimators.WeightedLeastSquares()
+        self.wls_.fit(y=self.effect_data_,
+                      v=self.variance_data_,
+                      X=self.design_matrix_
+                      )
+
+    def compute_contrast(
+        self,
+        con_val=None,
+    ):
+        outputs = {}
+        tmp = np.einsum('ij, jkl->kl', con_val, self.wls_.params_['inv_cov'])
+        outputs['effect_variance'] = np.einsum('ij, jk->k', con_val, tmp)
+        outputs['effect_size'] = np.einsum('ij, jk->k', con_val, self.wls_.params_['fe_params'])
+        outputs['stat'] = outputs['effect_size'] / outputs['effect_variance']**.5
+
+        outputs_array = {}
+        for image_type, image in outputs.items():
+            outputs_array[image_type] = np.zeros(self.nonzero_var_mask_.shape)
+            outputs_array[image_type][self.nonzero_var_mask_] = image
+
+        if self.is_cifti is False:
+            outputs_nifti = {}
+            for image_type, image in outputs_array.items():
+                outputs_nifti[image_type] = self.masker_.inverse_transform(image)
+            return outputs_nifti
+        else:
+            return outputs_array

setup.cfg

@@ -31,6 +31,7 @@ install_requires =
     tables>=3.2.1
     pybids~=0.15.4
     jinja2
+    pymare~=0.0.3
 
 [options.extras_require]
 duecredit = duecredit