Compute correlation matrix in signal voxels (#889)

* Compute correlation matrix in signal voxels

* Calculate correlation in good voxels

* Test nonaggressive mask here

* Update nimare/meta/ibma.py

---------

Co-authored-by: James Kent <[email protected]>

nimare/meta/ibma.py

@@ -92,6 +92,9 @@ def _preprocess_input(self, dataset):
         if isinstance(mask_img, str):
             mask_img = nib.load(mask_img)
 
+        # Reserve the key for the correlation matrix
+        self.inputs_["corr_matrix"] = None
+
         if self.aggressive_mask:
             # Ensure that protected values are not included among _required_inputs
             assert (
@@ -123,11 +126,12 @@ def _preprocess_input(self, dataset):
                 # Mask required input images using either the dataset's mask or the estimator's.
                 temp_arr = masker.transform(img4d)
 
-                # To save memory, we only save the original image array and perform masking later.
+                # To save memory, we only save the original image array and perform masking later
+                # in the estimator if self.aggressive_mask is True.
                 self.inputs_[name] = temp_arr
+
                 if self.aggressive_mask:
-                    # An intermediate step to mask out bad voxels.
-                    # Can be dropped once PyMARE is able to handle masked arrays or missing data.
+                    # Determine the good voxels here
                     nonzero_voxels_bool = np.all(temp_arr != 0, axis=0)
                     nonnan_voxels_bool = np.all(~np.isnan(temp_arr), axis=0)
                     good_voxels_bool = np.logical_and(nonzero_voxels_bool, nonnan_voxels_bool)
@@ -148,14 +152,11 @@ def _preprocess_input(self, dataset):
 
         # Further reduce image-based inputs to remove "bad" voxels
         # (voxels with zeros or NaNs in any studies)
-        if "aggressive_mask" in self.inputs_.keys():
+        if self.aggressive_mask:
             if n_bad_voxels := (
                 self.inputs_["aggressive_mask"].size - self.inputs_["aggressive_mask"].sum()
             ):
-                LGR.warning(
-                    f"Masking out {n_bad_voxels} additional voxels. "
-                    "The updated masker is available in the Estimator.masker attribute."
-                )
+                LGR.warning(f"Masking out {n_bad_voxels} additional voxels.")
 
 
 class Fishers(IBMAEstimator):
@@ -396,6 +397,25 @@ def _preprocess_input(self, dataset):
         self.inputs_["contrast_names"] = np.array([label_to_int[label] for label in labels])
         self.inputs_["num_contrasts"] = np.array([label_counts[label] for label in labels])
 
+        n_studies = len(self.inputs_["id"])
+        if n_studies != np.unique(self.inputs_["contrast_names"]).size:
+            # If all studies are not unique, we will need to correct for multiple contrasts
+            # Calculate correlation matrix on valid voxels
+            if self.aggressive_mask:
+                voxel_mask = self.inputs_["aggressive_mask"]
+                self.inputs_["corr_matrix"] = np.corrcoef(
+                    self.inputs_["z_maps"][:, voxel_mask],
+                    rowvar=True,
+                )
+            else:
+                self.inputs_["corr_matrix"] = np.zeros((n_studies, n_studies), dtype=float)
+                for bag in self.inputs_["data_bags"]["z_maps"]:
+                    study_bag = bag["study_mask"]
+                    self.inputs_["corr_matrix"][np.ix_(study_bag, study_bag)] = np.corrcoef(
+                        bag["values"],
+                        rowvar=True,
+                    )
+
     def _generate_description(self):
         description = (
             f"An image-based meta-analysis was performed with NiMARE {__version__} "
@@ -464,17 +484,12 @@ def _fit(self, dataset):
                 "will produce invalid results."
             )
 
-        corr = None
-        if self.inputs_["contrast_names"].size != np.unique(self.inputs_["contrast_names"]).size:
-            # If all studies are not unique, we will need to correct for multiple contrasts
-            corr = np.corrcoef(self.inputs_["z_maps"], rowvar=True)
-
         if self.aggressive_mask:
             voxel_mask = self.inputs_["aggressive_mask"]
 
             result_maps = self._fit_model(
                 self.inputs_["z_maps"][:, voxel_mask],
-                corr=corr,
+                corr=self.inputs_["corr_matrix"],
             )
 
             z_map, p_map, dof_map = tuple(
@@ -490,7 +505,9 @@ def _fit(self, dataset):
                     z_map[bag["voxel_mask"]],
                     p_map[bag["voxel_mask"]],
                     dof_map[bag["voxel_mask"]],
-                ) = self._fit_model(bag["values"], bag["study_mask"], corr=corr)
+                ) = self._fit_model(
+                    bag["values"], bag["study_mask"], corr=self.inputs_["corr_matrix"]
+                )
 
         maps = {"z": z_map, "p": p_map, "dof": dof_map}
         description = self._generate_description()

nimare/reports/base.py

@@ -33,6 +33,7 @@
     from importlib_resources import files
 
 import jinja2
+import numpy as np
 import pandas as pd
 
 from nimare.meta.cbma.base import CBMAEstimator, PairwiseCBMAEstimator
@@ -50,7 +51,6 @@
     plot_mask,
     plot_static_brain,
 )
