Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
63 changes: 55 additions & 8 deletions brainstat/stats/SLM.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@
# type: ignore
""" Standard Linear regression models. """
import warnings
from math import sqrt
Expand All @@ -19,17 +18,23 @@
from brainstat.datasets import fetch_parcellation, fetch_template_surface
from brainstat.datasets.base import fetch_yeo_networks_metadata
from brainstat.mesh.utils import _mask_edges, mesh_edges
from brainstat.stats._linear_model import _linear_model
from brainstat.stats._multiple_comparisons import _fdr, _random_field_theory
from brainstat.stats._t_test import _t_test
from brainstat.stats.terms import FixedEffect, MixedEffect
from brainstat.stats.utils import apply_mask, undo_mask


class SLM:
"""Core Class for running BrainStat linear models"""

# Import class methods
from ._linear_model import _linear_model
from ._multiple_comparisons import _fdr, _random_field_theory
from ._t_test import _t_test
# Class methods are imported at module level (above) and assigned here
# rather than via in-class `from ... import ...`, which trips
# python/mypy#10521 and previously forced a file-level `# type: ignore`.
_linear_model = _linear_model
_fdr = _fdr
_random_field_theory = _random_field_theory
_t_test = _t_test

def __init__(
self,
Expand Down Expand Up @@ -89,7 +94,7 @@ def __init__(
"""
# Input arguments.
self.model = model
self.contrast = np.array(contrast)
self.contrast = self._coerce_contrast(contrast)

if isinstance(surf, str):
self.surf_name = surf
Expand All @@ -98,7 +103,7 @@ def __init__(
self.surf_name = None
self.surf = surf

self.mask = mask
self.mask = self._coerce_mask(mask)
self.correction = [correction] if isinstance(correction, str) else correction
self.niter = 1
self.thetalim = thetalim
Expand All @@ -110,12 +115,54 @@ def __init__(
# Error check
if self.surf is None:
if self.correction is not None and "rft" in self.correction:
raise ValueError("Random Field Theory corrections require a surface.")
raise ValueError(
"Random Field Theory corrections require a surface. For "
"volumetric data, drop 'rft' from `correction` and use "
"`correction='fdr'` (cluster-level RFT/TFCE on volumes is "
"not supported by BrainStat — see issue #342)."
)

# We have to initialize fit parameters for our unit tests here.
# TODO: remove this requirement.
self._reset_fit_parameters()

@staticmethod
def _coerce_contrast(contrast: ArrayLike) -> np.ndarray:
"""Convert ``contrast`` to a numeric ndarray.

Pandas Series and categorical/object inputs are common foot-guns —
np.dot then raises a cryptic "can't multiply sequence by non-int"
error inside _t_test. Cast to float here and surface a clear
TypeError if that fails (see issue #342).
"""
if isinstance(contrast, FixedEffect):
return contrast # type: ignore[return-value]
arr = np.asarray(contrast)
if arr.dtype.kind in "biufc":
return arr
try:
return arr.astype(float)
except (ValueError, TypeError) as exc:
raise TypeError(
f"`contrast` must be numeric; got dtype {arr.dtype}. If you "
"passed a pandas Series with categorical/object dtype, cast "
"it to float first (e.g. df['col'].astype(float))."
) from exc

@staticmethod
def _coerce_mask(mask: Optional[ArrayLike]) -> Optional[np.ndarray]:
"""Convert ``mask`` to a 1-D boolean ndarray.

Volumetric workflows commonly produce 0/1 *integer* masks (e.g.
from `nib.load(...).get_fdata().astype(int)`). The downstream
boolean-indexing path then degenerates into fancy indexing, which
triggers shape-mismatch errors in `_fdr` (see issue #342). Cast to
bool here so callers don't have to remember.
"""
if mask is None:
return None
return np.asarray(mask).astype(bool)

def fit(self, Y: np.ndarray) -> None:
"""Fits the SLM model

Expand Down
1 change: 0 additions & 1 deletion brainstat/stats/_multiple_comparisons.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@
# type: ignore
"""Multiple comparison corrections."""
import copy
import math
Expand Down
50 changes: 50 additions & 0 deletions brainstat/tests/test_SLM.py
Original file line number Diff line number Diff line change
Expand Up @@ -180,3 +180,53 @@ def test_volumetric_input():

slm = SLM(model, contrast, surf=mask_image)
slm.fit(data)


def test_volumetric_contrast_coercion():
"""Pandas Series / categorical contrasts must be coerced rather than
crashing inside _t_test (issue #342)."""
model = FixedEffect(np.ones(3))
slm = SLM(model, pd.Series([1.0, 0.5, -0.5]))
assert slm.contrast.dtype == np.float64

with pytest.raises(TypeError, match="must be numeric"):
SLM(model, pd.Series(["a", "b", "c"]))


def test_volumetric_int_mask_is_coerced_to_bool():
"""Volumetric workflows commonly hand SLM a 0/1 *int* mask
(issue #342). It should be coerced rather than degenerating into
fancy indexing inside _fdr."""
model = FixedEffect(np.ones(3))
int_mask = np.array([1, 0, 1, 1, 0, 1], dtype=int)
slm = SLM(model, np.ones(3), mask=int_mask)
assert slm.mask.dtype == np.bool_
assert slm.mask.sum() == 4


def test_volumetric_rft_without_surface_explains_alternatives():
"""RFT on a volume is unsupported; the error should point users at FDR
rather than just refusing (issue #342)."""
model = FixedEffect(np.ones(3))
with pytest.raises(ValueError, match="fdr"):
SLM(model, np.ones(3), correction="rft")


def test_volumetric_fdr_with_int_mask_runs_end_to_end():
"""The full int-mask + FDR + non-trivial contrast path that the issue
reporters originally hit should now run without error (issue #342)."""
n_subjects, n_voxels = 8, 64
rng = np.random.default_rng(0)
data = rng.standard_normal((n_subjects, n_voxels))
# Half the voxels in the mask, expressed as ints (the foot-gun shape).
int_mask = np.zeros(n_voxels, dtype=int)
int_mask[: n_voxels // 2] = 1
model = FixedEffect(rng.standard_normal(n_subjects))
contrast = pd.Series(rng.standard_normal(n_subjects)) # pandas, on purpose
slm = SLM(model, contrast, mask=int_mask, correction="fdr")
slm.fit(data)
assert slm.t.shape[-1] == n_voxels
# Inside the mask the t-values must be finite; outside they're masked
# out and undo_mask fills them with NaN.
assert np.all(np.isfinite(slm.t[..., int_mask.astype(bool)]))
assert slm.Q.shape[-1] == n_voxels
Loading