Implement cluster mass-based cluster-level Monte Carlo correction for…

… CBMA algorithms (#609)

* Implement cluster mass-based approach.

* Update base.py

* Update base.py

* Use desc keyword for both size and mass.

* Update base.py

* Update examples.

* Fix.

* Run black.

* Improve documentation.

docs/outputs.rst

@@ -15,7 +15,7 @@ Here is the basic naming convention for statistical maps:
 
 .. code-block:: Text
 
-   [value]_desc-[description]_level-[cluster|voxel]_corr-[FWE|FDR]_method-[method].nii.gz
+   <value>[_desc-<label>][_level-<cluster|voxel>][_corr-<FWE|FDR>][_method-<label>].nii.gz
 
 
 First, the ``value`` represents type of data in the map (e.g., z-statistic, t-statistic).

examples/02_meta-analyses/01_plot_cbma.py

@@ -104,7 +104,7 @@
 corr = FWECorrector(method="montecarlo", n_iters=10, n_cores=1)
 cres = corr.transform(ale.results)
 plot_stat_map(
-    cres.get_map("logp_level-cluster_corr-FWE_method-montecarlo"),
+    cres.get_map("logp_desc-size_level-cluster_corr-FWE_method-montecarlo"),
     cut_coords=[0, 0, -8],
     draw_cross=False,
     cmap="RdBu_r",

examples/02_meta-analyses/05_plot_cbma_subtraction_conjunction.py

@@ -53,7 +53,7 @@
 related_corrected_results = corr.transform(related_results)
 
 fig, axes = plt.subplots(figsize=(12, 10), nrows=2)
-img = knowledge_corrected_results.get_map("z_level-cluster_corr-FWE_method-montecarlo")
+img = knowledge_corrected_results.get_map("z_desc-size_level-cluster_corr-FWE_method-montecarlo")
 plot_stat_map(
     img,
     cut_coords=4,
@@ -66,7 +66,7 @@
     figure=fig,
 )
 
-img2 = related_corrected_results.get_map("z_level-cluster_corr-FWE_method-montecarlo")
+img2 = related_corrected_results.get_map("z_desc-size_level-cluster_corr-FWE_method-montecarlo")
 plot_stat_map(
     img2,
     cut_coords=4,
@@ -85,7 +85,7 @@
 # -----------------------------------------------------------------------------
 
 jknife = Jackknife(
-    target_image="z_level-cluster_corr-FWE_method-montecarlo",
+    target_image="z_desc-size_level-cluster_corr-FWE_method-montecarlo",
     voxel_thresh=None,
 )
 knowledge_cluster_table, knowledge_cluster_img = jknife.transform(knowledge_corrected_results)

examples/02_meta-analyses/07_plot_cbma_combine_estimators_and_kernels.py

@@ -125,7 +125,7 @@
         corr = FWECorrector(method="montecarlo", n_iters=10, n_cores=1)
         cres = corr.transform(ale.results)
         plot_stat_map(
-            cres.get_map("logp_level-cluster_corr-FWE_method-montecarlo"),
+            cres.get_map("logp_desc-size_level-cluster_corr-FWE_method-montecarlo"),
             cut_coords=[0, 0, -8],
             draw_cross=False,
             cmap="RdBu_r",

nimare/correct.py

@@ -165,7 +165,7 @@ class FWECorrector(Corrector):
     To determine what methods are available for the Estimator you're using,
     check the Estimator's documentation. Estimators have special methods
     following the naming convention correct_[correction-type]_[method]
-    (e.g., :func:`~nimare.meta.ale.ALE.correct_fwe_montecarlo`).
+    (e.g., :func:`~nimare.meta.cbma.ale.ALE.correct_fwe_montecarlo`).
     """
 
     _correction_method = "fwe"

nimare/diagnostics.py

@@ -52,7 +52,7 @@ class Jackknife(NiMAREBase):
 
     def __init__(
         self,
-        target_image="z_level-cluster_corr-FWE_method-montecarlo",
+        target_image="z_desc-size_level-cluster_corr-FWE_method-montecarlo",
         voxel_thresh=None,
         n_cores=1,
     ):

nimare/meta/cbma/base.py

@@ -472,9 +472,9 @@ def _correct_fwe_montecarlo_permutation(self, params):
 
         Returns
         -------
-        (iter_max value, iter_max_cluster)
-            A 2-tuple of floats giving the maximum voxel-wise value, and maximum
-            cluster size for the permuted dataset.
+        (iter_max value, iter_max_cluster, iter_max_mass)
+            A 3-tuple of floats giving the maximum voxel-wise value, maximum cluster size,
+            and maximum cluster mass for the permuted dataset.
         """
         iter_xyz, iter_df, conn, voxel_thresh = params
 
@@ -488,21 +488,32 @@ def _correct_fwe_montecarlo_permutation(self, params):
 
         del iter_ma_maps
 
+        # Voxel-level inference
         iter_max_value = np.max(iter_ss_map)
+
+        # Cluster-level inference
         iter_ss_map = self.masker.inverse_transform(iter_ss_map).get_fdata().copy()
         iter_ss_map[iter_ss_map <= voxel_thresh] = 0
-
         labeled_matrix = ndimage.measurements.label(iter_ss_map, conn)[0]
+        clust_vals, clust_sizes = np.unique(labeled_matrix, return_counts=True)
+        assert clust_vals[0] == 0
+
+        # Cluster mass-based inference
+        iter_max_mass = 0
+        for unique_val in clust_vals[1:]:
+            ss_vals = iter_ss_map[labeled_matrix == unique_val] - voxel_thresh
+            iter_max_mass = np.maximum(iter_max_mass, np.sum(ss_vals))
 
-        del iter_ss_map
+        del iter_ss_map, labeled_matrix
 
-        _, clust_sizes = np.unique(labeled_matrix, return_counts=True)
+        # Cluster size-based inference
         clust_sizes = clust_sizes[1:]  # First cluster is zeros in matrix
         if clust_sizes.size:
             iter_max_cluster = np.max(clust_sizes)
         else:
             iter_max_cluster = 0
-        return iter_max_value, iter_max_cluster
+
+        return iter_max_value, iter_max_cluster, iter_max_mass
 
     def correct_fwe_montecarlo(
         self, result, voxel_thresh=0.001, n_iters=10000, n_cores=-1, vfwe_only=False
@@ -511,30 +522,80 @@ def correct_fwe_montecarlo(
 
         Only call this method from within a Corrector.
 
+        .. versionchanged:: 0.0.11
+
+            * Rename ``*_level-cluster`` maps to ``*_desc-size_level-cluster``.
+            * Add new ``*_desc-mass_level-cluster`` maps that use cluster mass-based inference.
+
         Parameters
         ----------
         result : :obj:`~nimare.results.MetaResult`
-            Result object from a KDA meta-analysis.
+            Result object from a CBMA meta-analysis.
         voxel_thresh : :obj:`float`, optional
             Cluster-defining p-value threshold. Default is 0.001.
         n_iters : :obj:`int`, optional
-            Number of iterations to build the vFWE and cFWE null distributions.
-            Default is 10000.
+            Number of iterations to build the voxel-level, cluster-size, and cluster-mass FWE
+            null distributions. Default is 10000.
         n_cores : :obj:`int`, optional
             Number of cores to use for parallelization.
             If <=0, defaults to using all available cores. Default is -1.
+        vfwe_only : :obj:`bool`, optional
+            If True, only calculate the voxel-level FWE-corrected maps. Voxel-level correction
+            can be performed very quickly if the Estimator's ``null_method`` was "montecarlo".
+            If this is set to True and the original null method was not montecarlo, an exception
+            will be raised.
+            Default is False.
 
         Returns
         -------
         images : :obj:`dict`
             Dictionary of 1D arrays corresponding to masked images generated by
             the correction procedure. The following arrays are generated by
-            this method: 'vthresh', 'logp_level-cluster', and 'logp_level-voxel'.
+            this method:
+
+            -   ``logp_desc-size_level-cluster``: Cluster-level FWE-corrected ``-log10(p)`` map
+                based on cluster size. This was previously simply called "logp_level-cluster".
+                This array is **not** generated if ``vfwe_only`` is ``True``.
+            -   ``logp_desc-mass_level-cluster``: Cluster-level FWE-corrected ``-log10(p)`` map
+                based on cluster mass. According to [1]_ and [2]_, cluster mass-based inference is
+                more powerful than cluster size.
+                This array is **not** generated if ``vfwe_only`` is ``True``.
+            -   ``logp_level-voxel``: Voxel-level FWE-corrected ``-log10(p)`` map.
+                Voxel-level correction is generally more conservative than cluster-level
+                correction, so it is only recommended for very large meta-analyses
+                (i.e., hundreds of studies), per [3]_.
+
+        Notes
+        -----
+        If ``vfwe_only`` is ``False``, this method adds three new keys to the
+        ``null_distributions_`` attribute:
+
+            -   ``values_level-voxel_corr-fwe_method-montecarlo``: The maximum summary statistic
+                value from each Monte Carlo iteration. An array of shape (n_iters,).
+            -   ``values_desc-size_level-cluster_corr-fwe_method-montecarlo``: The maximum cluster
+                size from each Monte Carlo iteration. An array of shape (n_iters,).
+            -   ``values_desc-mass_level-cluster_corr-fwe_method-montecarlo``: The maximum cluster
+                mass from each Monte Carlo iteration. An array of shape (n_iters,).
 
         See Also
         --------
         nimare.correct.FWECorrector : The Corrector from which to call this method.
 
+        References
+        ----------
+        .. [1] Bullmore, E. T., Suckling, J., Overmeyer, S., Rabe-Hesketh, S., Taylor, E., &
+               Brammer, M. J. (1999). Global, voxel, and cluster tests, by theory and permutation,
+               for a difference between two groups of structural MR images of the brain.
+               IEEE transactions on medical imaging, 18(1), 32-42. doi: 10.1109/42.750253
+        .. [2] Zhang, H., Nichols, T. E., & Johnson, T. D. (2009).
+               Cluster mass inference via random field theory. Neuroimage, 44(1), 51-61.
+               doi: 10.1016/j.neuroimage.2008.08.017
+        .. [3] Eickhoff, S. B., Nichols, T. E., Laird, A. R., Hoffstaedter, F., Amunts, K.,
+               Fox, P. T., ... & Eickhoff, C. R. (2016).
+               Behavior, sensitivity, and power of activation likelihood estimation characterized
+               by massive empirical simulation. Neuroimage, 137, 70-85.
+               doi: 10.1016/j.neuroimage.2016.04.072
+
         Examples
         --------
         >>> meta = MKDADensity()
@@ -576,8 +637,7 @@ def correct_fwe_montecarlo(
             iter_df = self.inputs_["coordinates"].copy()
             iter_dfs = [iter_df] * n_iters
 
-            # Find number of voxels per cluster (includes 0, which is empty space in
-            # the matrix)
+            # Define connectivity matrix for cluster labeling
             conn = np.zeros((3, 3, 3), int)
             conn[:, :, 1] = 1
             conn[:, 1, :] = 1
@@ -601,26 +661,51 @@ def correct_fwe_montecarlo(
                         )
                     )
 
-            fwe_voxel_max, fwe_clust_max = zip(*perm_results)
+            fwe_voxel_max, fwe_cluster_size_max, fwe_cluster_mass_max = zip(*perm_results)
 
             # Cluster-level FWE
+            # Extract the summary statistics in voxel-wise (3D) form, threshold, and cluster-label
             thresh_stat_values = self.masker.inverse_transform(stat_values).get_fdata()
             thresh_stat_values[thresh_stat_values <= ss_thresh] = 0
             labeled_matrix, _ = ndimage.measurements.label(thresh_stat_values, conn)
 
-            _, idx, sizes = np.unique(labeled_matrix, return_inverse=True, return_counts=True)
-            # first cluster has value 0 (i.e., all non-zero voxels in brain), so replace
-            # with 0, which gives us a p-value of 1.
-            sizes[0] = 0
-            p_vals = null_to_p(sizes, fwe_clust_max, "upper")
-            p_cfwe_map = p_vals[np.reshape(idx, labeled_matrix.shape)]
+            cluster_labels, idx, cluster_sizes = np.unique(
+                labeled_matrix,
+                return_inverse=True,
+                return_counts=True,
+            )
+            assert cluster_labels[0] == 0
+
+            # Cluster mass-based inference
+            cluster_masses = np.zeros(cluster_labels.shape)
+            for i_val in cluster_labels:
+                if i_val == 0:
+                    cluster_masses[i_val] = 0
+
+                cluster_mass = np.sum(thresh_stat_values[labeled_matrix == i_val] - ss_thresh)
+                cluster_masses[i_val] = cluster_mass
 
-            p_cfwe_values = np.squeeze(
-                self.masker.transform(nib.Nifti1Image(p_cfwe_map, self.masker.mask_img.affine))
+            p_cmfwe_vals = null_to_p(cluster_masses, fwe_cluster_mass_max, "upper")
+            p_cmfwe_map = p_cmfwe_vals[np.reshape(idx, labeled_matrix.shape)]
+
+            p_cmfwe_values = np.squeeze(
+                self.masker.transform(nib.Nifti1Image(p_cmfwe_map, self.masker.mask_img.affine))
             )
-            logp_cfwe_values = -np.log10(p_cfwe_values)
-            logp_cfwe_values[np.isinf(logp_cfwe_values)] = -np.log10(np.finfo(float).eps)
-            z_cfwe_values = p_to_z(p_cfwe_values, tail="one")
+            logp_cmfwe_values = -np.log10(p_cmfwe_values)
+            logp_cmfwe_values[np.isinf(logp_cmfwe_values)] = -np.log10(np.finfo(float).eps)
+            z_cmfwe_values = p_to_z(p_cmfwe_values, tail="one")
+
+            # Cluster size-based inference
+            cluster_sizes[0] = 0  # replace background's "cluster size" with zeros
+            p_csfwe_vals = null_to_p(cluster_sizes, fwe_cluster_size_max, "upper")
+            p_csfwe_map = p_csfwe_vals[np.reshape(idx, labeled_matrix.shape)]
+
+            p_csfwe_values = np.squeeze(
+                self.masker.transform(nib.Nifti1Image(p_csfwe_map, self.masker.mask_img.affine))
+            )
+            logp_csfwe_values = -np.log10(p_csfwe_values)
+            logp_csfwe_values[np.isinf(logp_csfwe_values)] = -np.log10(np.finfo(float).eps)
+            z_csfwe_values = p_to_z(p_csfwe_values, tail="one")
 
             # Voxel-level FWE
             LGR.info("Using null distribution for voxel-level FWE correction.")
@@ -629,27 +714,32 @@ def correct_fwe_montecarlo(
                 "values_level-voxel_corr-fwe_method-montecarlo"
             ] = fwe_voxel_max
             self.null_distributions_[
-                "values_level-cluster_corr-fwe_method-montecarlo"
-            ] = fwe_clust_max
+                "values_desc-size_level-cluster_corr-fwe_method-montecarlo"
+            ] = fwe_cluster_size_max
+            self.null_distributions_[
+                "values_desc-mass_level-cluster_corr-fwe_method-montecarlo"
+            ] = fwe_cluster_mass_max
 
         z_vfwe_values = p_to_z(p_vfwe_values, tail="one")
         logp_vfwe_values = -np.log10(p_vfwe_values)
         logp_vfwe_values[np.isinf(logp_vfwe_values)] = -np.log10(np.finfo(float).eps)
 
         if vfwe_only:
-            # Write out unthresholded value images
+            # Return unthresholded value images
             images = {
                 "logp_level-voxel": logp_vfwe_values,
                 "z_level-voxel": z_vfwe_values,
             }
 
         else:
-            # Write out unthresholded value images
+            # Return unthresholded value images
             images = {
                 "logp_level-voxel": logp_vfwe_values,
                 "z_level-voxel": z_vfwe_values,
-                "logp_level-cluster": logp_cfwe_values,
-                "z_level-cluster": z_cfwe_values,
+                "logp_desc-size_level-cluster": logp_csfwe_values,
+                "z_desc-size_level-cluster": z_csfwe_values,
+                "logp_desc-mass_level-cluster": logp_cmfwe_values,
+                "z_desc-mass_level-cluster": z_cmfwe_values,
             }
 
         return images

nimare/tests/test_meta_ale.py

@@ -109,16 +109,27 @@ def test_ALE_approximate_null_unit(testdata_cbma, tmp_path_factory):
     corr = FWECorrector(method="montecarlo", voxel_thresh=0.001, n_iters=5, n_cores=-1)
     cres = corr.transform(meta.results)
     assert isinstance(cres, nimare.results.MetaResult)
-    assert "z_level-cluster_corr-FWE_method-montecarlo" in cres.maps.keys()
+    assert "z_desc-size_level-cluster_corr-FWE_method-montecarlo" in cres.maps.keys()
+    assert "z_desc-mass_level-cluster_corr-FWE_method-montecarlo" in cres.maps.keys()
     assert "z_level-voxel_corr-FWE_method-montecarlo" in cres.maps.keys()
-    assert "logp_level-cluster_corr-FWE_method-montecarlo" in cres.maps.keys()
+    assert "logp_desc-size_level-cluster_corr-FWE_method-montecarlo" in cres.maps.keys()
+    assert "logp_desc-mass_level-cluster_corr-FWE_method-montecarlo" in cres.maps.keys()
     assert "logp_level-voxel_corr-FWE_method-montecarlo" in cres.maps.keys()
     assert isinstance(
-        cres.get_map("z_level-cluster_corr-FWE_method-montecarlo", return_type="image"),
+        cres.get_map("z_desc-size_level-cluster_corr-FWE_method-montecarlo", return_type="image"),
         nib.Nifti1Image,
     )
     assert isinstance(
-        cres.get_map("z_level-cluster_corr-FWE_method-montecarlo", return_type="array"), np.ndarray
+        cres.get_map("z_desc-size_level-cluster_corr-FWE_method-montecarlo", return_type="array"),
+        np.ndarray,
+    )
+    assert isinstance(
+        cres.get_map("z_desc-mass_level-cluster_corr-FWE_method-montecarlo", return_type="image"),
+        nib.Nifti1Image,
+    )
+    assert isinstance(
+        cres.get_map("z_desc-mass_level-cluster_corr-FWE_method-montecarlo", return_type="array"),
+        np.ndarray,
     )
 
     # Bonferroni FWE
@@ -182,16 +193,27 @@ def test_ALE_montecarlo_null_unit(testdata_cbma, tmp_path_factory):
     corr = FWECorrector(method="montecarlo", voxel_thresh=0.001, n_iters=5, n_cores=-1)
     cres = corr.transform(meta.results)
     assert isinstance(cres, nimare.results.MetaResult)
-    assert "z_level-cluster_corr-FWE_method-montecarlo" in cres.maps.keys()
+    assert "z_desc-size_level-cluster_corr-FWE_method-montecarlo" in cres.maps.keys()
+    assert "z_desc-mass_level-cluster_corr-FWE_method-montecarlo" in cres.maps.keys()
     assert "z_level-voxel_corr-FWE_method-montecarlo" in cres.maps.keys()
-    assert "logp_level-cluster_corr-FWE_method-montecarlo" in cres.maps.keys()
+    assert "logp_desc-size_level-cluster_corr-FWE_method-montecarlo" in cres.maps.keys()
+    assert "logp_desc-mass_level-cluster_corr-FWE_method-montecarlo" in cres.maps.keys()
     assert "logp_level-voxel_corr-FWE_method-montecarlo" in cres.maps.keys()
     assert isinstance(
-        cres.get_map("z_level-cluster_corr-FWE_method-montecarlo", return_type="image"),
+        cres.get_map("z_desc-size_level-cluster_corr-FWE_method-montecarlo", return_type="image"),
+        nib.Nifti1Image,
+    )
+    assert isinstance(
+        cres.get_map("z_desc-size_level-cluster_corr-FWE_method-montecarlo", return_type="array"),
+        np.ndarray,
+    )
+    assert isinstance(
+        cres.get_map("z_desc-mass_level-cluster_corr-FWE_method-montecarlo", return_type="image"),
         nib.Nifti1Image,
     )
     assert isinstance(
-        cres.get_map("z_level-cluster_corr-FWE_method-montecarlo", return_type="array"), np.ndarray
+        cres.get_map("z_desc-mass_level-cluster_corr-FWE_method-montecarlo", return_type="array"),
+        np.ndarray,
     )
 
     # Bonferroni FWE