[BUG] Trouble running arithmetic in second_level_contrast of non_parametric_inference

[andreifoldes]

Is there an existing issue for this?

• I have searched the existing issues

Operating system

• Linux
• Mac
• Windows

Operating system version

Rocky Linux release 8.8 (Green Obsidian)

Python version

• 3.12
• 3.11
• 3.10
• 3.9
• 3.8

nilearn version

0.10.2

Expected behavior

I would like to a run a paired t-test using nilearn.glm.second_level.non_parametric_inference. The documentation states that

Basically one can use the name of the conditions as they appear in the design matrix of the fitted model combined with operators +- and combined with numbers with operators +-*/. "

design_matrix_df
H1 H2 Subject 1 Subject 2
0 1.0 0.0 1.0 0.0
1 0.0 1.0 1.0 0.0
2 1.0 0.0 0.0 1.0
3 0.0 1.0 0.0 1.0

out_dict = non_parametric_inference(
    searchlight_results,
    design_matrix=design_matrix_df,
    second_level_contrast="H1",
    model_intercept=True,
    n_perm=500,  # 500 for the sake of time. Ideally, this should be 10,000.
    two_sided_test=False,
    smoothing_fwhm=None,
    n_jobs=-1,
    verbose=3,
    threshold=0.01,
)

out_dict = non_parametric_inference(
    searchlight_results,
    design_matrix=design_matrix_df,
    second_level_contrast="H2",
    model_intercept=True,
    n_perm=500,  # 500 for the sake of time. Ideally, this should be 10,000.
    two_sided_test=False,
    smoothing_fwhm=None,
    n_jobs=-1,
    verbose=3,
    threshold=0.01,
)

out_dict = non_parametric_inference(
    searchlight_results,
    design_matrix=design_matrix_df,
    second_level_contrast="H1-H2",
    model_intercept=True,
    n_perm=500,  # 500 for the sake of time. Ideally, this should be 10,000.
    two_sided_test=False,
    smoothing_fwhm=None,
    n_jobs=-1,
    verbose=3,
    threshold=0.01,
)

Current behavior & error messages

I was able to run second_level_contrast="H1", second_level_contrast="H2", but not second_level_contrast="H1-H2"

{
	"name": "ValueError",
	"message": "\"H1-H2\" is not a valid contrast name",
	"stack": "---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[80], line 1
----> 1 out_dict = non_parametric_inference(
      2     searchlight_results,
      3     design_matrix=design_matrix_df,
      4     second_level_contrast=\"H1-H2\",
      5     model_intercept=True,
      6     n_perm=500,  # 500 for the sake of time. Ideally, this should be 10,000.
      7     two_sided_test=False,
      8     smoothing_fwhm=None,
      9     n_jobs=-1,
     10     verbose=3,
     11     threshold=0.01,
     12 )

File ~/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py:969, in non_parametric_inference(second_level_input, confounds, design_matrix, second_level_contrast, first_level_contrast, mask, smoothing_fwhm, model_intercept, n_perm, two_sided_test, random_state, n_jobs, verbose, threshold, tfce)
    <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=962'>963</a>     sys.stderr.write(
    <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=963'>964</a>         \"\
Computation of second level model done in \"
    <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=964'>965</a>         f\"{time.time() - t0} seconds\
\"
    <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=965'>966</a>     )
    <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=967'>968</a> # Check and obtain the contrast
--> <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=968'>969</a> contrast = _get_contrast(second_level_contrast, design_matrix)
    <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=969'>970</a> # Get first-level effect_maps
    <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=970'>971</a> effect_maps = _infer_effect_maps(second_level_input, first_level_contrast)

File ~/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py:287, in _get_contrast(second_level_contrast, design_matrix)
    <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=284'>285</a>         contrast = second_level_contrast
    <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=285'>286</a>     else:
--> <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=286'>287</a>         raise ValueError(
    <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=287'>288</a>             f'\"{second_level_contrast}\" is not a valid contrast name'
    <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=288'>289</a>         )
    <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=289'>290</a> else:
    <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=290'>291</a>     # Check contrast definition
    <a href='file:///home/saptaf1/anaconda3/envs/spy310/lib/python3.10/site-packages/nilearn/glm/second_level/second_level.py?line=291'>292</a>     if second_level_contrast is None:

ValueError: \"H1-H2\" is not a valid contrast name"
}

Steps and code to reproduce bug

import numpy as np
import nibabel as nib

def create_random_nii(filename, shape=(100, 100, 100), affine=np.eye(4)):
    data = np.random.rand(*shape)
    img = nib.Nifti1Image(data, affine)
    nib.save(img, filename)
    return filename

# Create four random .nii.gz files and add them to a list
nii_files = [create_random_nii(f'dataset{i}.nii.gz') for i in range(1, 5)]

print(nii_files)

# Number of subjects and conditions
num_subjects = len(searchlight_results_h1)
num_conditions = 2

def generate_file_names(num_conditions, num_participants):
    files = []
    for i in range(1, num_participants + 1):
        for j in range(1, num_conditions + 1):
            files.append(f'sub-{i}_cond-{j}.nii.gz')
    return files

