my hunch is that given that there Nilearn has no non-sphericity correction at the second level, it may be better to compute cross-session fixed effects for each participant and then do your second level model on this.
The following script (run on the main branch of nilearn)
- creates dummy data,
- runs a 2 sample t-tests where sessions are modeled at the group level and
- runs a 1 sample t-test where cross-sessions are computed at the subject level and only summary stats are brought at the group level
- compare outputs
On this dummy data, the results equivalent because the data for all subjects and conditions are randomly generated, but it may show some difference on real data.
Relevant nilearn example: Simple example of two-runs fMRI model fitting - Nilearn
from pathlib import Path
import numpy as np
import pandas as pd
from nilearn._utils.data_gen import create_fake_bids_dataset
from nilearn.glm.contrasts import compute_fixed_effects
from nilearn.glm.first_level import first_level_from_bids
from nilearn.glm.second_level import SecondLevelModel
from nilearn.plotting import (
plot_design_matrix,
plot_event,
show,
plot_img_comparison,
plot_bland_altman,
)
# generate dummy data
base_dir = Path(__file__).parent / "tmp"
base_dir.mkdir(exist_ok=True)
n_sub = 3
n_runs = 2
first_level_contrast = "c0 - c1"
bids_dataset = create_fake_bids_dataset(
base_dir=base_dir,
n_sub=n_sub,
n_ses=2,
tasks=["main"],
n_runs=[n_runs],
n_voxels=30,
)
# model for session 1
(models_ses_1, run_imgs_ses_1, events_ses_1, _) = first_level_from_bids(
base_dir / "bids_dataset",
task_label="main",
space_label="MNI",
verbose=1,
img_filters=[("ses", "01")],
)
fig = plot_event(events_ses_1[0])
show()
for model, imgs, event in zip(models_ses_1, run_imgs_ses_1, events_ses_1):
model.fit(imgs, events=event)
# model for session 2
(models_ses_2, run_imgs_ses_2, events_ses_2, _) = first_level_from_bids(
base_dir / "bids_dataset",
task_label="main",
space_label="MNI",
verbose=1,
img_filters=[("ses", "02")],
)
for model, imgs, event in zip(models_ses_2, run_imgs_ses_2, events_ses_2):
model.fit(imgs, events=event)
# second level: 2 samples T-test
# ------------------------------
second_level_input = models_ses_1 + models_ses_2
design_matrix = np.vstack(
[
np.hstack(
[
np.ones((n_sub, 1)),
np.zeros((n_sub, 1)),
np.asarray([range(1, n_sub + 1)]).T,
]
),
np.hstack(
[
np.zeros((n_sub, 1)),
np.ones((n_sub, 1)),
np.asarray([range(1, n_sub + 1)]).T,
]
),
]
)
design_matrix = pd.DataFrame(design_matrix, columns=["ses_1", "ses_2", "subject"])
fig = plot_design_matrix(design_matrix)
show()
slm = SecondLevelModel()
slm.fit(second_level_input=second_level_input, design_matrix=design_matrix)
z_score_2s_ttest = slm.compute_contrast(
second_level_contrast="ses_2 - ses_1",
first_level_contrast=first_level_contrast,
)
# second level: 1 samples T-test
# ------------------------------
# compute fixed effects for each subject
all_constrasts = []
for flm_ses_1, flm_ses_2 in zip(models_ses_1, models_ses_2):
summary_statistics_ses_1 = flm_ses_1.compute_contrast(
first_level_contrast,
output_type="all",
)
summary_statistics_ses_2 = flm_ses_2.compute_contrast(
first_level_contrast,
output_type="all",
)
contrast_imgs = [
summary_statistics_ses_1["effect_size"],
summary_statistics_ses_2["effect_size"],
]
variance_imgs = [
summary_statistics_ses_1["effect_variance"],
summary_statistics_ses_2["effect_variance"],
]
fixed_fx_contrast, _, _, _ = compute_fixed_effects(
contrast_imgs, variance_imgs, return_z_score=True
)
all_constrasts.append(fixed_fx_contrast)
slm = SecondLevelModel()
design_matrix = np.ones((n_sub, 1))
design_matrix = pd.DataFrame(design_matrix, columns=["intercept"])
fig = plot_design_matrix(design_matrix)
show()
slm.fit(second_level_input=all_constrasts, design_matrix=design_matrix)
z_score_1s_ttest = slm.compute_contrast(
second_level_contrast="intercept",
)
# Compare results
plot_img_comparison(
z_score_1s_ttest,
z_score_2s_ttest,
ref_label="1 sample t-test",
src_label="2 samples t-test",
)
plot_bland_altman(
z_score_1s_ttest,
z_score_2s_ttest,
ref_label="1 sample t-test",
src_label="2 samples t-test",
)
show()