Note

Go to the end to download the full example code. or to run this example in your browser via Binder

Voxel-Based Morphometry on OASIS dataset

This example uses voxel-based morphometry (VBM) to study the relationship between aging, sex, and gray matter density.

The data come from the OASIS project. If you use it, you need to agree with the data usage agreement available on the website.

It has been run through a standard VBM pipeline (using SPM8 and NewSegment) to create VBM maps, which we study here.

VBM analysis of aging

We run a standard GLM analysis to study the association between age and gray matter density from the VBM data. We use only 100 subjects from the OASIS dataset to limit the memory usage.

Note that more power would be obtained from using a larger sample of subjects.

See also

For more information see the dataset description.

Load Oasis dataset

from nilearn import plotting
from nilearn.datasets import (
    fetch_icbm152_2009,
    fetch_icbm152_brain_gm_mask,
    fetch_oasis_vbm,
)

n_subjects = 100  # more subjects requires more memory

oasis_dataset = fetch_oasis_vbm(
    n_subjects=n_subjects,
)
gray_matter_map_filenames = oasis_dataset.gray_matter_maps
age = oasis_dataset.ext_vars["age"].astype(float)

[fetch_oasis_vbm] Dataset found in /home/runner/nilearn_data/oasis1

Sex is encoded as ‘M’ or ‘F’. Hence, we make it a binary variable.

sex = oasis_dataset.ext_vars["mf"] == "F"

Print basic information on the dataset.

print(
    "First gray-matter anatomy image (3D) is located at: "
    f"{oasis_dataset.gray_matter_maps[0]}"
)
print(
    "First white-matter anatomy image (3D) is located at: "
    f"{oasis_dataset.white_matter_maps[0]}"
)

First gray-matter anatomy image (3D) is located at: /home/runner/nilearn_data/oasis1/OAS1_0001_MR1/mwrc1OAS1_0001_MR1_mpr_anon_fslswapdim_bet.nii.gz
First white-matter anatomy image (3D) is located at: /home/runner/nilearn_data/oasis1/OAS1_0001_MR1/mwrc2OAS1_0001_MR1_mpr_anon_fslswapdim_bet.nii.gz

Get a mask image: A mask of the cortex of the ICBM template.

gm_mask = fetch_icbm152_brain_gm_mask()

[fetch_icbm152_brain_gm_mask] Dataset found in
/home/runner/nilearn_data/icbm152_2009

Resample the mask, since this mask has a different resolution.

from nilearn.image import resample_to_img

mask_img = resample_to_img(
    gm_mask,
    gray_matter_map_filenames[0],
    interpolation="nearest",
    copy_header=True,
    force_resample=True,
)

Analyze data

First, we create an adequate design matrix with three columns: ‘age’, ‘sex’, and ‘intercept’.

import numpy as np
import pandas as pd

intercept = np.ones(n_subjects)
design_matrix = pd.DataFrame(
    np.vstack((age, sex, intercept)).T,
    columns=["age", "sex", "intercept"],
)

from matplotlib import pyplot as plt

Let’s plot the design matrix.

fig, ax1 = plt.subplots(1, 1, figsize=(4, 8))
ax = plotting.plot_design_matrix(design_matrix, axes=ax1)
ax.set_ylabel("maps")
fig.suptitle("Second level design matrix")
Second level design matrix
Text(0.5, 0.98, 'Second level design matrix')

Next, we specify and fit the second-level model when loading the data and also smooth a little bit to improve statistical behavior.

from nilearn.glm.second_level import SecondLevelModel

second_level_model = SecondLevelModel(
    smoothing_fwhm=2.0,
    mask_img=mask_img,
    n_jobs=2,
    minimize_memory=False,
    verbose=1,
)
second_level_model.fit(
    gray_matter_map_filenames,
    design_matrix=design_matrix,
)

