Skip to content

Volumetric mixed effects #342

Description

@danjgale

I am trying to run a mixed effects model using Brainstat with volumetric data (NIFTIs). I was under the impression that Brainstat can do analyses on volumes (as per the companion paper and the documentation landing page). However, I'm not really sure how to go about this, as the tutorials and documentation are surface-centric.

Is the statistics module strictly surface data after all? If not, then what's the way to go about using SLM with volumetric data?

My code resembles this following basic example, with the ??? indicating where I am confused with regards to volumetric input:

from nilearn.maskers import NiftiMasker
from brainstat.stats.terms import FixedEffect, MixedEffect
from brainstat.stats.SLM import SLM

mask = # a binary brain mask NIFTI file

design = # a pandas dataframe of the design matrix with columns 'time', 'condition', and 'subject'

time = FixedEffect(design['time'])
condition = FixedEffect(design['condition'])
subject = MixedEffect(design['subject'])
model = condition + time + subject

contrast_condition = design['condition']
slm = SLM(
    model,
    contrast_condition,
    surf=???,
    mask=???,
    correction=["fdr", "rft"],
    cluster_threshold=0.01,
)

# get input voxel array required by .fit()
masker = NiftiMasker(mask_img=mask)
voxels = masker.fit_transform(imgs)

# pass in voxel array, which I think I should be doing?
slm.fit(voxels)

The last line raises an internal type error if I set surf=mask in SLM. If I take out surf and mask completely in SLM, I still get internal type errors. What am I missing in order for it to run on volumetric data, unless I am totally mistaken?

Example of the type error I get regardless:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[22], line 1
----> 1 slm.fit(voxels)

File ~/miniconda/envs/nhp-model2/lib/python3.10/site-packages/brainstat/stats/SLM.py:154, in SLM.fit(self, Y)
    152     Y_masked = Y.copy()
    153 self._linear_model(Y_masked)
--> 154 self._t_test()
    155 if self.mask is not None:
    156     self._unmask()

File ~/miniconda/envs/nhp-model2/lib/python3.10/site-packages/brainstat/stats/_t_test.py:45, in _t_test(self)
     42         sys.exit("Contrast is not estimable :-(")
     44 else:
---> 45     c = np.dot(pinvX, self.contrast)
     46     r = self.contrast - np.dot(self.X, c)
     48     if np.square(np.ravel(r, "F")).sum() / np.square(
     49         np.ravel(self.contrast, "F")
     50     ).sum() > np.spacing(1):

File <__array_function__ internals>:180, in dot(*args, **kwargs)

TypeError: can't multiply sequence by non-int of type 'float'

