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 and gray matter density.

The data comes from the Open Access Series of Imaging Studies (OASIS) project. If you use it, you need to agree with the data usage agreement available on the website.

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

See also

For more information see the dataset description.

Predictive modeling analysis

We run a Support Vector Regression (SVR) and ANOVA using the nilearn DecoderRegressor to predict age from the VBM data. We use a subset of the subjects from the OASIS dataset to limit the memory usage.

Important

Note that for an actual predictive modeling study of aging, the study should be ran on the full set of subjects.

Also, all parameters should be selected by nested cross-validation. Indeed, things like the smoothing applied to the data and the number of features selected by the ANOVA step can impact significantly the prediction score.

# Use a single variable to control the verbosity of the script.
verbose = 1

# Several of Nilearn's estimators (like the DecoderRegressor we use here)
# accept a ``n_jobs=<some_high_value>``
# to take advantage of a multi-core system.
n_jobs = 2

Load Oasis dataset

We fetch the data and split it into training set and test set.

from sklearn.model_selection import train_test_split

from nilearn.datasets import fetch_oasis_vbm

n_subjects = 200  # more subjects require more memory

oasis_dataset = fetch_oasis_vbm(n_subjects=n_subjects, verbose=verbose)

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

gray_matter_map_filenames = oasis_dataset.gray_matter_maps
age = oasis_dataset.ext_vars["age"].to_numpy()
gm_imgs_train, gm_imgs_test, age_train, age_test = train_test_split(
    gray_matter_map_filenames, age, train_size=0.6, random_state=0
)

[fetch_oasis_vbm] Dataset directory found: /home/runner/work/nilearn/nilearn/nilearn_data/oasis1
First gray-matter anatomy image (3D) is located at: /home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0001_MR1/mwrc1OAS1_0001_MR1_mpr_anon_fslswapdim_bet.nii.gz

Preprocess data

The voxels with too low between-subject variance are removed using sklearn.feature_selection.VarianceThreshold.

Then we convert the data back to the mask image in order to use it for decoding process.

from sklearn.feature_selection import VarianceThreshold

from nilearn.maskers import NiftiMasker

nifti_masker = NiftiMasker(
    standardize=None,
    smoothing_fwhm=2,
    memory="nilearn_cache",  # cache options
    verbose=verbose,
)
gm_maps_masked = nifti_masker.fit_transform(gm_imgs_train)

variance_threshold = VarianceThreshold(threshold=0.01)
variance_threshold.fit_transform(gm_maps_masked)

mask = nifti_masker.inverse_transform(variance_threshold.get_support())

