8.3.4. Voxel-Based Morphometry on Oasis dataset with Space-Net prior

Predicting age from gray-matter concentration maps from OASIS dataset. Note that age is a continuous variable, we use the regressor here, and not the classification object.

See also the SpaceNet documentation: SpaceNet: decoding with spatial structure for better maps.

# Authors: DOHMATOB Elvis
#          FRITSCH Virgile

8.3.4.1. Load the Oasis VBM dataset

import numpy as np
from nilearn import datasets
n_subjects = 200  # increase this number if you have more RAM on your box
dataset_files = datasets.fetch_oasis_vbm(n_subjects=n_subjects)
age = dataset_files.ext_vars['age'].astype(float)
age = np.array(age)
gm_imgs = np.array(dataset_files.gray_matter_maps)


# Split data into training set and test set
from sklearn.utils import check_random_state
from sklearn.cross_validation import train_test_split
rng = check_random_state(42)
gm_imgs_train, gm_imgs_test, age_train, age_test = train_test_split(
    gm_imgs, age, train_size=.6, random_state=rng)

# Sort test data for better visualization (trend, etc.)
perm = np.argsort(age_test)[::-1]
age_test = age_test[perm]
gm_imgs_test = gm_imgs_test[perm]

8.3.4.2. Fit the SpaceNet and predict with it

from nilearn.decoding import SpaceNetRegressor

# To save time (because these are anat images with many voxels), we include
# only the 5-percent voxels most correlated with the age variable to fit.
# Also, we set memory_level=2 so that more of the intermediate computations
# are cached. Also, you may pass and n_jobs=<some_high_value> to the
# SpaceNetRegressor class, to take advantage of a multi-core system.
#
# Also, here we use a graph-net penalty but more beautiful results can be
# obtained using the TV-l1 penalty, at the expense of longer runtimes.
decoder = SpaceNetRegressor(memory="nilearn_cache", penalty="graph-net",
                            screening_percentile=5., memory_level=2)
decoder.fit(gm_imgs_train, age_train)  # fit
coef_img = decoder.coef_img_
y_pred = decoder.predict(gm_imgs_test).ravel()  # predict
mse = np.mean(np.abs(age_test - y_pred))
print('Mean square error (MSE) on the predicted age: %.2f' % mse)

Out:

[NiftiMasker.fit] Loading data from [/home/varoquau/nilearn_data/oasis1/OAS1_0003_MR1/mwrc1OAS1_0003_MR1_mpr_anon_fslswapdim_bet.nii.gz, /home/varoquau/nilearn_data/oasis1/OAS1_0086_MR1/mwrc1OAS1_0086_MR1_mpr_anon_fslswapdim_bet.nii.gz,
[NiftiMasker.fit] Computing the mask
[NiftiMasker.fit] Resampling mask
[NiftiMasker.transform_single_imgs] Loading data from Nifti1Image(
shape=(91, 109, 91, 120),
affine=array([[  -2.,    0.,    0.,   90.],
       [   0.,    2.,    0., -126.],
       [   0.,    0.,    2.,  -72.],
       [   0.,    0.,    0.,    1.]])
)
[NiftiMasker.transform_single_imgs] Extracting region signals
[NiftiMasker.transform_single_imgs] Cleaning extracted signals
________________________________________________________________________________
[Memory] Calling nilearn.decoding.space_net.path_scores...
path_scores(<function _graph_net_squared_loss at 0x7fbcd366d2a8>, array([[ 0., ...,  0.],
       ...,
       [ 0., ...,  0.]]), array([ 73., ...,  62.]), array([[[ True, ...,  True],
        ...,
        [ True, ...,  True]],

       ...,
       [[ True, ...,  True],
        ...,
        [ True, ...,  True]]], dtype=bool),
None, [0.5], array([ 15, ..., 119]), array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14]), { 'max_iter': 1000, 'tol': 0.0001}, debias=False, verbose=1, n_alphas=10, eps=0.001, key=(0, 0), screening_percentile=1.2652228933482088, is_classif=False)
_____________________________________________________path_scores - 39.2s, 0.7min
________________________________________________________________________________
[Memory] Calling nilearn.decoding.space_net.path_scores...
path_scores(<function _graph_net_squared_loss at 0x7fbcd366d2a8>, array([[ 0., ...,  0.],
       ...,
       [ 0., ...,  0.]]), array([ 73., ...,  62.]), array([[[ True, ...,  True],
        ...,
        [ True, ...,  True]],

       ...,
       [[ True, ...,  True],
        ...,
        [ True, ...,  True]]], dtype=bool),
None, [0.5], array([  0, ..., 119]), array([15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]), { 'max_iter': 1000, 'tol': 0.0001}, debias=False, verbose=1, n_alphas=10, eps=0.001, key=(0, 1), screening_percentile=1.2652228933482088, is_classif=False)
_____________________________________________________path_scores - 26.4s, 0.4min
________________________________________________________________________________
[Memory] Calling nilearn.decoding.space_net.path_scores...
path_scores(<function _graph_net_squared_loss at 0x7fbcd366d2a8>, array([[ 0., ...,  0.],
       ...,
       [ 0., ...,  0.]]), array([ 73., ...,  62.]), array([[[ True, ...,  True],
        ...,
        [ True, ...,  True]],

       ...,
       [[ True, ...,  True],
        ...,
        [ True, ...,  True]]], dtype=bool),
None, [0.5], array([  0, ..., 119]), array([30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44]), { 'max_iter': 1000, 'tol': 0.0001}, debias=False, verbose=1, n_alphas=10, eps=0.001, key=(0, 2), screening_percentile=1.2652228933482088, is_classif=False)
_____________________________________________________path_scores - 39.2s, 0.7min
________________________________________________________________________________
[Memory] Calling nilearn.decoding.space_net.path_scores...
path_scores(<function _graph_net_squared_loss at 0x7fbcd366d2a8>, array([[ 0., ...,  0.],
       ...,
       [ 0., ...,  0.]]), array([ 73., ...,  62.]), array([[[ True, ...,  True],
        ...,
        [ True, ...,  True]],

       ...,
       [[ True, ...,  True],
        ...,
        [ True, ...,  True]]], dtype=bool),
None, [0.5], array([  0, ..., 119]), array([45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59]), { 'max_iter': 1000, 'tol': 0.0001}, debias=False, verbose=1, n_alphas=10, eps=0.001, key=(0, 3), screening_percentile=1.2652228933482088, is_classif=False)
_____________________________________________________path_scores - 39.1s, 0.7min
________________________________________________________________________________
[Memory] Calling nilearn.decoding.space_net.path_scores...
path_scores(<function _graph_net_squared_loss at 0x7fbcd366d2a8>, array([[ 0., ...,  0.],
       ...,
       [ 0., ...,  0.]]), array([ 73., ...,  62.]), array([[[ True, ...,  True],
        ...,
        [ True, ...,  True]],

       ...,
       [[ True, ...,  True],
        ...,
        [ True, ...,  True]]], dtype=bool),
None, [0.5], array([  0, ..., 119]), array([60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74]), { 'max_iter': 1000, 'tol': 0.0001}, debias=False, verbose=1, n_alphas=10, eps=0.001, key=(0, 4), screening_percentile=1.2652228933482088, is_classif=False)
_____________________________________________________path_scores - 39.2s, 0.7min
________________________________________________________________________________
[Memory] Calling nilearn.decoding.space_net.path_scores...
path_scores(<function _graph_net_squared_loss at 0x7fbcd366d2a8>, array([[ 0., ...,  0.],
       ...,
       [ 0., ...,  0.]]), array([ 73., ...,  62.]), array([[[ True, ...,  True],
        ...,
        [ True, ...,  True]],

       ...,
       [[ True, ...,  True],
        ...,
        [ True, ...,  True]]], dtype=bool),
None, [0.5], array([  0, ..., 119]), array([75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89]), { 'max_iter': 1000, 'tol': 0.0001}, debias=False, verbose=1, n_alphas=10, eps=0.001, key=(0, 5), screening_percentile=1.2652228933482088, is_classif=False)
_____________________________________________________path_scores - 32.5s, 0.5min
________________________________________________________________________________
[Memory] Calling nilearn.decoding.space_net.path_scores...
path_scores(<function _graph_net_squared_loss at 0x7fbcd366d2a8>, array([[ 0., ...,  0.],
       ...,
       [ 0., ...,  0.]]), array([ 73., ...,  62.]), array([[[ True, ...,  True],
        ...,
        [ True, ...,  True]],

       ...,
       [[ True, ...,  True],
        ...,
        [ True, ...,  True]]], dtype=bool),
None, [0.5], array([  0, ..., 119]), array([ 90,  91,  92,  93,  94,  95,  96,  97,  98,  99, 100, 101, 102,
       103, 104]),
{ 'max_iter': 1000, 'tol': 0.0001}, debias=False, verbose=1, n_alphas=10, eps=0.001, key=(0, 6), screening_percentile=1.2652228933482088, is_classif=False)
_____________________________________________________path_scores - 28.3s, 0.5min
________________________________________________________________________________
[Memory] Calling nilearn.decoding.space_net.path_scores...
path_scores(<function _graph_net_squared_loss at 0x7fbcd366d2a8>, array([[ 0., ...,  0.],
       ...,
       [ 0., ...,  0.]]), array([ 73., ...,  62.]), array([[[ True, ...,  True],
        ...,
        [ True, ...,  True]],

       ...,
       [[ True, ...,  True],
        ...,
        [ True, ...,  True]]], dtype=bool),
None, [0.5], array([  0, ..., 104]), array([105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117,
       118, 119]),
{ 'max_iter': 1000, 'tol': 0.0001}, debias=False, verbose=1, n_alphas=10, eps=0.001, key=(0, 7), screening_percentile=1.2652228933482088, is_classif=False)
_____________________________________________________path_scores - 39.2s, 0.7min
Time Elapsed: 310.384 seconds, 5 minutes.
[NiftiMasker.transform_single_imgs] Loading data from Nifti1Image(
shape=(91, 109, 91, 80),
affine=array([[  -2.,    0.,    0.,   90.],
       [   0.,    2.,    0., -126.],
       [   0.,    0.,    2.,  -72.],
       [   0.,    0.,    0.,    1.]])
)
[NiftiMasker.transform_single_imgs] Extracting region signals
[NiftiMasker.transform_single_imgs] Cleaning extracted signals
Mean square error (MSE) on the predicted age: 12.49

Visualize the resulting maps

from nilearn.plotting import plot_stat_map
# weights map
background_img = gm_imgs[0]
plot_stat_map(coef_img, background_img, title="graph-net weights",
              display_mode="z", cut_coords=1)
../../_images/sphx_glr_plot_oasis_vbm_space_net_001.png

8.3.4.3. Visualize the quality of predictions

import matplotlib.pyplot as plt
plt.figure()
plt.suptitle("graph-net: Mean Absolute Error %.2f years" % mse)
linewidth = 3
ax1 = plt.subplot('211')
ax1.plot(age_test, label="True age", linewidth=linewidth)
ax1.plot(y_pred, '--', c="g", label="Predicted age", linewidth=linewidth)
ax1.set_ylabel("age")
plt.legend(loc="best")
ax2 = plt.subplot("212")
ax2.plot(age_test - y_pred, label="True age - predicted age",
         linewidth=linewidth)
ax2.set_xlabel("subject")
plt.legend(loc="best")

plt.show()
../../_images/sphx_glr_plot_oasis_vbm_space_net_002.png

Total running time of the script: ( 5 minutes 18.978 seconds)

Generated by Sphinx-Gallery