2.3.1. Second-level fMRI model: a two-sample test

Full step-by-step example of fitting a GLM to perform a second level analysis in experimental data and visualizing the results.

More specifically:

  1. A sample of n=16 visual activity fMRIs are downloaded.

  2. A two-sample t-test is applied to the brain maps in order to see the effect of the contrast difference across subjects.

The contrast is between responses to vertical versus horizontal checkerboards than are retinotopically distinct. At the individual level, these stimuli are sometimes used to map the borders of primary visual areas. At the group level, such a mapping is not possible. Yet, we may observe some significant effects in these areas.

import pandas as pd
from nilearn import plotting
from nilearn.datasets import fetch_localizer_contrasts

2.3.1.1. Fetch dataset

We download a list of left vs right button press contrasts from a localizer dataset.

n_subjects = 16
sample_vertical = fetch_localizer_contrasts(
    ["vertical checkerboard"], n_subjects, get_tmaps=True)
sample_horizontal = fetch_localizer_contrasts(
    ["horizontal checkerboard"], n_subjects, get_tmaps=True)

# What remains implicit here is that there is a one-to-one
# correspondence between the two samples: the first image of both
# samples comes from subject S1, the second from subject S2 etc.

Out:

/home/kshitij/.programs/anaconda3/envs/nistats-py36-latest/lib/python3.6/site-packages/numpy/lib/npyio.py:2278: VisibleDeprecationWarning: Reading unicode strings without specifying the encoding argument is deprecated. Set the encoding, use None for the system default.
  output = genfromtxt(fname, **kwargs)

2.3.1.2. Estimate second level model

We define the input maps and the design matrix for the second level model and fit it.

second_level_input = sample_vertical['cmaps'] + sample_horizontal['cmaps']

model the effect of conditions (sample 1 vs sample 2)

import numpy as np
condition_effect = np.hstack(([1] * n_subjects, [- 1] * n_subjects))

model the subject effect: each subject is observed in sample 1 and sample 2

subject_effect = np.vstack((np.eye(n_subjects), np.eye(n_subjects)))
subjects = ['S%02d' % i for i in range(1, n_subjects + 1)]

Assemble those in a design matrix

design_matrix = pd.DataFrame(
    np.hstack((condition_effect[:, np.newaxis], subject_effect)),
    columns=['vertical vs horizontal'] + subjects)

plot the design_matrix

from nistats.reporting import plot_design_matrix
plot_design_matrix(design_matrix)
../../_images/sphx_glr_plot_second_level_two_sample_test_001.png

formally specify the analysis model and fit it

from nistats.second_level_model import SecondLevelModel
second_level_model = SecondLevelModel().fit(
    second_level_input, design_matrix=design_matrix)

Out:

/home/kshitij/.programs/anaconda3/envs/nistats-py36-latest/lib/python3.6/site-packages/nilearn/_utils/cache_mixin.py:232: DeprecationWarning: The 'cachedir' attribute has been deprecated in version 0.12 and will be removed in version 0.14.
Use os.path.join(memory.location, 'joblib') attribute instead.
  if (memory.cachedir is None and memory_level is not None
/home/kshitij/.programs/anaconda3/envs/nistats-py36-latest/lib/python3.6/site-packages/numpy/lib/function_base.py:3250: RuntimeWarning: Invalid value encountered in median
  r = func(a, **kwargs)

Estimating the contrast is very simple. We can just provide the column name of the design matrix.

z_map = second_level_model.compute_contrast('vertical vs horizontal',
                                            output_type='z_score')

Out:

/home/kshitij/.programs/anaconda3/envs/nistats-py36-latest/lib/python3.6/site-packages/numpy/core/fromnumeric.py:83: RuntimeWarning: invalid value encountered in reduce
  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)

We threshold the second level contrast and plot it

threshold = 3.1  # correponds to  p < .001, uncorrected
display = plotting.plot_glass_brain(
    z_map, threshold=threshold, colorbar=True, plot_abs=False,
    title='vertical vs horizontal checkerboard (unc p<0.001')
../../_images/sphx_glr_plot_second_level_two_sample_test_002.png

Out:

/home/kshitij/.programs/anaconda3/envs/nistats-py36-latest/lib/python3.6/site-packages/scipy/ndimage/measurements.py:272: DeprecationWarning: In future, it will be an error for 'np.bool_' scalars to be interpreted as an index
  return _nd_image.find_objects(input, max_label)

Unsurprisingly, we see activity in the primary visual cortex, both positive and negative.

Total running time of the script: ( 0 minutes 2.730 seconds)

Gallery generated by Sphinx-Gallery