Nilearn second-level from list of first-levels only uses the last subject's cmaps?

Software Support
,
1

Hi neurostars,

I am familiar with conducting GLM analysis in SPM for past projects. I have been transitioning my new project pipeline to nilearn and have been setting up a first- and second-level GLM but have run into some confusion about how to use nilearn’s GLM modules correctly.

For simplicity let’s say I have a series of preprocessed BOLD timeseries for different subjects in niftis_all, and each subject’s design matrix in designmats_all.

niftis_all = [...] # list where each element is a file path to preprocessed BOLD
designmats_all = [...] # list where each element is the result of make_first_level_design_matrix()

I have been struggling to find the appropriate method to fit each subject’s first-level, extract each subject’s contrast map (cmap) from their first-level, and then fit an intercept second-level to the first-level cmaps.

After a lot of tweaking, the only way I have found that will yield second-level results that appear reasonable is:

# SET UP FIRST-LEVELS
level1_base = nilearn.glm.first_level.FirstLevelModel(...) # various first-level args
level1s_all = [level1_base] * len(niftis_all) # create list with 1 empty first-level for each subject
# empty list to contain cmaps later
cmaps_all = []
# trivial example contrast on the design matrix
example_contrast = np.array([1, 0, 0, 0])

# FIT FIRST-LEVELS
for model, imgs, designmat in zip(model, niftis_all, all_subs_designmats):
    
    # fit this subject's GLM
    # inplace modification seems to work as expected
    model.fit(
        run_imgs=imgs,
        design_matrices=designmat
    )
    # THIS STEP APPEARS TO BE CRITICAL?
    # IMMEDIATELY compute this subject's first-level contrast map
    # BEFORE fitting the next subject's GLM
    z_map = model.compute_contrast(con_level1s_all, output_type="z_score")
    cmaps.append(z_map)


# SET UP AND FIT SECOND-LEVEL
design_matrix_level2 = pd.DataFrame(
    [1] * len(glm_level1s_all),
    columns=["intercept"]
)

glm_level2 = SecondLevelModel()
# input first-level contrast maps into second-level, not GLM objs
glm_level2 = glm_level2.fit(
    con_level1s_all,
    design_matrix=design_matrix_level2
)

# save report as I was doing
report_level2 = glm_level2.generate_report(
    contrasts=["intercept"],
    # example args
    alpha=.005,
    display_mode="x"
)

report_level2.save_as_html('report_level2.html')

The other first-level → second-level methods I tried seem to return funky results for me that look something like below. The spurious second-levels tend to look like this, with massive and patchy z-values. I think this pattern suggests that the second-level has accidentally been fit on repeated copies of one subject’s cmap (the final one?), as opposed to fitting across cmaps from different subjects.

image
image1930×354 45.1 KB
Here are the methods I’ve tried:

  • FAILS LIKE IMAGE: Fitting second-level model directly on the first-level model objects, as shown in BIDS dataset first- and second-level analysis - Nilearn
  • FAILS LIKE IMAGE: Fitting first-level models in a loop, next extracting first-level cmaps in a second loop, finally inputting the resulting list of cmaps into the second-level model
  • WORKS: the current strategy, extracting first-level cmaps inside the loop where first-level models are fit, as shown in Surface-based dataset first and second level analysis of a dataset - Nilearn . While I’m doing a volumetric analysis and not a surface-based analysis, this was the guide document that suggested to me the strategy of saving out each subject’s first-level cmap inside the first-level loop immediately after fitting their model.

Have I been doing something incorrectly to get results inconsistent with the guides? I’m glad I found one method that worked but for future edification I would like to understand what might have been going wrong with my other methods.

All of the examples in “Examples > GLM”, including the ones referenced in “User guide > Second level models: Second level analysis”, demonstrate second-level analysis on cmaps read from file, not estimated from first-levels in the same session, so I wasn’t really able to use them to adjudicate my issue.

Thank you for your help!

2

Thx for reporting.

I advise you to create the intermediate contrast maps to disk and to perform the second-level analysis on those. This is equivalent to working with the GLM objects, but it will help you introspect the problems.

LMK if this fails.

HTH,

3

I went ahead and wrote the intermediate contrast maps to file in both scenarios and then read them back in for visual inspection and can now confirm:

WORKS: 1 overall loop to fit each subject’s first-level model and immediately compute their contrast map before advancing to the next loop iteration

image
image1579×1580 564 KB

FAILS: Loop 1 to fit first-level models, loop 2 over list of fit model objects to compute first-level contrast maps

image
image1580×1580 580 KB

I will make sure to use method 1 going forward!

I don’t totally understand why method 2 is failing–could it be that fitting first-levels in a loop actually is overwriting the previous subject’s first-level every time? Or is it something else funky about using a list of fit first-levels as the iterator in a for loop?

4

To investigate what’s wrong here, I need to reproduce what you did. But

level1s_all = [level1_base] * len(niftis_all) # create list with 1 empty first-level for each subject

Seems to create replicas of the same model, leading to the behavior you observe. What happens if you explicitly create models in a loop ?

HTH,

Bertrand

5

Aha!! I modified my code to do like the following:

from nilearn.glm.first_level import FirstLevelModel

niftis_all = [...]
designmats_all = [...]

level1s_all = []
cmaps_all = []
example_contrast = np.array([1, 0, 0, 0])

for imgs, designmat in zip(niftis_all, designmats_all):
    
    model = FirstLevelModel(...).fit(
        run_imgs=imgs,
        design_matrices=designmat
    )

    level1s_all.append(model)
    

for model in level1s_all:
   z_map = model.compute_contrast(example_contrast, output_type="z_score")
   cmaps_all.append(z_map)

for i in range(len(cmaps_all)):
    cmaps_all[i].to_filename(f'temp/img_{i}.nii.gz')

And reading the resulting first-level cmaps back in does look okay! Each subject’s cmap looks different again.

image
image790×790 212 KB

So it would seem the thing I was doing before that was creating unexpected behavior was indeed initializing the list of models like:

level1_base = nilearn.glm.first_level.FirstLevelModel(...) # various first-level args
level1s_all = [level1_base] * len(niftis_all) # THIS ONE IS PROBLEMATIC

and when fitting an identically spec’d first-level model to a series of subjects, the proper move is NOT to pre-create the list of models but append the model objects to a list as they’re fit? I assume it’s something about how the logic of list_repeated = [item_to_repeat] * n_repeats works in Python memory?? It’s at least making more sense to me now why that approach creates unexpected results!

6

Yes, that’s it !

Best,

Bertrand

7

A more thorough option than just saving the contrast maps is using the nilearn.glm.save_glm_to_bids - Nilearn function: if you give a contrast or list of contrast and a path it will not only save the contrast image for you but also a whole bunch of details about the model.

In your case this may help to then inspect the design matrix of your model: unless you have exactly the same design across subjects (unlikely if you include motion parameters in the model), then you can quickly see if models are all the same and this points to the way you the model is constructed.

1 Like