search

nilearn / nilearn

Machine learning for NeuroImaging in Python
http://nilearn.github.io
Other
1.17k stars 594 forks source link

[API] Contrast images are systematically 4D images, even if they're only one-dimensional #4413

Closed bthirion closed 3 weeks ago

bthirion commented 4 months ago

Is there an existing issue for this?

Operating system

Operating system version

Python version

nilearn version

0.10.4

Expected behavior

Regarding the way we handle effect_size maps (aka beta maps) in GLM contrasts, these are always 4 dimensional even if the contrast is one dimensional. I'm unsure whether this is a good behavior. Handling such image as a 3D image would make more sense wouldn't it ? Note: my impression is that it changed recently. Would need to heck when. Note: this is to open a discussion, I don't have a very strong opinion on the matter.

Current behavior & error messages

For instance after running https://nilearn.github.io/stable/auto_examples/00_tutorials/plot_single_subject_single_run.html

In [5]: eff_map.shape
Out[5]: (64, 64, 64, 1)

Steps and code to reproduce bug

# Paste your code here
Remi-Gau commented 3 months ago

Code to reproduce

from pathlib import Path

import numpy as np

from nilearn._utils.data_gen import write_fake_fmri_data_and_design
from nilearn.glm.first_level import FirstLevelModel

output_dir = Path() / "tmp"
output_dir.mkdir(exist_ok=True, parents=True)

# %%
# one run
shapes, rk = [(7, 8, 9, 15)], 3
mask, fmri_data, design_matrices = write_fake_fmri_data_and_design(
    shapes, rk, file_path=output_dir
)
single_run_model = FirstLevelModel(mask_img=mask).fit(
    fmri_data, design_matrices=design_matrices
)
for output_type in ["z_score", "p_value", "stat", "effect_size", "effect_variance"]:
    output = single_run_model.compute_contrast(
        contrast_def=np.array([1, 0, 0]), stat_type="t", output_type=output_type
    )
    print(f"{output_type}={output.shape}")

# %%
# Two runs
shapes, rk = ((7, 8, 7, 15), (7, 8, 7, 16)), 3
mask, fmri_data, design_matrices = write_fake_fmri_data_and_design(
    shapes, rk, file_path=output_dir
)
multi_run_model = FirstLevelModel(mask_img=mask).fit(
    fmri_data, design_matrices=design_matrices
)
for output_type in ["z_score", "p_value", "stat", "effect_size", "effect_variance"]:
    output = multi_run_model.compute_contrast(
        contrast_def=np.array([1, 0, 0]), stat_type="t", output_type=output_type
    )
    print(f"{output_type}={output.shape}")
Remi-Gau commented 3 months ago

This is a "new" behavior indeed.

on 0.10.1 we used to get:

on 0.10.1✔ 24-06-06 - 16:35:20
╰─(base) ⠠⠵ /home/remi/miniconda3/bin/python /home/remi/github/nilearn/nilearn/tmp.py
z_score=(7, 8, 9)
p_value=(7, 8, 9)
stat=(7, 8, 9)
effect_size=(7, 8, 9)
effect_variance=(7, 8, 9)
/home/remi/github/nilearn/nilearn/nilearn/glm/first_level/first_level.py:706: UserWarning: One contrast given, assuming it for all 2 runs
  warn(f'One contrast given, assuming it for all {int(n_runs)} runs')
z_score=(7, 8, 7)
p_value=(7, 8, 7)
stat=(7, 8, 7)
effect_size=(7, 8, 7)
effect_variance=(7, 8, 7)

but from 0.10.2 we then got

on 0.10.2✔ 24-06-06 - 16:37:09
╰─(base) ⠠⠵ /home/remi/miniconda3/bin/python /home/remi/github/nilearn/nilearn/tmp.py
z_score=(7, 8, 9)
p_value=(7, 8, 9)
stat=(7, 8, 9)
effect_size=(7, 8, 9, 1)
effect_variance=(7, 8, 9)
/home/remi/github/nilearn/nilearn/nilearn/glm/first_level/first_level.py:799: UserWarning: One contrast given, assuming it for all 2 runs
  warn(f"One contrast given, assuming it for all {int(n_runs)} runs")
z_score=(7, 8, 7)
p_value=(7, 8, 7)
stat=(7, 8, 7)
effect_size=(7, 8, 7, 1)
effect_variance=(7, 8, 7)
Remi-Gau commented 3 months ago

Slightly modified version to check t and F contrast:

output_dir = Path() / "tmp"
output_dir.mkdir(exist_ok=True, parents=True)

# %%
shapes, rk = [(7, 8, 9, 15)], 3
mask, fmri_data, design_matrices = write_fake_fmri_data_and_design(
    shapes, rk, file_path=output_dir
)
single_run_model = FirstLevelModel(mask_img=mask).fit(
    fmri_data, design_matrices=design_matrices
)
t_contrast = np.array([1, 0, 0])
F_contrast = np.eye(rk)[:2]
for output_type in ["z_score", "p_value", "stat", "effect_size", "effect_variance"]:
    output = single_run_model.compute_contrast(
        contrast_def=t_contrast, output_type=output_type
    )
    print(f"t contrast: {output_type}={output.shape}")

    output = single_run_model.compute_contrast(
        contrast_def=F_contrast, output_type=output_type
    )
    print(f"F contrast: {output_type}={output.shape}")

0.10.1