files = generate_file_names(num_conditions, num_subjects)

# Initialize the design matrix
design_matrix = np.zeros((len(files), num_conditions + num_subjects))

# Iterate over the files
for i, file in enumerate(files):
    # Split the file name to get the subject and condition
    parts = file.split('_')
    subject = int(parts[0].split('-')[1])
    condition = int(parts[1].split('-')[1].split('.')[0])
    # Assign a value of 1 for the corresponding condition
    design_matrix[i, condition-1] = 1
    # Assign a value of 1 for the corresponding subject
    design_matrix[i, num_conditions + subject - 1] = 1

print("Design matrix:")
print(design_matrix)

# turn design matrix into a dataframe
import pandas as pd

# Convert design matrix to dataframe
design_matrix_df = pd.DataFrame(design_matrix)
# name first column H1 and second column H2 and rest of columns subject numbers
design_matrix_df.columns = ['H1', 'H2'] + [f'Subject {i}' for i in range(1, num_subjects + 1)]

out_dict = non_parametric_inference(
    nii_files,
    design_matrix=design_matrix_df,
    second_level_contrast="H1-H2",
    model_intercept=True,
    n_perm=500,  # 500 for the sake of time. Ideally, this should be 10,000.
    two_sided_test=False,
    smoothing_fwhm=None,
    n_jobs=-1,
    verbose=3,
    threshold=0.01,
)

[bthirion]

Thx for posting. Indeed this function only supports contrasts corresponding to design matrix columns and not arbitrary combinations. I think that it should call https://github.com/nilearn/nilearn/blob/main/nilearn/glm/contrasts.py#L20
The reason why it does not do so atm maybe that the added complexity would be large on the following computations ---this has to be checked.
Before this gets fixed, if it does, the best thing you can do is to compute H1-H2 in each individual and do formally inference on it. This is equivalent to the paired test.

[Remi-Gau]

isn't that a duplicate of issue #3807 ?

[ymzayek]

Yes it's the same issue. I've looked into it but didn't seem as simple as calling the expression_to_contrast_vector function to get the contrast values because non_parametric_inference calls permuted_ols which would have to be modified to support the computation IIUC

[Remi-Gau]

I think that this resonates with the conclusion I came to after a "shallow dive".

I would prefer to keep only issue for this. Which one though?

[andreifoldes]

Feel free to close this! I should have found the previous one during search! Sorry.

[Remi-Gau]

With over 250 issues open even the maintainers may be struggling to find potential duplicates, so honestly I would not expect contributors to easily find them. 😉

[Remi-Gau]

@andreifoldes
I will close this one in favor of the previous issue
I tagged you in it so hopefully you should get notifications when there is progress on this front.

[Remi-Gau]

Adapted code to reproduce and test

import numpy as np
import nibabel as nib
from nilearn.glm.second_level import non_parametric_inference

def create_random_nii(filename, shape=(20, 20, 20), affine=np.eye(4)):
    data = np.random.rand(*shape)
    img = nib.Nifti1Image(data, affine)
    nib.save(img, filename)
    return filename

# Number of subjects and conditions
num_subjects = 10
num_conditions = 2

# Create four random .nii.gz files and add them to a list
nii_files = [create_random_nii(f'dataset{i}.nii.gz') for i in range(num_subjects * num_conditions)]

print(nii_files)

def generate_file_names(num_conditions, num_participants):
    files = []
    for i in range(1, num_participants + 1):
        for j in range(1, num_conditions + 1):
            files.append(f'sub-{i}_cond-{j}.nii.gz')
    return files

files = generate_file_names(num_conditions, num_subjects)

# Initialize the design matrix
design_matrix = np.zeros((len(files), num_conditions + num_subjects))

# Iterate over the files
for i, file in enumerate(files):
    # Split the file name to get the subject and condition
    parts = file.split('_')
    subject = int(parts[0].split('-')[1])
    condition = int(parts[1].split('-')[1].split('.')[0])
    # Assign a value of 1 for the corresponding condition
    design_matrix[i, condition-1] = 1
    # Assign a value of 1 for the corresponding subject
    design_matrix[i, num_conditions + subject - 1] = 1

print("Design matrix:")
print(design_matrix)

# turn design matrix into a dataframe
import pandas as pd

# Convert design matrix to dataframe
design_matrix_df = pd.DataFrame(design_matrix)
# name first column H1 and second column H2 and rest of columns subject numbers
design_matrix_df.columns = ['H1', 'H2'] + [f'Subject {i}' for i in range(1, num_subjects + 1)]

out_dict = non_parametric_inference(
    nii_files,
    design_matrix=design_matrix_df,
    second_level_contrast="H1-H2",
    model_intercept=True,
    n_perm=100,  # 500 for the sake of time. Ideally, this should be 10,000.
    two_sided_test=False,
    smoothing_fwhm=None,
    n_jobs=-1,
    verbose=3,
    threshold=0.01,
)