Activity

  1. hechlera commented on Jul 27, 2023

    @hechlera

    I am also struggeling with making brainstat work for volumetric input. I do not get the error shown above, but multiple problems keep me from getting valid results.

    copes = # first level contrast parameter maps
    
    firstlvl_data = np.array([cope.get_fdata().flatten(order="F") for cope in copes]) # flattening with Fortran order, as this is specified in the function docs
    mask_data = img_mni_gm.get_fdata().flatten(order="F").astype(int)
    
    term_conf = FixedEffect(df_behav.memory_confidence)
    term_sub = MixedEffect(df_behav.subject)
    
    mixed_model = term_conf + term_sub
    
    slm_mixed = SLM(
        mixed_model,
        df_behav.memory_confidence,
        mask=mask_data,
        # skipping surf argument
        correction=["fdr", "rft"]
    )
    
    slm_mixed.fit(firstlvl_data)
    
    • as noted above, the argument terminology is highly confusing for volumetric data, even given the function help. The use of "surf" especially is unclear

    • using rft correction will throw an error that it into works for surface data (so I guess no cluster correction for volumetric data?)
      ValueError: Random Field Theory corrections require a surface.

    • during FDR correction, a shape mismatch error is thrown even though the input data and mask have the same format
      shape of data: (42, 1082035) shape of mask: (1082035,)

    ---------------------------------------------------------------------------
    ValueError                                Traceback (most recent call last)
    /home/tumnic/ahechler/ebits_notebooks/bold/2nd_level_nilearn.ipynb Cell 43 in ()
          [4]mixed_model = term_conf + term_sub
          [6]slm_mixed = SLM(
          [7]mixed_model,
          [8]df_behav.memory_confidence,
          [9]mask=mask_data,
         [10]correction=["fdr"]
    ---> [13]slm_mixed.fit(firstlvl_data)
    
    File [~/.conda/envs/.../lib/python3.10/site-packages/brainstat/stats/SLM.py:158), in SLM.fit(self, Y)
        156     self._unmask()
        157 if self.correction is not None:
    --> 158     self.multiple_comparison_corrections(student_t_test)
    
    File [~/.conda/envs/.../lib/python3.10/site-packages/brainstat/stats/SLM.py:234), in SLM.multiple_comparison_corrections(self, student_t_test)
        231 def multiple_comparison_corrections(self, student_t_test: bool) -> None:
        232     "Performs multiple comparisons corrections. If a (one-sided) student-t
        233     test was run, then make it two-tailed if requested."
    --> 234     P1, Q1 = self._run_multiple_comparisons()
        236     if self.two_tailed and student_t_test:
        237         self.t = -self.t
    ...
         71 Q = np.ones((self.mask.shape[0]))
    ---> 72 Q[self.mask] = np.squeeze(Q_tmp[0, :])
         74 return Q
    
    ValueError: shape mismatch: value array of shape (126006,) could not be broadcast to indexing result of shape (1082035,)
    
    • Even if I leave out the correction, I get an array of t-values that includes two 0 values and every other value is NaN
      slm_mixed.t[~np.isnan(slm_mixed.t)]
      array([0., 0.])

    Brainstat looks like an amazing tool to properly move from bash/FSL to python, especially since nilearn does not offer the mixed model approach. I am of course totally open to brainstat specializing on surface data, but currently it is not clear to me to which extent this (should) work.

  2. ageoly-git commented on Nov 16, 2023

    @ageoly-git

    Hello Experts!

    Have there been any updates on using the slm_mixed on volumetric NIFTI data? I understand that both the 'rft' feature as well as all of the plotting seems to rely on having a surface that contains both vertices and faces. The foundational paper does say that volumetric analysis is possible but so far this capability is rather elusive.

    For context, I am working with cortical thickness data in the Oasis-TRT template space (from ANTs / Mindboggle).

    Best,
    Andrew

  3. lconcha commented on May 28, 2024

    @lconcha

    Like the previous posts, I also assume that BrainStat can perform volume analyses, since SurfStat was prepared for the task, and BrainStat is just awesome.

    But, also like the others, I cannot get it to work. I cannot figure out how data must be input, since slm.fit seems to expect surface-based data. I'd be very happy if someone who got it to run on volumes could pitch in.

  4. Bella0709 commented on Jan 22, 2025

    @Bella0709

    Can VBM data be analyzed in Brainstat, or is it limited to analyzing surface data? Currently, most studies focus on analyzing surface data.

  5. zihuaihuai commented on May 4, 2026

    @zihuaihuai
    Collaborator

    Thanks for the reports, @danjgale @hechlera @ageoly-git @lconcha @Bella0709 — pulling all of these together since you're hitting variations of the same set of issues.

    Short version: BrainStat's SLM class was written for surface data and its volumetric path is, in practice, only partially functional today. The companion paper describes volumetric capability (inherited from SurfStat), but the Python implementation of the multiple-comparisons machinery, the masking utilities and the plotting helpers all assume a 2D subjects × vertices array on a surface mesh. There is no maintained, end-to-end volumetric tutorial.

    What works today (with caveats):

    1. Pass surf=None and feed a 2D (n_subjects, n_voxels) array of voxel values into SLM.fit. Use Fortran ordering (order="F") when flattening the NIfTIs, as the docstring states.
    2. Coefficient/t/p maps that come out of the linear-model + t-test step are valid; you can inverse_transform them back into NIfTI space with NiftiMasker.
    3. FDR correction will run if your mask is consistent — see (b) below.

    Why each of you hit an error:

    a. TypeError: can't multiply sequence by non-int of type 'float' (@danjgale): the contrast argument is being passed as a pandas Series of object/categorical dtype rather than a numeric column. Cast with np.asarray(design['condition'], dtype=float) before passing it in.
    b. shape mismatch: (126006,) could not be broadcast to (1082035,) (@hechlera): the FDR step computes Q over the masked data and then writes back into an array sized by mask.shape[0]. If mask is a flat int array rather than a boolean array — and especially if it has been "applied" before being passed in — the indexing fails. Pass mask as a 1D boolean array of length n_voxels (same length as the columns of Y), and let SLM do the masking.
    c. ValueError: Random Field Theory corrections require a surface. (@hechlera): expected. RFT is implemented only for surfaces. For volumes, use correction="fdr" only.
    d. All-NaN slm.t (@hechlera): symptom of the dtype problem above; once the contrast is numeric and the mask is the right shape, the t-map populates.

    Volumetric cluster-level inference (RFT/TFCE) is not supported. SurfStat had it for volumes; BrainStat's port does not. Today, the realistic options are:

    1. Run BrainStat for the GLM and FDR map.
    2. Use nilearn.glm.second_level if you want voxel-level cluster correction at the second level.
    3. Use FSL randomise / SPM for permutation-based cluster correction.

    Mixed effects with MixedEffect(subject) does run on volumetric data once (a)–(d) above are sorted, since it's just additional variance components in the design matrix; it doesn't depend on the surface.

    I'll leave this issue open as a tracking item for a proper volumetric tutorial + a documented SLM volumetric mode — the right fix is probably a small mode="volume" switch that disables RFT and ships with an example notebook. PRs welcome.

  6. zihuaihuai commented on May 4, 2026

    @zihuaihuai
    Collaborator

    Partially resolved in #381 — every SLM(...) invocation pattern that was posted in this thread now runs end-to-end on volumetric input.

    Concretely:

    • contrast: pandas Series / categorical inputs are coerced to float in SLM.__init__ instead of producing the cryptic TypeError: can't multiply sequence by non-int of type 'float' (or all-NaN t-values) downstream.
    • mask: integer 0/1 masks (the natural shape of nib.load(...).get_fdata().astype(int)) are coerced to bool, fixing the shape (126006,) could not be broadcast to indexing result of shape (1082035,) error in _fdr.
    • The "Random Field Theory corrections require a surface" error now explicitly tells you to use correction='fdr' for volumes.

    The ambitious half of this issue — cluster-level RFT / TFCE on volumes — is still not implemented. SurfStat had it; BrainStat's port doesn't, and a faithful re-port is multi-week work. Leaving this issue open as the tracker for that and labelling enhancement / help wanted. PRs welcome.

  7. added a commit that references this issue on May 4, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    PythonOn the Python implementationenhancementNew feature or requesthelp wantedExtra attention is needed

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions