Add spatial CBMR estimator and tutorial#1080
Conversation
Reviewer's GuideIntroduce a new spatially varying CBMR (SpatialCBMR) framework, including approximate and full torch backends, inference utilities, tests, and an example tutorial, and wire it into NiMARE’s public meta-analysis API. File-Level Changes
Tips and commandsInteracting with Sourcery
Customizing Your ExperienceAccess your dashboard to:
Getting Help
|
![sourcery-ai[bot]](https://avatars.githubusercontent.com/in/48477?s=60&v=4){kind=small}
There was a problem hiding this comment.
Hey - I've found 3 issues, and left some high level feedback:
- The approximate backend builds a full Kronecker preconditioner with
np.kron(moderator_info_inv, basis_info_inv), which will become extremely large for realistic numbers of moderators/bases; consider using a linear-operator style preconditioner (e.g., solving with the Kronecker factors via reshape/MatMul) instead of materializing the full matrix. - Both
_fit_fulland_fit_approximateaccept adatasetargument that is never used; consider dropping the parameter or renaming it to_datasetto make it clear it is intentionally unused and avoid confusion.
Prompt for AI Agents
Please address the comments from this code review:
## Overall Comments
- The approximate backend builds a full Kronecker preconditioner with `np.kron(moderator_info_inv, basis_info_inv)`, which will become extremely large for realistic numbers of moderators/bases; consider using a linear-operator style preconditioner (e.g., solving with the Kronecker factors via reshape/MatMul) instead of materializing the full matrix.
- Both `_fit_full` and `_fit_approximate` accept a `dataset` argument that is never used; consider dropping the parameter or renaming it to `_dataset` to make it clear it is intentionally unused and avoid confusion.
## Individual Comments
### Comment 1
<location path="nimare/meta/utils.py" line_range="129-135" />
<code_context>
+ return result.x[:n_bases], result.x[n_bases:]
+
+
+def _compute_spatial_cbmr_preconditioner(moderators, bases, mean_moderator, mean_basis, damping):
+ """Build an approximate Kronecker preconditioner for the gradient step."""
+ moderator_info = moderators.T @ (moderators * mean_moderator)
+ basis_info = bases.T @ (bases * mean_basis)
+ moderator_info_inv = np.linalg.pinv(moderator_info + damping * np.eye(moderators.shape[1]))
+ basis_info_inv = np.linalg.pinv(basis_info + damping * np.eye(bases.shape[1]))
+ return np.kron(moderator_info_inv, basis_info_inv)
+
+
</code_context>
<issue_to_address>
**issue (performance):** Preconditioner forms a potentially huge Kronecker matrix, which may not scale for typical moderator/basis sizes.
Because the preconditioner is a dense (n_moderators * n_bases)² matrix, even moderate settings (e.g., 10–20 moderators, 50–100 bases) lead to 10^8–10^10 entries, which is likely infeasible in memory and slow to apply. Since it’s only used as `preconditioner @ gradient`, consider keeping `moderator_info_inv` and `basis_info_inv` separate and implementing an implicit Kronecker–vector product instead of materializing `np.kron(...)`.
</issue_to_address>
### Comment 2
<location path="nimare/meta/models.py" line_range="43-52" />
<code_context>
+class SpatialCBMRModel(torch.nn.Module):
</code_context>
<issue_to_address>
**suggestion (performance):** Forward path densifies group foci, which can be very memory-intensive for realistic voxel grids.
`_prepare_torch_inputs` converts `foci_by_experiment_voxel[group]` to dense tensors before passing them here. For whole-brain masks (tens of thousands of voxels) and many experiments, this structure is extremely sparse, so densifying per group can cause large memory overhead and slower kernels. If you plan to support larger datasets, consider keeping foci sparse with a sparse-aware loss (e.g., summing over nonzero entries plus a background term), or at least adding a warning/config flag to control this densification.
Suggested implementation:
```python
class SpatialCBMRModel(torch.nn.Module):
"""Torch log-Poisson model for spatially varying CBMR.
This model is used by :class:`~nimare.meta.spatial_cbmr.SpatialCBMREstimator`.
For experiment ``m`` and voxel ``v`` in group ``g``, the linear predictor is
``B(v) @ alpha_g + Z_m @ beta_g @ B(v).T``.
Parameters
----------
groups : :obj:`list` of :obj:`str`
Ordered group names.
densify_foci : :obj:`bool`, optional
If ``True``, group-wise foci are converted to dense voxel × experiment
tensors before being passed through the model, which is simpler but can be
very memory-intensive for realistic whole-brain voxel grids.
If ``False``, inputs are expected to remain sparse and the estimator/model
should use a sparse-aware loss that operates on nonzero entries plus a
background term, avoiding per-group densification.
```
To fully implement this configuration instead of only documenting it, you should:
1. Update the `SpatialCBMRModel` constructor to accept the new flag, store it, and optionally warn when densifying large inputs:
- Change the `__init__` signature to include a `densify_foci: bool = True` (or similarly named) keyword argument.
- Store it on the instance (e.g., `self.densify_foci = densify_foci`).
- Optionally, when `densify_foci` is ``True`` and the voxel × experiment tensor size exceeds a heuristic threshold, issue a `warnings.warn` about potential memory usage.
2. Wire the flag into the data-preparation pipeline:
- In the `_prepare_torch_inputs` logic that currently converts `foci_by_experiment_voxel[group]` to dense tensors, branch on this flag:
- If `densify_foci` is ``True``, keep the current behavior.
- If `densify_foci` is ``False``, keep the inputs sparse (e.g., in a sparse `torch` layout or as index/value pairs) and adjust what is passed to `SpatialCBMRModel.forward`.
3. Make `forward` sparse-aware when `densify_foci` is ``False``:
- Add a code path in `SpatialCBMRModel.forward` that operates directly on sparse inputs (e.g., summing over nonzero entries plus a background term) without calling `.to_dense()` anywhere.
- Ensure the loss computation is compatible with sparse tensors and does not densify internally.
4. Propagate the configuration from the estimator:
- In `SpatialCBMREstimator` (or whatever constructs `SpatialCBMRModel`), add a matching `densify_foci` argument and pass it into the model constructor so users can control the behavior from the public API.
</issue_to_address>
### Comment 3
<location path="nimare/meta/spatial_cbmr.py" line_range="918-927" />
<code_context>
+ @classmethod
</code_context>
<issue_to_address>
**issue:** Chi-square / Wald computations rely on `np.linalg.solve` without guarding against singular or ill-conditioned covariance blocks.
In `_chi_square_log_intensity` and `_compute_spatial_coefficient_statistics`, `np.linalg.solve` is applied voxel-wise to contrast covariance matrices that can be singular or ill-conditioned due to B-spline multicollinearity and weakly informed moderators/groups. This risks `LinAlgError` or numerically unstable statistics. Consider adding a small voxel-wise ridge term to `contrast_cov` or falling back to `np.linalg.pinv` (optionally with a warning or by masking voxels with excessive condition numbers).
</issue_to_address>Help me be more useful! Please click 👍 or 👎 on each comment and I'll use the feedback to improve your reviews.
| def _compute_spatial_cbmr_preconditioner(moderators, bases, mean_moderator, mean_basis, damping): | ||
| """Build an approximate Kronecker preconditioner for the gradient step.""" | ||
| moderator_info = moderators.T @ (moderators * mean_moderator) | ||
| basis_info = bases.T @ (bases * mean_basis) | ||
| moderator_info_inv = np.linalg.pinv(moderator_info + damping * np.eye(moderators.shape[1])) | ||
| basis_info_inv = np.linalg.pinv(basis_info + damping * np.eye(bases.shape[1])) | ||
| return np.kron(moderator_info_inv, basis_info_inv) |
There was a problem hiding this comment.
issue (performance): Preconditioner forms a potentially huge Kronecker matrix, which may not scale for typical moderator/basis sizes.
Because the preconditioner is a dense (n_moderators * n_bases)² matrix, even moderate settings (e.g., 10–20 moderators, 50–100 bases) lead to 10^8–10^10 entries, which is likely infeasible in memory and slow to apply. Since it’s only used as preconditioner @ gradient, consider keeping moderator_info_inv and basis_info_inv separate and implementing an implicit Kronecker–vector product instead of materializing np.kron(...).
| class SpatialCBMRModel(torch.nn.Module): | ||
| """Torch log-Poisson model for spatially varying CBMR. | ||
|
|
||
| This model is used by :class:`~nimare.meta.spatial_cbmr.SpatialCBMREstimator`. | ||
| For experiment ``m`` and voxel ``v`` in group ``g``, the linear predictor is | ||
| ``B(v) @ alpha_g + Z_m @ beta_g @ B(v).T``. | ||
|
|
||
| Parameters | ||
| ---------- | ||
| groups : :obj:`list` of :obj:`str` |
There was a problem hiding this comment.
suggestion (performance): Forward path densifies group foci, which can be very memory-intensive for realistic voxel grids.
_prepare_torch_inputs converts foci_by_experiment_voxel[group] to dense tensors before passing them here. For whole-brain masks (tens of thousands of voxels) and many experiments, this structure is extremely sparse, so densifying per group can cause large memory overhead and slower kernels. If you plan to support larger datasets, consider keeping foci sparse with a sparse-aware loss (e.g., summing over nonzero entries plus a background term), or at least adding a warning/config flag to control this densification.
Suggested implementation:
class SpatialCBMRModel(torch.nn.Module):
"""Torch log-Poisson model for spatially varying CBMR.
This model is used by :class:`~nimare.meta.spatial_cbmr.SpatialCBMREstimator`.
For experiment ``m`` and voxel ``v`` in group ``g``, the linear predictor is
``B(v) @ alpha_g + Z_m @ beta_g @ B(v).T``.
Parameters
----------
groups : :obj:`list` of :obj:`str`
Ordered group names.
densify_foci : :obj:`bool`, optional
If ``True``, group-wise foci are converted to dense voxel × experiment
tensors before being passed through the model, which is simpler but can be
very memory-intensive for realistic whole-brain voxel grids.
If ``False``, inputs are expected to remain sparse and the estimator/model
should use a sparse-aware loss that operates on nonzero entries plus a
background term, avoiding per-group densification.To fully implement this configuration instead of only documenting it, you should:
-
Update the
SpatialCBMRModelconstructor to accept the new flag, store it, and optionally warn when densifying large inputs:- Change the
__init__signature to include adensify_foci: bool = True(or similarly named) keyword argument. - Store it on the instance (e.g.,
self.densify_foci = densify_foci). - Optionally, when
densify_fociisTrueand the voxel × experiment tensor size exceeds a heuristic threshold, issue awarnings.warnabout potential memory usage.
- Change the
-
Wire the flag into the data-preparation pipeline:
- In the
_prepare_torch_inputslogic that currently convertsfoci_by_experiment_voxel[group]to dense tensors, branch on this flag:- If
densify_fociisTrue, keep the current behavior. - If
densify_fociisFalse, keep the inputs sparse (e.g., in a sparsetorchlayout or as index/value pairs) and adjust what is passed toSpatialCBMRModel.forward.
- If
- In the
-
Make
forwardsparse-aware whendensify_fociisFalse:- Add a code path in
SpatialCBMRModel.forwardthat operates directly on sparse inputs (e.g., summing over nonzero entries plus a background term) without calling.to_dense()anywhere. - Ensure the loss computation is compatible with sparse tensors and does not densify internally.
- Add a code path in
-
Propagate the configuration from the estimator:
- In
SpatialCBMREstimator(or whatever constructsSpatialCBMRModel), add a matchingdensify_fociargument and pass it into the model constructor so users can control the behavior from the public API.
- In
| @classmethod | ||
| def _validate_method(cls, method): | ||
| """Validate and normalize an inference standard-error method.""" | ||
| if isinstance(method, str): | ||
| if method.lower() == "fi": | ||
| return "FI" | ||
| if method.lower() == "sandwich": | ||
| return "sandwich" | ||
| raise ValueError("method must be one of {'sandwich', 'FI'}.") | ||
|
|
There was a problem hiding this comment.
issue: Chi-square / Wald computations rely on np.linalg.solve without guarding against singular or ill-conditioned covariance blocks.
In _chi_square_log_intensity and _compute_spatial_coefficient_statistics, np.linalg.solve is applied voxel-wise to contrast covariance matrices that can be singular or ill-conditioned due to B-spline multicollinearity and weakly informed moderators/groups. This risks LinAlgError or numerically unstable statistics. Consider adding a small voxel-wise ridge term to contrast_cov or falling back to np.linalg.pinv (optionally with a warning or by masking voxels with excessive condition numbers).
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #1080 +/- ##
==========================================
+ Coverage 85.50% 85.96% +0.45%
==========================================
Files 56 56
Lines 11248 11849 +601
==========================================
+ Hits 9618 10186 +568
- Misses 1630 1663 +33 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
jdkent
left a comment
jdkent
left a comment
There was a problem hiding this comment.
Should SpatialCBMR be its own class? I think this should be a variant/flag of the original CBMR class, but I'm open to discussion if there is some technical barrier making it difficult to incorporate into the original CBMR class.
(How much code is reimplemented/borrowed from the original class/how much new code is being introduced?)
jdkent
left a comment
jdkent
left a comment
There was a problem hiding this comment.
Okay, taking a deeper look, my opinion on the user consuming API hasn't changed, I still want a single CBMREstimator class. with the effect being specified as follows:
CBMREstimator(moderator_effect="voxelwise|global")
the moderator effect can be voxelwise or global.
However, my opinion on the underlying code has changed. There is a significant difference in the calculation of likelihood. There can be more separation between CBMR and SpatialCBMR in the code. I'm still unsure on how different Result and Inference classes should be, and if you're adding a whole new class with duplicate methods (but different implementations of those methods), then you should define a base class with those methods stated, and then inherit that base class for CBMR and SpatialCBMR classes.
| # We simulate a coordinate-based Studyset with the same structure used in the | ||
| # CBMR tutorial: studies have reported foci, sample sizes, diagnosis labels, | ||
| # drug-status labels, and continuous study-level moderators. Spatial CBMR can be | ||
| # computationally heavier than standard CBMR because it keeps experiment-by-voxel |
There was a problem hiding this comment.
Are the experiment-by-voxel matrices kept in sparse/dense format? I've found some big speedups/memory reduction when using Compressed Sparsed Row (CSR) matrices.
There was a problem hiding this comment.
The other CBMR example should just be added to, a whole new example is not needed.
Thanks for the detailed feedback! I agree and have updated the implementation accordingly. Now, the user-facing API still uses a single CBMREstimator interface with: CBMREstimator(moderator_effect="voxelwise|global") |
jdkent
left a comment
jdkent
left a comment
There was a problem hiding this comment.
Still working through the code, nice work!
I have some comments on structure to keep it in line with the philosophy of NiMARE.
| if scipy.sparse.issparse(foci): | ||
| residual_scale, correction_factor = cls._sandwich_correction_scale( | ||
| correction, | ||
| bread_inverse, | ||
| moderators, | ||
| bases, | ||
| mean, | ||
| ) | ||
| meat_matrix = cls._sandwich_meat_matrix_sparse_response( | ||
| moderators, | ||
| bases, | ||
| foci, | ||
| mean, | ||
| meat, | ||
| residual_scale=residual_scale, | ||
| ) | ||
| else: | ||
| y = cls._as_dense_response(foci) | ||
| residuals = np.nan_to_num(y - mean, nan=0.0, posinf=0.0, neginf=0.0) | ||
| residuals, correction_factor = cls._apply_sandwich_correction( | ||
| correction, | ||
| bread_inverse, | ||
| moderators, | ||
| bases, | ||
| mean, | ||
| residuals, | ||
| ) | ||
| meat_matrix = cls._sandwich_meat_matrix(moderators, bases, residuals, meat) |
There was a problem hiding this comment.
when will the response be dense?
| * ``moderator_effect="global"`` estimates one scalar coefficient per moderator. | ||
| This is the classic CBMR model and answers whether a study-level covariate has | ||
| an overall effect on the spatial intensity function. |
There was a problem hiding this comment.
| * ``moderator_effect="global"`` estimates one scalar coefficient per moderator. | |
| This is the classic CBMR model and answers whether a study-level covariate has | |
| an overall effect on the spatial intensity function. | |
| * ``moderator_effect="global"`` estimates one scalar coefficient per moderator. This assumes the effect of the moderator has a uniform impact across the entire brain. |
Is this accurate to say? I want to make the description more applicable to the users of the algorithm, which will be focused on analyzing the brain.
There was a problem hiding this comment.
I also want to avoid using the word "study-level covariate", since the covariate/moderator could represent something about the study, or about the specific analysis within the study.
| * ``moderator_effect="voxelwise"`` estimates a smooth map for each moderator. | ||
| This spatially varying model asks where the study-level covariate effect is | ||
| stronger or weaker over voxels. |
There was a problem hiding this comment.
| * ``moderator_effect="voxelwise"`` estimates a smooth map for each moderator. | |
| This spatially varying model asks where the study-level covariate effect is | |
| stronger or weaker over voxels. | |
| * ``moderator_effect="voxelwise"`` estimates a scalar coefficient _per voxel_ per moderator. | |
| This assumes the effect of the moderator differentially impacts voxels throughout the brain. This is likely a more accurate assumption, but requires a lot more data for estimation. |
| For a more detailed introduction to the elements of a coordinate-based meta-regression, | ||
| see the | ||
| `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>`_. |
There was a problem hiding this comment.
It looks like these links were tangential, so it looks good to remove those.
| @@ -1,427 +1,315 @@ | |||
| """ | |||
| r""" | |||
There was a problem hiding this comment.
other documents do not use the raw text marker, unless you have a specific reason you need raw text and another solution will not work, I would remove this marker.
| # The same voxelwise inference helpers can use inverse Fisher information rather | ||
| # than the sandwich estimator. |
There was a problem hiding this comment.
why would a user want to do this? Explain what they gain, what different information do they get if they use inverse fisher information versus the sandwich estimator. Does the inverse Fisher information have more assumptions about the model? Is it more computationally efficient than the sandwich estimator.
| def sv_moderator_names(self): | ||
| """Return spatially varying moderator map names.""" | ||
| return tuple(name for name in self.maps if name.startswith("svModerator_")) | ||
|
|
There was a problem hiding this comment.
instead of spatially varying, the consistent language used more throughout NiMARE is "voxelwise".
There was a problem hiding this comment.
for public methods (that will be consumed by people using the tool), I want to keep docstrings of the functions that desribe the parameters:
Parameters
----------
group_contrasts : bool, dict, list, tuple, str, or None, optional
Group homogeneity or comparison specification. Use ``False`` to skip group inference.
moderator_contrasts : bool, dict, list, tuple, str, or None, optional
Moderator effect or comparison specification. Use ``False`` to skip moderator
inference.
device : str, optional
Compute device to use for inference. Defaults to the device recorded on the fitted
estimator.
"""
for internal methods that are only called by other methods within the class (and are not meant to be seen by users, add a leading underscore.
| _required_inputs = {"coordinates": ("coordinates", None)} | ||
| _group_column = "_cbmr_group" | ||
| _valid_moderator_effects = ("global", "voxelwise") | ||
| _valid_backends = ("full", "approximate") |
There was a problem hiding this comment.
is the difference between full and approximate in the documentation/in the docstring?
| @staticmethod | ||
| def _get_spatial_cbmr_approximate_solver(): | ||
| """Return the approximate solver used by the voxelwise backend.""" | ||
| return fit_spatial_cbmr_approximate |
There was a problem hiding this comment.
This should be a propery/attribute, not a method.
Store sparse experiment-by-voxel focus matrices during CBMR preprocessing so experiment-level focus counts remain aligned with grouped experiment metadata. Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Closes # .
Changes proposed in this pull request:
Summary by Sourcery
Introduce spatially varying coordinate-based meta-regression (Spatial CBMR) to NiMARE, including core estimator, inference utilities, and example usage.
New Features:
Enhancements:
Tests: