[REF] CBMR internal and public API refactor (#1005)

* Optimize CBMR fit, inference, and device handling

Profile and refactor the CBMR hot paths while preserving numerical
behavior within tolerance.

- speed up summary table construction and inference result copying
- cache tensorized inputs, log-intensity values, and covariance terms
- vectorize multi-contrast GLH paths and add analytic Poisson Fisher info
- trim the synthetic CBMR fixture and expand regression coverage
- remove stale CBMR/model code, fix lint issues, and support CUDA devices

* update how splines are computed

* redesign the CBMR interface

* respond to review

* add useful helpers/abstractions

* fix documentation drift

* fix line length/style

* fix the documentation for cbmr

Author: jdkent

examples/02_meta-analyses/11_plot_cbmr.py

@@ -1,12 +1,12 @@
-"""
+"""Coordinate-based meta-regression algorithms.
 
 .. _metas_cbmr:
 
 ===========================================
 Coordinate-based meta-regression algorithms
 ===========================================
 
-A tour of Coordinate-based meta-regression (CBMR) algorithms in NiMARE
+A tour of Coordinate-based meta-regression (CBMR) algorithms in NiMARE.
 
 CBMR is a generative framework to approximate smooth activation intensity function
 and investigate the effect of study-level moderators (e.g., year of pubilication,
@@ -24,6 +24,7 @@
 `online course <https://www.coursera.org/lecture/functional-mri-2/module-3-meta-analysis-Vd4zz>`_
 or a `brief overview <https://libguides.princeton.edu/neuroimaging_meta>`_.
 """
+
 import numpy as np
 import scipy
 from nilearn.plotting import plot_stat_map
@@ -75,7 +76,7 @@
 ###############################################################################
 # Estimation of group-specific spatial intensity functions
 # -----------------------------------------------------------------------------
-# CBMR can generate estimation of group-specific spatial internsity
+# CBMR can generate estimates of group-specific spatial intensity
 # functions for multiple groups simultaneously, with different group-specific
 # spatial regression coefficients.
 #
@@ -84,14 +85,14 @@
 # of moderators (shared by all groups).
 #
 # Note that study-level moderators can only have global effects instead of localized
-# effects within CBMR framework. In the scenario that there're multiple subgroups
+# effects within CBMR framework. In the scenario that there are multiple subgroups
 # within a group (e.g., indexed as subgroup-1 to subgroup-n, but one or more of them
 # don't have enough number of studies to be inferred as a separate group). Using
 # categorical encoding, CBMR can interpret the subgroups as categorical moderators
 # for each study (either 0 or 1), and estimate the global activation intensity
 # associated with each subgroup (comparing to the average).
 
-from nimare.meta.cbmr import CBMREstimator
+from nimare.meta import CBMREstimator
 
 studyset = StandardizeField(fields=["sample_sizes", "avg_age"]).transform(studyset)
 
@@ -106,7 +107,7 @@
     model=models.PoissonEstimator,
     penalty=False,
     lr=1e-1,
-    tol=1e3,   # a reasonable choice is 1e-2, 1e3 is for speed
+    tol=1e3,  # a reasonable choice is 1e-2, 1e3 is for speed
     device="cpu",  # "cuda" if you have GPU
 )
 results = cbmr.fit(dataset=studyset)
@@ -163,16 +164,10 @@
 ###############################################################################
 # Generalized Linear Hypothesis (GLH) testing for spatial homogeneity
 # -----------------------------------------------------------------------------
-# In the most basic scenario of spatial homogeneity test, contrast matrix `t_con_groups`
-# can be generated by `create_contrast` function, with group names specified.
-from nimare.meta.cbmr import CBMRInference
-
-inference = CBMRInference(device="cuda")
-inference.fit(result=results)
-t_con_groups = inference.create_contrast(
-    ["SchizophreniaYes", "SchizophreniaNo", "DepressionYes", "DepressionNo"], source="groups"
-)
-contrast_result = inference.transform(t_con_groups=t_con_groups)
+# In the most basic scenario of spatial homogeneity testing, the fitted CBMR result can run
+# inference directly. The available groups and moderators are discoverable from the result.
+print(results.describe_inference_inputs())
+contrast_result = results.test_groups()
 
 ###############################################################################
 # Now that we have done spatial homogeneity tests, we can plot the z-score maps.
@@ -231,7 +226,7 @@
 # Areas with significant p-values are highlighted (under significance level :math:`0.05`).
 
 ###############################################################################
-# Perform fasle discovery rate (FDR) correction on spatial homogeneity test
+# Perform false discovery rate (FDR) correction on spatial homogeneity test
 # -----------------------------------------------------------------------------
 # The default FDR correction method is "indep", using Benjamini-Hochberg(BH) procedure.
 from nimare.correct import FDRCorrector
@@ -295,18 +290,14 @@
 ###############################################################################
 # GLH testing for group comparisons among any two groups
 # -----------------------------------------------------------------------------
-# In the most basic scenario of group comparison test, contrast matrix `t_con_groups`
-# can be generated by `create_contrast` function, with `contrast_name` specified as
-# "group1-group2".
-t_con_groups = inference.create_contrast(
+# Pairwise group comparisons can also be expressed more directly with tuples.
+contrast_result = results.compare_groups(
     [
-        "SchizophreniaYes-SchizophreniaNo",
-        "SchizophreniaNo-DepressionNo",
-        "DepressionYes-DepressionNo",
-    ],
-    source="groups",
+        ("SchizophreniaYes", "SchizophreniaNo"),
+        ("SchizophreniaNo", "DepressionNo"),
+        ("DepressionYes", "DepressionNo"),
+    ]
 )
-contrast_result = inference.transform(t_con_groups=t_con_groups, t_con_moderators=False)
 
 ###############################################################################
 # Now that we have done group comparison tests,
@@ -358,22 +349,24 @@
 ###############################################################################
 # GLH testing with contrast matrix specified
 # -----------------------------------------------------------------------------
-# CBMR supports more flexible GLH test by specifying a contrast matrix.
-# For example, group comparison test `2xgroup_0-1xgroup_1-1xgroup_2` can be
-# represented as `t_con_group=[2, -1, -1, 0]`, as an input in `compute_contrast`
-# function. Multiple independent GLH tests can be conducted simultaneously by
-# including multiple contrast vectors/matrices in `t_con_group`.
+# CBMR supports more flexible GLH tests by specifying contrast vectors or matrices
+# directly through the result-level `infer` API. For example, the group comparison
+# `2xgroup_0-1xgroup_1-1xgroup_2` can be represented as
+# `group_contrasts=[[2, -1, -1, 0]]`. Multiple independent GLH tests can be
+# conducted simultaneously by including multiple contrast vectors or matrices in
+# `group_contrasts`.
 #
-# CBMR also allows simultaneous GLH tests (consisting of multiple contrast vectors)
-# when it's represented as one of elements in `t_con_group` (datatype: list).
+# CBMR also allows simultaneous GLH tests consisting of multiple contrast vectors,
+# represented as one element of `group_contrasts`.
 # Only if all of null hypotheses are rejected at voxel level, p-values are significant.
-# For example, `t_con_group=[[1,-1,0,0], [1,0,-1,0], [0,0,1,-1]]` is used for testing
-# the equality of spatial intensity estimation among all of four groups (finding the
-# consistent activation regions). Note that only :math:`n-1` contrast vectors are necessary
-# for testing the equality of :math:`n` groups.
-
-contrast_result = inference.transform(
-    t_con_groups=[[[1, -1, 0, 0], [1, 0, -1, 0], [0, 0, 1, -1]]], t_con_moderators=False
+# For example, `[[1, -1, 0, 0], [1, 0, -1, 0], [0, 0, 1, -1]]` tests the equality
+# of spatial intensity estimates across all four groups (finding consistent activation
+# regions). Note that only :math:`n-1` contrast vectors are necessary for testing the
+# equality of :math:`n` groups.
+
+contrast_result = results.infer(
+    group_contrasts=[[[1, -1, 0, 0], [1, 0, -1, 0], [0, 0, 1, -1]]],
+    moderator_contrasts=False,
 )
 
 ###############################################################################
@@ -395,11 +388,9 @@
 ###############################################################################
 # GLH testing for study-level moderators
 # -----------------------------------------------------------------------------
-# CBMR framework can estimate global study-level moderator effects,
-# and allows inference on the existence of m.
-contrast_name = results.estimator.moderators
-t_con_moderators = inference.create_contrast(contrast_name, source="moderators")
-contrast_result = inference.transform(t_con_moderators=t_con_moderators)
+# The CBMR framework can estimate global study-level moderator effects and allows
+# inference on whether those moderator effects differ from zero.
+contrast_result = results.test_moderators()
 print(contrast_result.tables["moderators_regression_coef"])
 print(
     "P-values of moderator effects `sample_sizes` is {}".format(
@@ -420,18 +411,18 @@
 # a chosen subtype, spatial intensity estimations of the other :math:`4` subtypes of
 # schizophrenia are moderatored globally.
 
-t_con_moderators = inference.create_contrast(
-    ["standardized_sample_sizes-standardized_avg_age"], source="moderators"
+contrast_result = results.compare_moderators(
+    [("standardized_sample_sizes", "standardized_avg_age")]
 )
-contrast_result = inference.transform(t_con_moderators=t_con_moderators)
 print(
-    "P-values of difference in two moderator effectors (`sample_size-avg_age`) is {}".format(
+    "P-values of the difference between two moderator effects (`sample_size-avg_age`) is {}".format(
         contrast_result.tables["p_standardized_sample_sizes-standardized_avg_age"]
     )
 )
 
 ###############################################################################
 # CBMR also allows flexible contrasts between study-level covariates.
-# For example, we can write `contrast_name` (an input to `create_contrast`
-# function) as `standardized_sample_sizes-standardized_avg_age` when exploring
-# if the moderator effects of `sample_sizes` and `avg_age` are equivalent.
+# For example, we can express the comparison
+# `standardized_sample_sizes-standardized_avg_age` directly through
+# `results.compare_moderators(...)` when exploring whether the moderator effects
+# of `sample_sizes` and `avg_age` are equivalent.

nimare/estimator.py

@@ -12,14 +12,16 @@
 
 
 class Estimator(NiMAREBase):
-    """Estimators take in Datasets and return MetaResults.
+    """Estimators take in collections and return fitted result objects.
 
     All Estimators must have a ``_fit`` method implemented, which applies algorithm-specific
-    methods to a Dataset and returns a dictionary of arrays to be converted into a MetaResult.
+    methods to a collection and returns a dictionary of arrays to be converted into a fitted
+    result object.
 
     Users will interact with the ``_fit`` method by calling the user-facing ``fit`` method.
     ``fit`` takes in a ``Dataset``, calls ``_collect_inputs``, then ``_preprocess_input``,
-    then ``_fit``, and finally converts the dictionary returned by ``_fit`` into a ``MetaResult``.
+    then ``_fit``, and finally converts the dictionary returned by ``_fit`` into a result
+    object via ``_make_result``.
 
     .. warning::
         Support for :class:`~nimare.dataset.Dataset` inputs is deprecated and will be removed in
@@ -51,7 +53,7 @@ def _collect_inputs(self, dataset, drop_invalid=True):
 
         Parameters
         ----------
-        dataset : :obj:`~nimare.dataset.Dataset`
+        dataset : :obj:`~nimare.nimads.Studyset` or :obj:`~nimare.dataset.Dataset`
         drop_invalid : :obj:`bool`, default=True
             Whether to automatically drop any studies in the Dataset without valid data or not.
             Default is True.
@@ -111,8 +113,8 @@ def _preprocess_input(self, dataset):
 
         Parameters
         ----------
-        dataset : :obj:`~nimare.dataset.Dataset`
-            The Dataset
+        dataset : :obj:`~nimare.nimads.Studyset` or :obj:`~nimare.dataset.Dataset`
+            The collection to preprocess.
 
         .. warning::
             Support for :class:`~nimare.dataset.Dataset` inputs is deprecated and will be removed
@@ -129,21 +131,29 @@ def _fit(self, dataset):
         """
         pass
 
+    def _make_result(self, dataset, maps=None, tables=None, description=""):
+        """Construct the fitted result object for this estimator.
+
+        Subclasses may override this to return a specialized ``MetaResult`` subclass.
+        """
+        masker = getattr(self, "masker", None) or dataset.masker
+        return MetaResult(self, mask=masker, maps=maps, tables=tables, description=description)
+
     def fit(self, dataset, drop_invalid=True):
-        """Fit Estimator to a Studyset-backed collection.
+        """Fit Estimator to a collection.
 
         Parameters
         ----------
         dataset : :obj:`~nimare.nimads.Studyset` or :obj:`~nimare.dataset.Dataset`
             Collection object to analyze.
         drop_invalid : :obj:`bool`, optional
             Whether to automatically ignore any studies without the required data or not.
-            Default is False.
+            Default is True.
 
         Returns
         -------
         :obj:`~nimare.results.MetaResult`
-            Results of Estimator fitting.
+            Result of Estimator fitting. Subclasses may return a ``MetaResult`` subclass.
 
         Attributes
         ----------
@@ -166,9 +176,4 @@ def fit(self, dataset, drop_invalid=True):
         if not self.generate_description:
             description = ""
 
-        if hasattr(self, "masker") and self.masker is not None:
-            masker = self.masker
-        else:
-            masker = dataset.masker
-
-        return MetaResult(self, mask=masker, maps=maps, tables=tables, description=description)
+        return self._make_result(dataset, maps=maps, tables=tables, description=description)

nimare/generate.py

@@ -22,6 +22,13 @@
 CONTRAST_OF_INTEREST = {"animal": "as-Animal"}
 
 
+def _studyset_from_dataset(dataset):
+    """Convert a Dataset into a Studyset."""
+    from nimare.nimads import Studyset
+
+    return Studyset.from_dataset(dataset)
+
+
 def create_coordinate_dataset(
     foci=1,
     foci_percentage="100%",
@@ -145,8 +152,6 @@ def create_coordinate_studyset(
         Generated foci in xyz (mm) coordinates.
     studyset : :class:`~nimare.nimads.Studyset`
     """
-    from nimare.nimads import Studyset
-
     ground_truth_foci, dataset = create_coordinate_dataset(
         foci=foci,
         foci_percentage=foci_percentage,
@@ -157,8 +162,7 @@ def create_coordinate_studyset(
         seed=seed,
         space=space,
     )
-    studyset = Studyset.from_dataset(dataset)
-    return ground_truth_foci, studyset
+    return ground_truth_foci, _studyset_from_dataset(dataset)
 
 
 def create_neurovault_dataset(
@@ -241,16 +245,14 @@ def create_neurovault_studyset(
     :obj:`~nimare.nimads.Studyset`
         Studyset object containing experiment information from NeuroVault.
     """
-    from nimare.nimads import Studyset
-
     dataset = create_neurovault_dataset(
         collection_ids=collection_ids,
         contrasts=contrasts,
         img_dir=img_dir,
         map_type_conversion=map_type_conversion,
         **dset_kwargs,
     )
-    return Studyset.from_dataset(dataset)
+    return _studyset_from_dataset(dataset)
 
 
 def _create_source(foci, sample_sizes, space="MNI"):

nimare/io.py

@@ -1927,7 +1927,13 @@ def convert_neurovault_to_dataset(
     dataset_dict = {}
     for coll_name, nv_coll in collection_ids.items():
         nv_url = f"https://neurovault.org/api/collections/{nv_coll}/images/?format=json"
-        images = requests.get(nv_url).json()
+        images = _get_json_with_retries(
+            nv_url,
+            error_message=(
+                f"Failed to retrieve NeuroVault collection {nv_coll}. "
+                "The NeuroVault API returned an invalid or incomplete response."
+            ),
+        )
         if "Not found" in images.get("detail", ""):
             raise ValueError(
                 f"Collection {nv_coll} not found. "
@@ -1963,7 +1969,8 @@ def convert_neurovault_to_dataset(
                 )
 
                 if not filename.exists():
-                    r = requests.get(img_dict["file"])
+                    r = requests.get(img_dict["file"], timeout=30)
+                    r.raise_for_status()
                     with open(filename, "wb") as f:
                         f.write(r.content)
 
@@ -2004,6 +2011,31 @@ def convert_neurovault_to_dataset(
     return dataset
 
 
+def _get_json_with_retries(url, *, error_message, max_attempts=3, timeout=30):
+    """Fetch JSON from a remote endpoint with limited retries for transient failures."""
+    last_error = None
+    for _ in range(max_attempts):
+        try:
+            response = requests.get(url, timeout=timeout)
+            try:
+                payload = response.json()
+            except json.JSONDecodeError as exc:
+                if response.status_code == 404:
+                    return {"detail": "Not found"}
+                last_error = exc
+                continue
+
+            if response.status_code == 404:
+                return payload if isinstance(payload, dict) else {"detail": "Not found"}
+
+            response.raise_for_status()
+            return payload
+        except requests.RequestException as exc:
+            last_error = exc
+
+    raise ValueError(f"{error_message} Last error: {last_error}") from last_error
+
+
 def _resolve_sample_size(sample_sizes):
     """Choose modal sample_size if there are multiple sample_sizes to choose from."""
     sample_size_counts = Counter(sample_sizes)