Add Workflow and CBMAWorkflow classes. Support pairwise CBMA workflows

docs/api.rst

@@ -334,7 +334,9 @@ For more information about fetching data from the internet, see :ref:`fetching t
 
    workflows.ale_sleuth_workflow
    workflows.macm_workflow
-   workflows.cbma_workflow
+   workflows.base.Workflow
+   workflows.cbma.CBMAWorkflow
+   workflows.cbma.PairwiseCBMAWorkflow
 
 :mod:`nimare.reports`: NiMARE report
 --------------------------------------------------

examples/02_meta-analyses/10_plot_cbma_workflow.py

@@ -20,7 +20,7 @@
 from nimare.dataset import Dataset
 from nimare.reports.base import run_reports
 from nimare.utils import get_resource_path
-from nimare.workflows import cbma_workflow
+from nimare.workflows.cbma import CBMAWorkflow
 
 ###############################################################################
 # Load Dataset
@@ -32,25 +32,26 @@
 ###############################################################################
 # Run CBMA Workflow
 # -----------------------------------------------------------------------------
-# The CBMA workflow function runs the following steps:
+# The fit method of a CBMA workflow class runs the following steps:
 #
 # 1. Runs a meta-analysis using the specified method (default: ALE)
 # 2. Applies a corrector to the meta-analysis results (default: FWECorrector, montecarlo)
 # 3. Generates cluster tables and runs diagnostics on the corrected results (default: Jackknife)
 #
-# All in one function call!
+# All in one call!
 #
-# result = cbma_workflow(dset)
+# result = CBMAWorkflow().fit(dset)
 #
 # For this example, we use an FDR correction because the default corrector (FWE correction with
 # Monte Carlo simulation) takes a long time to run due to the high number of iterations that
 # are required
-result = cbma_workflow(dset, corrector="fdr")
+workflow = CBMAWorkflow(corrector="fdr")
+result = workflow.fit(dset)
 
 ###############################################################################
 # Plot Results
 # -----------------------------------------------------------------------------
-# The CBMA workflow function returns a :class:`~nimare.results.MetaResult` object,
+# The fit method of the CBMA workflow class returns a :class:`~nimare.results.MetaResult` object,
 # where you can access the corrected results of the meta-analysis and diagnostics tables.
 #
 # Corrected map:

nimare/diagnostics.py

@@ -12,6 +12,8 @@
 from tqdm.auto import tqdm
 
 from nimare.base import NiMAREBase
+from nimare.meta.cbma.base import PairwiseCBMAEstimator
+from nimare.meta.ibma import IBMAEstimator
 from nimare.utils import _check_ncores, get_masker, mm2vox, tqdm_joblib
 
 LGR = logging.getLogger(__name__)
@@ -109,19 +111,10 @@ def transform(self, result):
                 correspond to positive and negative tails.
                 If no clusters are found, this list will be empty.
         """
+        pairwaise_estimators = issubclass(type(result.estimator), PairwiseCBMAEstimator)
         masker = result.estimator.masker
         diag_name = self.__class__.__name__
 
-        none_contribution_table = False
-        if not hasattr(result.estimator, "dataset"):
-            LGR.warning(
-                "MetaResult was not generated by an Estimator with a `dataset` attribute. "
-                "This may be because the Estimator was a pairwise Estimator. The "
-                "Jackknife/FocusCounter method does not currently work with pairwise Estimators. "
-                "The ``contribution_table`` will be returned as ``None``."
-            )
-            none_contribution_table = True
-
         # Collect the thresholded cluster map
         if self.target_image in result.maps:
             target_img = result.get_map(self.target_image, return_type="image")
@@ -160,7 +153,7 @@ def transform(self, result):
 
         # Define bids-like names for tables and maps
         image_name = "_".join(self.target_image.split("_")[1:])
-        image_name = "_" + image_name if image_name else image_name
+        image_name = f"_{image_name}" if image_name else image_name
         clusters_table_name = f"{self.target_image}_tab-clust"
         contribution_table_name = f"{self.target_image}_diag-{diag_name}_tab-counts"
         label_map_names = (
@@ -170,7 +163,7 @@ def transform(self, result):
         )
 
         # Check number of clusters
-        if (n_clusters == 0) or none_contribution_table:
+        if n_clusters == 0:
             result.tables[clusters_table_name] = clusters_table
             result.tables[contribution_table_name] = None
             result.maps[label_map_names[0]] = None
@@ -180,7 +173,12 @@ def transform(self, result):
 
         # Use study IDs in inputs_ instead of dataset, because we don't want to try fitting the
         # estimator to a study that might have been filtered out by the estimator's criteria.
-        meta_ids = result.estimator.inputs_["id"]
+        # Use only id1 for pairwise estimators.
+        meta_ids = (
+            result.estimator.inputs_["id1"]
+            if pairwaise_estimators
+            else result.estimator.inputs_["id"]
+        )
         rows = list(meta_ids)
 
         contribution_tables = []
@@ -204,8 +202,14 @@ def transform(self, result):
 
             contribution_tables.append(contribution_table.reset_index())
 
-        # Concat PositiveTail and NegativeTail tables
-        contribution_table = pd.concat(contribution_tables, ignore_index=True, sort=False)
+        if pairwaise_estimators or len(label_maps) == 1:
+            # Only export PositiveTail table for pairwise estimators
+            contribution_table = contribution_tables[0]
+        else:
+            # Merge PositiveTail and NegativeTail tables
+            contribution_table = pd.merge(
+                contribution_tables[0], contribution_tables[1], on="id", sort=False
+            )
 
         # Save tables and maps to result
         diag_tables_dict = {
@@ -228,14 +232,18 @@ def transform(self, result):
 class Jackknife(Diagnostics):
     """Run a jackknife analysis on a meta-analysis result.
 
+    .. versionchanged:: 0.1.1
+
+        * Support for pairwise meta-analyses.
+
     .. versionchanged:: 0.0.14
 
         * New parameter: `cluster_threshold`.
         * Return clusters table.
 
     .. versionchanged:: 0.0.13
 
-        Change cluster neighborhood from faces+edges to faces, to match Nilearn.
+        * Change cluster neighborhood from faces+edges to faces, to match Nilearn.
 
     .. versionadded:: 0.0.11
 
@@ -246,11 +254,6 @@ class Jackknife(Diagnostics):
     statistic for all experiments *except* the target experiment, dividing the resulting test
     summary statistics by the summary statistics from the original meta-analysis, and finally
     averaging the resulting proportion values across all voxels in each cluster.
-
-    Warnings
-    --------
-    Pairwise meta-analyses, like ALESubtraction and MKDAChi2, are not yet supported in this
-    method.
     """
 
     def _transform(self, expid, label_map, result):
@@ -272,11 +275,11 @@ def _transform(self, expid, label_map, result):
         """
         # We need to copy the estimator because it will otherwise overwrite the original version
         # with one missing a study in its inputs.
+        pairwaise_estimators = issubclass(type(result.estimator), PairwiseCBMAEstimator)
         estimator = copy.deepcopy(result.estimator)
 
-        dset = estimator.dataset
+        all_ids = estimator.inputs_["id1"] if pairwaise_estimators else estimator.inputs_["id"]
         original_masker = estimator.masker
-        all_ids = estimator.inputs_["id"]
 
         # Mask using a labels masker, so that we can easily get the mean value for each cluster
         cluster_masker = input_data.NiftiLabelsMasker(label_map)
@@ -288,15 +291,21 @@ def _transform(self, expid, label_map, result):
             target_value_map = "est"
         elif "stat" in result.maps:
             target_value_map = "stat"
+        elif "z_desc-consistency" in result.maps:
+            target_value_map = "z_desc-consistency"
         else:
             target_value_map = "z"
 
         stat_values = result.get_map(target_value_map, return_type="array")
 
         # Fit Estimator to all studies except the target study
         other_ids = [id_ for id_ in all_ids if id_ != expid]
-        temp_dset = dset.slice(other_ids)
-        temp_result = estimator.fit(temp_dset)
+        if pairwaise_estimators:
+            temp_dset = estimator.dataset1.slice(other_ids)
+            temp_result = estimator.fit(temp_dset, estimator.dataset2)
+        else:
+            temp_dset = estimator.dataset.slice(other_ids)
+            temp_result = estimator.fit(temp_dset)
 
         # Collect the target values (e.g., ALE values) from the N-1 meta-analysis
         temp_stat_img = temp_result.get_map(target_value_map, return_type="image")
@@ -317,6 +326,10 @@ def _transform(self, expid, label_map, result):
 class FocusCounter(Diagnostics):
     """Run a focus-count analysis on a coordinate-based meta-analysis result.
 
+    .. versionchanged:: 0.1.1
+
+        * Support for pairwise meta-analyses.
+
     .. versionchanged:: 0.0.14
 
         * New parameter: `cluster_threshold`.
@@ -337,9 +350,6 @@ class FocusCounter(Diagnostics):
     Warnings
     --------
     This method only works for coordinate-based meta-analyses.
-
-    Pairwise meta-analyses, like ALESubtraction and MKDAChi2, are not yet supported in this
-    method.
     """
 
     def _transform(self, expid, label_map, result):
@@ -359,11 +369,18 @@ def _transform(self, expid, label_map, result):
         stat_prop_values : 1D :obj:`numpy.ndarray`
             1D array with the contribution of `expid` in each cluster of `label_map`.
         """
+        if issubclass(type(result.estimator), IBMAEstimator):
+            raise ValueError("This method only works for coordinate-based meta-analyses.")
+
         affine = label_map.affine
         label_arr = label_map.get_fdata()
         clust_ids = sorted(list(np.unique(label_arr)[1:]))
 
-        coordinates_df = result.estimator.inputs_["coordinates"]
+        coordinates_df = (
+            result.estimator.inputs_["coordinates1"]
+            if issubclass(type(result.estimator), PairwiseCBMAEstimator)
+            else result.estimator.inputs_["coordinates"]
+        )
         coords = coordinates_df.loc[coordinates_df["id"] == expid]
         ijk = mm2vox(coords[["x", "y", "z"]], affine)
 

nimare/reports/base.py

@@ -29,7 +29,7 @@
 import jinja2
 from pkg_resources import resource_filename as pkgrf
 
-from nimare.meta.cbma.base import CBMAEstimator
+from nimare.meta.cbma.base import CBMAEstimator, PairwiseCBMAEstimator
 from nimare.reports.figures import (
     gen_table,
     plot_clusters,
@@ -76,6 +76,7 @@
     "fdr": "False discovery rate (FDR) correction",
     "method": "Method",
     "alpha": "Alpha",
+    "prior": "Prior",
 }
 
 SUMMARY_TEMPLATE = """\
@@ -93,6 +94,7 @@
 
 ESTIMATOR_TEMPLATE = """\
 <ul class="elem-desc">
+<li>Estimator: {est_name}</li>
 <li>Kernel Transformer: {kernel_transformer}{ker_params_text}</li>
 {est_params_text}
 </ul>
@@ -133,7 +135,10 @@ def _gen_est_summary(obj, out_filename):
         est_params_text.append(f"<li>{PARAMETERS_DICT[k]}: {v}</li>")
     est_params_text = "".join(est_params_text)
 
+    est_name = obj.__class__.__name__
+
     summary_text = ESTIMATOR_TEMPLATE.format(
+        est_name=est_name,
         kernel_transformer=str(params_dict["kernel_transformer"]),
         ker_params_text=ker_params_text,
         est_params_text=est_params_text,
@@ -359,19 +364,25 @@ def __init__(
         self.fig_dir = self.out_dir / "figures"
         self.fig_dir.mkdir(parents=True, exist_ok=True)
 
+        dataset = (
+            self.results.estimator.dataset1
+            if issubclass(type(self.results.estimator), PairwiseCBMAEstimator)
+            else self.results.estimator.dataset
+        )
+
         # Generate summary text
-        _gen_summary(self.results.estimator.dataset, self.fig_dir / "preliminary_summary.html")
+        _gen_summary(dataset, self.fig_dir / "preliminary_summary.html")
 
         # Plot mask
         plot_mask(
-            self.results.estimator.dataset.masker.mask_img,
+            dataset.masker.mask_img,
             self.fig_dir / "preliminary_figure-mask.png",
         )
 
         if issubclass(type(self.results.estimator), CBMAEstimator):
             # Plot coordinates for CBMA estimators
             plot_coordinates(
-                self.results.estimator.dataset.coordinates,
+                dataset.coordinates,
                 self.fig_dir / "preliminary_figure-static.png",
                 self.fig_dir / "preliminary_figure-interactive.html",
                 self.fig_dir / "preliminary_figure-legend.png",

nimare/reports/figures.py

@@ -81,7 +81,7 @@ def _reorder_matrix(mat, row_labels, col_labels, reorder):
 
     # Order columns
     col_linkage_matrix = linkage(mat.T, method=reorder)
-    col_ordered_linkage = optimal_leaf_ordering(col_linkage_matrix, mat)
+    col_ordered_linkage = optimal_leaf_ordering(col_linkage_matrix, mat.T)
     col_index = leaves_list(col_ordered_linkage)
 
     # Make sure labels is an ndarray and copy it
@@ -124,7 +124,6 @@ def plot_static_brain(img, out_filename, threshold=1e-06):
         black_bg=False,
         draw_cross=False,
         threshold=threshold,
-        vmax=4,
         display_mode="mosaic",
     )
     fig.savefig(out_filename, dpi=300)
@@ -251,7 +250,7 @@ def plot_interactive_brain(img, out_filename, threshold=1e-06):
     _check_extention(out_filename, [".html"])
 
     template = datasets.load_mni152_template(resolution=1)
-    html_view = view_img(img, bg_img=template, black_bg=False, threshold=threshold, vmax=4)
+    html_view = view_img(img, bg_img=template, black_bg=False, threshold=threshold)
     html_view.save_as_html(out_filename)