[ENH] Surfaced-based GLM needs the function "nilearn.glm.compute_fixed_effects" to accept 1-D array contrast and variance

[Brandon-YuHu]

Is there an existing issue for this?

• [X] I have searched the existing issues

Describe your proposed enhancement in detail.

I was trying to do surface-based GLM on my data. But my data have several sessions, and there is no function available to compute the fixed effects given the input contrast is a 1-D numpy.array instead of Nifti1Images. Your need to have an function to do that with accepting input other than Nifti1Images when surface-based analysis is implemented.

Benefits to the change

Your example of surface-based GLM is only for one session data. Your functions for surface-based analysis cannot deal with data with more than one session when there is need to compute fixed effects.

Pseudocode for the new behavior, if applicable

# insert your code below

[bthirion]

Indeed. Thx for posting. Do you think that you can contribute it ? Best,

[Remi-Gau]

FYI I could adapt the mutli run GLM example with the new surface based object to compute fixed effects:

import numpy as np
import pandas as pd

from nilearn import surface
from nilearn.datasets import func
from nilearn.experimental.surface import SurfaceImage, load_fsaverage
from nilearn.glm.first_level import FirstLevelModel

data = func.fetch_fiac_first_level()
fmri_imgs = [data["func1"], data["func2"]]

fsaverage5 = load_fsaverage()

run_imgs = []
for volume_image in [data["func1"], data["func2"]]:
    texture_left = surface.vol_to_surf(
        volume_image, fsaverage5["pial"].parts["left"]
    )
    texture_right = surface.vol_to_surf(
        volume_image, fsaverage5["pial"].parts["right"]
    )
    run_imgs.append(
        SurfaceImage(
            mesh=fsaverage5["pial"],
            data={
                "left": texture_left.T,
                "right": texture_right.T,
            },
        )
    )

design_files = [data["design_matrix1"], data["design_matrix2"]]
design_matrices = [pd.DataFrame(np.load(df)["X"]) for df in design_files]

fmri_glm = FirstLevelModel(
    minimize_memory=True,
)

fmri_glm_multirun = fmri_glm.fit(run_imgs, design_matrices=design_matrices)

n_columns = design_matrices[0].shape[1]
contrasts = {
    "SStSSp_minus_DStDSp": np.array([[1, 0, 0, -1]]),
    "DStDSp_minus_SStSSp": np.array([[-1, 0, 0, 1]]),
    "DSt_minus_SSt": np.array([[-1, -1, 1, 1]]),
    "DSp_minus_SSp": np.array([[-1, 1, -1, 1]]),
    "DSt_minus_SSt_for_DSp": np.array([[0, -1, 0, 1]]),
    "DSp_minus_SSp_for_DSt": np.array([[0, 0, -1, 1]]),
    "Deactivation": np.array([[-1, -1, -1, -1, 4]]),
    "Effects_of_interest": np.eye(n_columns)[:5, :],  # An F-contrast
}

print("Computing contrasts...")
for index, (contrast_id, contrast_val) in enumerate(contrasts.items()):
    print(f"  Contrast {index + 1:02g} out of {len(contrasts)}: {contrast_id}")
    # Estimate the contasts.
    z_map = fmri_glm.compute_contrast(contrast_val, output_type="z_score")

I can try to adapt the compute_fixed_effects to allow it to accept SurfaceImage objects for now. Not sure we want support arrays as input though.

[Remi-Gau]

I can try to adapt the compute_fixed_effects to allow it to accept SurfaceImage objects for now. Not sure we want support arrays as input though.

Will require adapting SurfaceMasker.fit() to accept list of SurfaceImage: https://nilearn.github.io/dev/modules/generated/nilearn.maskers.NiftiMasker.html#nilearn.maskers.NiftiMasker.fit

[bthirion]

Indeed rather support SurfaceImages than arrays.

[bthirion]

FYI I could adapt the mutli run GLM example with the new surface based object to compute fixed effects:

import numpy as np import pandas as pd

from nilearn import surface from nilearn.datasets import func from nilearn.experimental.surface import SurfaceImage, load_fsaverage from nilearn.glm.first_level import FirstLevelModel

data = func.fetch_fiac_first_level() fmri_imgs = [data["func1"], data["func2"]]

fsaverage5 = load_fsaverage()

run_imgs = [] for volume_image in [data["func1"], data["func2"]]: texture_left = surface.vol_to_surf( volume_image, fsaverage5["pial"].parts["left"] ) texture_right = surface.vol_to_surf( volume_image, fsaverage5["pial"].parts["right"] ) run_imgs.append( SurfaceImage( mesh=fsaverage5["pial"], data={ "left": texture_left.T, "right": texture_right.T, }, ) )

design_files = [data["design_matrix1"], data["design_matrix2"]] design_matrices = [pd.DataFrame(np.load(df)["X"]) for df in design_files]

fmri_glm = FirstLevelModel( minimize_memory=True, )

fmri_glm_multirun = fmri_glm.fit(run_imgs, design_matrices=design_matrices)