Note that here the F contrast output is always 3 dimensional

╰─(base) ⠠⠵ /home/remi/miniconda3/bin/python /home/remi/github/nilearn/nilearn/tmp.py
t contrast: z_score=(7, 8, 9)
F contrast: z_score=(7, 8, 9)
t contrast: p_value=(7, 8, 9)
F contrast: p_value=(7, 8, 9)
t contrast: stat=(7, 8, 9)
F contrast: stat=(7, 8, 9)
t contrast: effect_size=(7, 8, 9)
F contrast: effect_size=(7, 8, 9)
t contrast: effect_variance=(7, 8, 9)
F contrast: effect_variance=(7, 8, 9)

on main

Note that here the F contrast output is always 4 dimensional for effect_size

╰─(base) ⠠⠵ /home/remi/miniconda3/bin/python /home/remi/github/nilearn/nilearn/tmp.py
t contrast: z_score=(7, 8, 9)
F contrast: z_score=(7, 8, 9)
t contrast: p_value=(7, 8, 9)
F contrast: p_value=(7, 8, 9)
t contrast: stat=(7, 8, 9)
F contrast: stat=(7, 8, 9)
t contrast: effect_size=(7, 8, 9, 1)
F contrast: effect_size=(7, 8, 9, 2)
t contrast: effect_variance=(7, 8, 9)
F contrast: effect_variance=(7, 8, 9)
Remi-Gau commented 3 months ago
Remi-Gau commented 3 months ago

Maybe this did it: https://github.com/nilearn/nilearn/pull/3203/files#r1629761075

TrevorKMDay commented 1 month ago

It looks like this issue is still being worked on, is there a workaround the team can suggest in the meantime?

Remi-Gau commented 1 month ago

double check the version but this bug was introduced in version 0.10.2, so working with nilearn==0.10.1 should 'fix' it

TrevorKMDay commented 1 month ago

Thanks for the quick response. I rolled back to 0.10.1, and that gave me a different error in the first-level modeling:

ValueError: solve: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (m,m),(m,n)->(m,n) (size 73161 is different from 1)

So I rolled back to 0.9.2, the last version I knew to be working (on Linux, I'm running on MacOS now) and get the same dimension mismatch error. Each of these times I created the conda environment from scratch. 

python                    3.10.14         h2469fbe_0_cpython    conda-forge
nibabel                   5.2.1              pyha770c72_0    conda-forge
nilearn                   0.9.2              pyhd8ed1ab_0    conda-forge
nipype                    1.8.6           py310hbe9552e_0    conda-forge

I am not sure whether this points to a package manager issue or maybe this bug is coming from a different package?

Remi-Gau commented 1 month ago

Hum... This is surprising. I am thinking that either this something to do with the data you are tryting to analyze or some interaction with an older version of nilearn and more recent versions of its dependencies. If you can do you think you could try to run one of the nilearn examples that covers what you are tying to do? If that fails then we can rule out the data issue.

TrevorKMDay commented 1 month ago

Sure, I will try that and let you know.

bthirion commented 1 month ago

Thx for the feedback. As a hotfix, can't you coerce the contast images to 3D manually (eg by simple reshaping) ?

TrevorKMDay commented 1 month ago

Unless I'm misunderstanding, it's buried in the methods

    sec_lvl_model = sec_lvl_model.fit(second_level_input=fixed_files,
                                      design_matrix=design_matrix)

    # Calculate z-statistic from second lvl map
    zstat_map = sec_lvl_model.compute_contrast(
        second_level_contrast='Intercept',
        second_level_stat_type='t',
        output_type='z_score',
    )

compute_contrast() is what throws the error, and fixed_files is just a list of strings pointing to the maps on disk, so there is no chance for my to intervene and reshape.

I could load and reshape, and then invoke compute_contrast(), I suppose?

bthirion commented 1 month ago

I could load and reshape, and then invoke compute_contrast(), I suppose?

I think so. Sorry for the mess

TrevorKMDay commented 1 month ago

Well, this is interesting. I pass a list of files using shell globbing, so the list is explicit. Under 0.10.1, it initially gives me a doubled list:

['foo/sub-1_stat-effect.nii.gz', 'foo//sub-1_stat-effect.nii.gz']

I de-duped those and looped over the list:

for i in fixed_files2:
    x = nil.image.load_img(i)
    print(x.shape)

This works fine and gives me (54, 64, 55) for each image. But as soon as I try to nil.image.load_img() more than one file at a time, it gives me the error below again:

Reference shape: 
(54, 66, 55)
Input shape:
(54, 65, 55, 1)

I can use a list comprehension to get the list of loaded images, but then that fails again in compute_contrast().

bthirion commented 1 month ago

@Remi-Gau Are you on this ? LMK if you need me to get into it.

Remi-Gau commented 1 month ago

I wanted to work on this in #4439 but I quickly got the feeling that you'd be faster and troubleshooting this one. I won't have time to work on this in the coming week though.

In any case we want to have this fixed in our next release.

TrevorKMDay commented 1 month ago

If it helps, I noticed a working file had pixdim values of [-1 3 3 3 0 0 0 0] but the failing file had [1 3 3 3 2.996 1 1 1 1] so I edited the header, thinking it was checking one of those values, but it didn't help.

(The 2.996 pixdim3 value is because fMRIPREP left it in native, but SPM conformed it to 3 mm isotropic).