[SecondLevelModel.fit] Fitting second level model. Take a deep breath.
[SecondLevelModel.fit] Loading data from Nifti1Image(
shape=(91, 109, 91),
affine=array([[  -2.,    0.,    0.,   90.],
       [   0.,    2.,    0., -126.],
       [   0.,    0.,    2.,  -72.],
       [   0.,    0.,    0.,    1.]])
)
[SecondLevelModel.fit] loading mask from Nifti1Image(
shape=(91, 109, 91),
affine=array([[  -2.,    0.,    0.,   90.],
       [   0.,    2.,    0., -126.],
       [   0.,    0.,    2.,  -72.],
       [   0.,    0.,    0.,    1.]])
)
/home/runner/work/nilearn/nilearn/examples/05_glm_second_level/plot_oasis.py:116: UserWarning: [NiftiMasker.fit] Generation of a mask has been requested (imgs != None) while a mask was given at masker creation. Given mask will be used.
  second_level_model.fit(
[SecondLevelModel.fit] Resampling mask
[SecondLevelModel.fit] Finished fit
[SecondLevelModel.fit]
Computation of second level model done in 0.28 seconds.

Second Level Model
In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.


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

z_map = second_level_model.compute_contrast(
    second_level_contrast=[1, 0, 0],
    output_type="z_score",
)

[SecondLevelModel.compute_contrast] Loading data from Nifti1Image(
shape=(91, 109, 91, 100),
affine=array([[  -2.,    0.,    0.,   90.],
       [   0.,    2.,    0., -126.],
       [   0.,    0.,    2.,  -72.],
       [   0.,    0.,    0.,    1.]])
)
[SecondLevelModel.compute_contrast] Smoothing images
[SecondLevelModel.compute_contrast] Extracting region signals
[SecondLevelModel.compute_contrast] Cleaning extracted signals

We threshold the second level contrast at FDR-corrected p < 0.05 and plot it.

from nilearn.glm import threshold_stats_img

_, threshold = threshold_stats_img(z_map, alpha=0.05, height_control="fdr")
print(f"The FDR=.05-corrected threshold is: {threshold:03g}")

fig = plt.figure(figsize=(5, 3))
display = plotting.plot_stat_map(
    z_map,
    threshold=threshold,
    display_mode="z",
    cut_coords=[-4, 26],
    figure=fig,
)
fig.suptitle("age effect on gray matter density (FDR = .05)")
plotting.show()
age effect on gray matter density (FDR = .05)
The FDR=.05-corrected threshold is: 2.40175

We can also study the effect of sex by computing the contrast, thresholding it and plot the resulting map.

z_map = second_level_model.compute_contrast(
    second_level_contrast="sex",
    output_type="z_score",
)
_, threshold = threshold_stats_img(z_map, alpha=0.05, height_control="fdr")
plotting.plot_stat_map(
    z_map,
    threshold=threshold,
    title="sex effect on gray matter density (FDR = .05)",
)
plot oasis
[SecondLevelModel.compute_contrast] Loading data from Nifti1Image(
shape=(91, 109, 91, 100),
affine=array([[  -2.,    0.,    0.,   90.],
       [   0.,    2.,    0., -126.],
       [   0.,    0.,    2.,  -72.],
       [   0.,    0.,    0.,    1.]])
)
[SecondLevelModel.compute_contrast] Smoothing images
[SecondLevelModel.compute_contrast] Extracting region signals
[SecondLevelModel.compute_contrast] Cleaning extracted signals

<nilearn.plotting.displays._slicers.OrthoSlicer object at 0x7f8eda1461f0>

Note that there does not seem to be any significant effect of sex on gray matter density on that dataset.

Saving model outputs to disk

It can be useful to quickly generate a portable, ready-to-view report with most of the pertinent information. We can do this by saving the output of the GLM to disk including an HTML report. This is easy to do if you have a fitted model and the list of contrasts, which we do here.

from pathlib import Path

from nilearn.interfaces.bids import save_glm_to_bids

output_dir = Path.cwd() / "results" / "plot_oasis"
output_dir.mkdir(exist_ok=True, parents=True)

icbm152_2009 = fetch_icbm152_2009()

second_level_model = save_glm_to_bids(
    second_level_model,
    contrasts=["age", "sex"],
    out_dir=output_dir / "derivatives" / "nilearn_glm",
    prefix="ageEffectOnGM",
    bg_img=icbm152_2009["t1"],
    alpha=0.05,
    height_control="fdr",
)

[fetch_icbm152_2009] Dataset found in /home/runner/nilearn_data/icbm152_2009
[save_glm_to_bids] Generating design matrices figures...
[save_glm_to_bids] Generating contrast matrices figures...
/home/runner/work/nilearn/nilearn/.tox/doc/lib/python3.9/site-packages/nilearn/_utils/plotting.py:154: UserWarning: constrained_layout not applied.  At least one axes collapsed to zero width or height.
  contrast_plot.figure.savefig(output["dir"] / contrast_fig)
[save_glm_to_bids] Saving contrast-level statistical maps...
[SecondLevelModel.compute_contrast] Loading data from Nifti1Image(
shape=(91, 109, 91, 100),
affine=array([[  -2.,    0.,    0.,   90.],
       [   0.,    2.,    0., -126.],
       [   0.,    0.,    2.,  -72.],
       [   0.,    0.,    0.,    1.]])
)
[SecondLevelModel.compute_contrast] Smoothing images
[SecondLevelModel.compute_contrast] Extracting region signals
[SecondLevelModel.compute_contrast] Cleaning extracted signals
[SecondLevelModel.compute_contrast] Loading data from Nifti1Image(
shape=(91, 109, 91, 100),
affine=array([[  -2.,    0.,    0.,   90.],
       [   0.,    2.,    0., -126.],
       [   0.,    0.,    2.,  -72.],
       [   0.,    0.,    0.,    1.]])
)
[SecondLevelModel.compute_contrast] Smoothing images
[SecondLevelModel.compute_contrast] Extracting region signals
[SecondLevelModel.compute_contrast] Cleaning extracted signals
/home/runner/work/nilearn/nilearn/examples/05_glm_second_level/plot_oasis.py:184: UserWarning: The given float value must not exceed 2.75855356998246e-06. But, you have given threshold=inf.
  second_level_model = save_glm_to_bids(
/home/runner/work/nilearn/nilearn/examples/05_glm_second_level/plot_oasis.py:184: UserWarning: The given float value must not exceed 0.0. But, you have given threshold=inf.
  second_level_model = save_glm_to_bids(
/home/runner/work/nilearn/nilearn/examples/05_glm_second_level/plot_oasis.py:184: UserWarning: Attention: No clusters with stat higher than inf
  second_level_model = save_glm_to_bids(
/home/runner/work/nilearn/nilearn/examples/05_glm_second_level/plot_oasis.py:184: UserWarning: The given float value must not exceed 0.00743128949229107. But, you have given threshold=inf.
  second_level_model = save_glm_to_bids(
/home/runner/work/nilearn/nilearn/examples/05_glm_second_level/plot_oasis.py:184: UserWarning: The given float value must not exceed 0.0. But, you have given threshold=inf.
  second_level_model = save_glm_to_bids(
/home/runner/work/nilearn/nilearn/examples/05_glm_second_level/plot_oasis.py:184: UserWarning: Attention: No clusters with stat higher than inf
  second_level_model = save_glm_to_bids(
[save_glm_to_bids] Saving model level statistical maps...
[save_glm_to_bids] Generating HTML...

View the generated files

files = sorted((output_dir / "derivatives" / "nilearn_glm").glob("**/*"))
print("\n".join([str(x.relative_to(output_dir)) for x in files]))

#  %%
# Generate a report and view it.
# If no new contrast is passed to ``generate_report``,
# the results saved to disk will be reused to generate the report.
report = second_level_model.generate_report(
    bg_img=icbm152_2009["t1"],
    plot_type="glass",
    alpha=0.05,
    height_control=None,
)
report

derivatives/nilearn_glm/dataset_description.json
derivatives/nilearn_glm/group
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-age_clusters.json
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-age_clusters.tsv
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-age_design.png
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-age_stat-effect_statmap.nii.gz
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-age_stat-p_statmap.nii.gz
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-age_stat-t_statmap.nii.gz
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-age_stat-variance_statmap.nii.gz
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-age_stat-z_statmap.nii.gz
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-sex_clusters.json
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-sex_clusters.tsv
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-sex_design.png
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-sex_stat-effect_statmap.nii.gz
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-sex_stat-p_statmap.nii.gz
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-sex_stat-t_statmap.nii.gz
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-sex_stat-variance_statmap.nii.gz
derivatives/nilearn_glm/group/ageEffectOnGM_contrast-sex_stat-z_statmap.nii.gz
derivatives/nilearn_glm/group/ageEffectOnGM_corrdesign.png
derivatives/nilearn_glm/group/ageEffectOnGM_design.png
derivatives/nilearn_glm/group/ageEffectOnGM_design.tsv
derivatives/nilearn_glm/group/ageEffectOnGM_mask.nii.gz
derivatives/nilearn_glm/group/ageEffectOnGM_report.html
derivatives/nilearn_glm/group/ageEffectOnGM_stat-errorts_statmap.nii.gz
derivatives/nilearn_glm/group/ageEffectOnGM_stat-rsquared_statmap.nii.gz
derivatives/nilearn_glm/group/ageEffectOnGM_statmap.json
[SecondLevelModel.generate_report] Generating contrast-level figures...
[SecondLevelModel.generate_report] Generating design matrices figures...
[SecondLevelModel.generate_report] Generating contrast matrices figures...

Statistical Report - Second Level Model Implement the :term:`General Linear Model<GLM>` for multiple subject :term:`fMRI` data.

Description

Data were analyzed using Nilearn (version= 0.11.2.dev324+gc95672634; RRID:SCR_001362).

At the group level, a mass univariate analysis was performed with a linear regression at each voxel of the brain.

Input images were smoothed with gaussian kernel (full-width at half maximum=2.0 mm).

Model details

Value
Parameter
smoothing_fwhm (mm) 2.0

Design Matrix

run 0

Plot of design matrix for 0.

correlation matrix

Plot of correlation of design matrix for run 0.

Contrasts

Plot of the contrast age (run 0).
Plot of the contrast sex (run 0).

Mask

Mask image

The mask includes 191002 voxels (21.2 %) of the image.

Statistical Maps

age

Stat map plot for the contrast: age
Cluster Table
Height control None
Threshold Z 3.09
Cluster size threshold (voxels) 0
Minimum distance (mm) 8.0
Cluster ID X Y Z Peak Stat Cluster Size (mm3)

sex

Stat map plot for the contrast: sex
Cluster Table
Height control None
Threshold Z 3.09
Cluster size threshold (voxels) 0
Minimum distance (mm) 8.0
Cluster ID X Y Z Peak Stat Cluster Size (mm3)

About

  • Date preprocessed:


Total running time of the script: (1 minutes 9.420 seconds)

Estimated memory usage: 1732 MB

Gallery generated by Sphinx-Gallery