[NiftiMasker.wrapped] Loading data from [
 mwrc1OAS1_0211_MR1...,
         ...
 mwrc1OAS1_0195_MR1...,
]
[NiftiMasker.wrapped] Computing mask
________________________________________________________________________________
[Memory] Calling nilearn.masking.compute_background_mask...
compute_background_mask([ '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0211_MR1/mwrc1OAS1_0211_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0052_MR1/mwrc1OAS1_0052_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0155_MR1/mwrc1OAS1_0155_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0221_MR1/mwrc1OAS1_0221_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0121_MR1/mwrc1OAS1_0121_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0070_MR1/mwrc1OAS1_0070..., verbose=0)
__________________________________________compute_background_mask - 1.5s, 0.0min
[NiftiMasker.wrapped] Resampling mask
________________________________________________________________________________
[Memory] Calling nilearn.image.resampling.resample_img...
resample_img(<nibabel.nifti1.Nifti1Image object at 0x7fa1b05df910>, target_affine=None, target_shape=None, copy=False, interpolation='nearest')
_____________________________________________________resample_img - 0.0s, 0.0min
[NiftiMasker.wrapped] Finished fit
________________________________________________________________________________
[Memory] Calling nilearn.maskers.nifti_masker.filter_and_mask...
filter_and_mask([ '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0211_MR1/mwrc1OAS1_0211_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0052_MR1/mwrc1OAS1_0052_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0155_MR1/mwrc1OAS1_0155_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0221_MR1/mwrc1OAS1_0221_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0121_MR1/mwrc1OAS1_0121_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0070_MR1/mwrc1OAS1_0070...,
<nibabel.nifti1.Nifti1Image object at 0x7fa1b05df910>, { 'clean_args': None,
  'clean_kwargs': {},
  'cmap': 'gray',
  'detrend': False,
  'dtype': None,
  'high_pass': None,
  'high_variance_confounds': False,
  'low_pass': None,
  'reports': True,
  'runs': None,
  'smoothing_fwhm': 2,
  'standardize': None,
  'standardize_confounds': True,
  't_r': None,
  'target_affine': None,
  'target_shape': None}, memory_level=1, memory=Memory(location=nilearn_cache/joblib), verbose=1, confounds=None, sample_mask=None, copy=True, sklearn_output_config=None)
[NiftiMasker.wrapped] Loading data from <nibabel.nifti1.Nifti1Image object at 0x7fa1d4a85a10>
[NiftiMasker.wrapped] Smoothing images
[NiftiMasker.wrapped] Extracting region signals
[NiftiMasker.wrapped] Cleaning extracted signals
__________________________________________________filter_and_mask - 5.3s, 0.1min
/home/runner/work/nilearn/nilearn/examples/02_decoding/plot_oasis_vbm.py:98: UserWarning:

Casting boolean input to <class 'numpy.int32'>

[NiftiMasker.inverse_transform] Computing image from signals
________________________________________________________________________________
[Memory] Calling nilearn.masking.unmask...
unmask(array([0, ..., 0], dtype=int32), <nibabel.nifti1.Nifti1Image object at 0x7fa1b05df910>)
___________________________________________________________unmask - 0.0s, 0.0min

Prediction pipeline with ANOVA and SVR using DecoderRegressor

In Nilearn we can benefit from the built-in DecoderRegressor object to do an ANOVA with SVR instead of manually defining the whole pipeline.

This estimator also uses cross validation to select best models and ensemble them.

To save time (because these are anat images with many voxels), we include only the 1-percent of the voxels most correlated with the age variable to fit.

We also want to set mask hyperparameter to be the mask we just obtained above.

We then fit and predict with the decoder and sort test data for better visualization (trend, etc.).

import numpy as np

from nilearn.decoding import DecoderRegressor

decoder = DecoderRegressor(
    mask=mask,
    scoring="neg_mean_absolute_error",
    screening_percentile=1,
    n_jobs=n_jobs,
    verbose=verbose,
)
decoder.fit(gm_imgs_train, age_train)

perm = np.argsort(age_test)[::-1]
age_test = age_test[perm]
gm_imgs_test = np.array(gm_imgs_test)[perm]
age_pred = decoder.predict(gm_imgs_test)

prediction_score = -np.mean(decoder.cv_scores_["beta"])

print(f"explained variance for the cross-validation: {prediction_score:f}")

[DecoderRegressor.fit] Mask volume = 1.21634e+06mm^3 = 1216.34cm^3
[DecoderRegressor.fit] Standard brain volume = 1.88299e+06mm^3
[DecoderRegressor.fit] Original screening-percentile: 1
[DecoderRegressor.fit] Corrected screening-percentile: 1.54807
[DecoderRegressor.fit] The decoding model will be trained on 1520 features.
[DecoderRegressor.fit] The decoding model will be trained on 1520 features.
[Parallel(n_jobs=2)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=2)]: Done  10 out of  10 | elapsed:    3.0s finished
explained variance for the cross-validation: 11.774939

Visualization

import matplotlib.pyplot as plt

from nilearn.image import mean_img
from nilearn.plotting import plot_stat_map, show

bg_filename = mean_img(gray_matter_map_filenames)
z_slice = 0

weight_img = decoder.coef_img_["beta"]

display = plot_stat_map(
    weight_img,
    title="SVM weights",
    bg_img=bg_filename,
    display_mode="z",
    cut_coords=[z_slice],
    figure=plt.figure(figsize=(5.5, 7.5), facecolor="k"),
)
show()
plot oasis vbm

Visualize the quality of predictions

linewidth = 3

plt.figure(figsize=(6, 4.5))
plt.suptitle(f"Decoder: Mean Absolute Error {prediction_score:.2f} years")
plt.plot(age_test, label="True age", linewidth=linewidth)
plt.plot(age_pred, "--", c="g", label="Predicted age", linewidth=linewidth)
plt.ylabel("age")
plt.xlabel("subject")
plt.legend(loc="best")

plt.figure(figsize=(6, 4.5))
plt.plot(
    age_test - age_pred, label="True age - predicted age", linewidth=linewidth
)
plt.xlabel("subject")
plt.legend(loc="best")

show()
  • Decoder: Mean Absolute Error 11.77 years
  • plot oasis vbm

Brain mapping with mass univariate

SVM weights are very noisy, partly because heavy smoothing is detrimental for the prediction here. A standard analysis using mass-univariate GLM (here permuted to have exact correction for multiple comparisons) gives a much clearer view of the important regions.

from nilearn.image import get_data
from nilearn.mass_univariate import permuted_ols

gm_maps_masked = NiftiMasker(
    standardize=None,
    memory="nilearn_cache",  # cache options
    verbose=verbose,
).fit_transform(gray_matter_map_filenames)
data = variance_threshold.fit_transform(gm_maps_masked)

output = permuted_ols(
    age,
    data,  # + intercept as a covariate by default
    # few permutations in the interest of time; 10000 would be better
    n_perm=2000,
    verbose=verbose,
    n_jobs=n_jobs,
    random_state=0,  # to ensure reproducible results
)

[NiftiMasker.wrapped] Loading data from [
 mwrc1OAS1_0001_MR1...,
         ...
 mwrc1OAS1_0227_MR1...,
]
[NiftiMasker.wrapped] Computing mask
________________________________________________________________________________
[Memory] Calling nilearn.masking.compute_background_mask...
compute_background_mask([ '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0001_MR1/mwrc1OAS1_0001_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0002_MR1/mwrc1OAS1_0002_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0003_MR1/mwrc1OAS1_0003_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0004_MR1/mwrc1OAS1_0004_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0005_MR1/mwrc1OAS1_0005_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0006_MR1/mwrc1OAS1_0006..., verbose=0)
__________________________________________compute_background_mask - 2.4s, 0.0min
[NiftiMasker.wrapped] Resampling mask
[NiftiMasker.wrapped] Finished fit
________________________________________________________________________________
[Memory] Calling nilearn.maskers.nifti_masker.filter_and_mask...
filter_and_mask([ '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0001_MR1/mwrc1OAS1_0001_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0002_MR1/mwrc1OAS1_0002_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0003_MR1/mwrc1OAS1_0003_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0004_MR1/mwrc1OAS1_0004_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0005_MR1/mwrc1OAS1_0005_MR1_mpr_anon_fslswapdim_bet.nii.gz',
  '/home/runner/work/nilearn/nilearn/nilearn_data/oasis1/OAS1_0006_MR1/mwrc1OAS1_0006...,
<nibabel.nifti1.Nifti1Image object at 0x7fa1d4d14310>, { 'clean_args': None,
  'clean_kwargs': {},
  'cmap': 'gray',
  'detrend': False,
  'dtype': None,
  'high_pass': None,
  'high_variance_confounds': False,
  'low_pass': None,
  'reports': True,
  'runs': None,
  'smoothing_fwhm': None,
  'standardize': None,
  'standardize_confounds': True,
  't_r': None,
  'target_affine': None,
  'target_shape': None}, memory_level=1, memory=Memory(location=nilearn_cache/joblib), verbose=1, confounds=None, sample_mask=None, copy=True, sklearn_output_config=None)
[NiftiMasker.wrapped] Loading data from <nibabel.nifti1.Nifti1Image object at 0x7fa1814065d0>
[NiftiMasker.wrapped] Extracting region signals
[NiftiMasker.wrapped] Cleaning extracted signals
__________________________________________________filter_and_mask - 7.0s, 0.1min
[Parallel(n_jobs=2)]: Using backend LokyBackend with 2 concurrent workers.
[Parallel(n_jobs=2)]: Done   2 out of   2 | elapsed:   34.5s finished

Show results

neg_log_pvals = output["logp_max_t"]
t_scores_original_data = output["t"]
signed_neg_log_pvals = neg_log_pvals * np.sign(t_scores_original_data)
signed_neg_log_pvals_unmasked = nifti_masker.inverse_transform(
    variance_threshold.inverse_transform(signed_neg_log_pvals)
)
threshold = -np.log10(0.1)  # 10% corrected

n_detections = (get_data(signed_neg_log_pvals_unmasked) > threshold).sum()

title = (
    "Negative $\\log_{10}$ p-values"
    "\n(Non-parametric + max-type correction)"
    f"\n{int(n_detections)} detections"
)

plot_stat_map(
    signed_neg_log_pvals_unmasked,
    threshold=threshold,
    title=title,
    bg_img=bg_filename,
    display_mode="z",
    cut_coords=[z_slice],
    figure=plt.figure(figsize=(5.5, 7.5), facecolor="k"),
)
show()
plot oasis vbm
[NiftiMasker.inverse_transform] Computing image from signals
________________________________________________________________________________
[Memory] Calling nilearn.masking.unmask...
unmask(array([[0., ..., 0.]]), <nibabel.nifti1.Nifti1Image object at 0x7fa1b05df910>)
___________________________________________________________unmask - 0.1s, 0.0min

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

Estimated memory usage: 5948 MB

Gallery generated by Sphinx-Gallery