n_columns = design_matrices[0].shape[1] contrasts = { "SStSSp_minus_DStDSp": np.array([[1, 0, 0, -1]]), "DStDSp_minus_SStSSp": np.array([[-1, 0, 0, 1]]), "DSt_minus_SSt": np.array([[-1, -1, 1, 1]]), "DSp_minus_SSp": np.array([[-1, 1, -1, 1]]), "DSt_minus_SSt_for_DSp": np.array([[0, -1, 0, 1]]), "DSp_minus_SSp_for_DSt": np.array([[0, 0, -1, 1]]), "Deactivation": np.array([[-1, -1, -1, -1, 4]]), "Effects_of_interest": np.eye(n_columns)[:5, :], # An F-contrast }

print("Computing contrasts...") for index, (contrast_id, contrast_val) in enumerate(contrasts.items()): print(f" Contrast {index + 1:02g} out of {len(contrasts)}: {contrast_id}")

Estimate the contasts.

z_map = fmri_glm.compute_contrast(contrast_val, output_type="z_score")

I can try to adapt the compute_fixed_effects to allow it to accept SurfaceImage objects for now. Not sure we want support arrays as input though.

This example runs correctly on main, so there is nothing to fix ? Or did I misunderstand ?

[Remi-Gau]

to reproduce the error

import numpy as np
import pandas as pd

from nilearn import surface
from nilearn.datasets import func
from nilearn.experimental.surface import SurfaceImage, load_fsaverage
from nilearn.glm.contrasts import compute_fixed_effects
from nilearn.glm.first_level import FirstLevelModel

data = func.fetch_fiac_first_level()
fmri_imgs = [data["func1"], data["func2"]]

fsaverage5 = load_fsaverage()

run_imgs = []
for volume_image in [data["func1"], data["func2"]]:
    texture_left = surface.vol_to_surf(volume_image, fsaverage5["pial"].parts["left"])
    texture_right = surface.vol_to_surf(volume_image, fsaverage5["pial"].parts["right"])
    run_imgs.append(
        SurfaceImage(
            mesh=fsaverage5["pial"],
            data={
                "left": texture_left.T,
                "right": texture_right.T,
            },
        )
    )

design_files = [data["design_matrix1"], data["design_matrix2"]]
design_matrices = [pd.DataFrame(np.load(df)["X"]) for df in design_files]

fmri_glm = FirstLevelModel(
    minimize_memory=True,
)

contrast_id = "DSt_minus_SSt"

contrast_val = [[-1, -1, 1, 1]]

fmri_glm_run_1 = fmri_glm.fit(run_imgs[0], design_matrices=design_matrices[0])
summary_statistics_run_1 = fmri_glm_run_1.compute_contrast(
    contrast_val,
    output_type="all",
)

fmri_glm_run_2 = fmri_glm.fit(run_imgs[1], design_matrices=design_matrices[1])

contrast_val = np.array([[-1, -1, 1, 1]])

summary_statistics_run_2 = fmri_glm_run_2.compute_contrast(
    contrast_val,
    output_type="all",
)

contrast_imgs = [
    summary_statistics_run_1["effect_size"],
    summary_statistics_run_2["effect_size"],
]
variance_imgs = [
    summary_statistics_run_1["effect_variance"],
    summary_statistics_run_2["effect_variance"],
]

fixed_fx_contrast, fixed_fx_variance, fixed_fx_stat = compute_fixed_effects(
    contrast_imgs,
    variance_imgs,
)

gives this

Traceback (most recent call last):
  File "/home/remi/github/nilearn/nilearn/tmp/2_runs_surf.py", line 91, in <module>
    fixed_fx_contrast, fixed_fx_variance, fixed_fx_stat = compute_fixed_effects(
                                                          ^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remi/github/nilearn/nilearn/nilearn/glm/contrasts.py", line 529, in compute_fixed_effects
    variances = masker.transform(variance_imgs)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remi/github/nilearn/nilearn/env/lib/python3.12/site-packages/sklearn/utils/_set_output.py", line 316, in wrapped
    data_to_wrap = f(self, X, *args, **kwargs)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remi/github/nilearn/nilearn/nilearn/maskers/base_masker.py", line 249, in transform
    return self.transform_single_imgs(
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remi/github/nilearn/nilearn/nilearn/maskers/nifti_masker.py", line 578, in transform_single_imgs
    data = self._cache(
           ^^^^^^^^^^^^
  File "/home/remi/github/nilearn/nilearn/env/lib/python3.12/site-packages/joblib/memory.py", line 312, in __call__
    return self.func(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remi/github/nilearn/nilearn/nilearn/maskers/nifti_masker.py", line 95, in _filter_and_mask
    temp_imgs = _utils.check_niimg(imgs)
                ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remi/github/nilearn/nilearn/nilearn/_utils/niimg_conversions.py", line 316, in check_niimg
    return ni.image.concat_imgs(
           ^^^^^^^^^^^^^^^^^^^^^
  File "/home/remi/github/nilearn/nilearn/nilearn/image/image.py", line 1548, in concat_imgs
    first_niimg = check_niimg(next(literator), ensure_ndim=ndim)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remi/github/nilearn/nilearn/nilearn/_utils/niimg_conversions.py", line 321, in check_niimg
    niimg = load_niimg(niimg, dtype=dtype)
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/remi/github/nilearn/nilearn/nilearn/_utils/niimg.py", line 129, in load_niimg
    raise TypeError(
TypeError: Data given cannot be loaded because it is not compatible with nibabel format:
<SurfaceImage (204...