-from nimare.stats import pearson
 
 PARAMETERS_DICT = {
     "kernel_transformer__fwhm": "FWHM",
@@ -282,13 +282,6 @@ def _gen_fig_summary(img_key, threshold, out_filename):
     (out_filename).write_text(summary_text, encoding="UTF-8")
 
 
-def _compute_similarities(maps_arr, ids_):
-    """Compute the similarity between maps."""
-    corrs = [pearson(img_map, maps_arr) for img_map in list(maps_arr)]
-
-    return pd.DataFrame(index=ids_, columns=ids_, data=corrs)
-
-
 def _gen_figures(results, img_key, diag_name, threshold, fig_dir):
     """Generate html and png objects for the report."""
     # Plot brain images if not empty
@@ -527,7 +520,31 @@ def __init__(
                     self.fig_dir / f"preliminary_dset-{dset_i+1}_figure-summarystats.html",
                 )
 
-                similarity_table = _compute_similarities(maps_arr, ids_)
+                # Compute similarity matrix
+                if self.results.estimator.inputs_["corr_matrix"] is None:
+                    n_studies = len(ids_)
+                    if self.results.estimator.aggressive_mask:
+                        voxel_mask = self.results.estimator.inputs_["aggressive_mask"]
+                        corr = np.corrcoef(
+                            self.results.estimator.inputs_["z_maps"][:, voxel_mask],
+                            rowvar=True,
+                        )
+                    else:
+                        corr = np.zeros((n_studies, n_studies), dtype=float)
+                        for bag in self.results.estimator.inputs_["data_bags"]["z_maps"]:
+                            study_bag = bag["study_mask"]
+                            corr[np.ix_(study_bag, study_bag)] = np.corrcoef(
+                                bag["values"],
+                                rowvar=True,
+                            )
+                else:
+                    corr = self.inputs_["corr_matrix"]
+
+                similarity_table = pd.DataFrame(
+                    index=ids_,
+                    columns=ids_,
+                    data=corr,
+                )
 
                 plot_heatmap(
                     similarity_table,

nimare/tests/test_meta_ibma.py

@@ -219,9 +219,10 @@ def test_ibma_resampling(testdata_ibma_resample, resample_kwargs):
     assert isinstance(results, nimare.results.MetaResult)
 
 
-def test_stouffers_multiple_contrasts(testdata_ibma_multiple_contrasts):
+@pytest.mark.parametrize("aggressive_mask", [True, False], ids=["aggressive", "liberal"])
+def test_stouffers_multiple_contrasts(testdata_ibma_multiple_contrasts, aggressive_mask):
     """Test Stouffer's correction with multiple contrasts."""
-    meta = ibma.Stouffers()
+    meta = ibma.Stouffers(aggressive_mask=aggressive_mask)
     results = meta.fit(testdata_ibma_multiple_contrasts)
 
     assert isinstance(results, nimare.results.MetaResult)