brain_age¶

Predict age from brain grey matter (regression). Aging is associated with is grey matter (GM) atrophy, each year, an adult lose 0.1% of GM. We will try to learn a predictor of the chronological age (true age) using GM measurements on the brain on a population of healthy control participants.

Such a predictor provides the expected brain age of a subject. Deviation from this expected brain age indicates acceleration or slowdown of the aging process which may be associated with a pathological neurobiological process or protective factor of aging.

Dataset¶

There are 357 samples in the training set and 90 samples in the test set.

Input data¶

Voxel-based_morphometry VBM using cat12 software which provides:

  • Regions Of Interest (rois) of Grey Matter (GM) scaled for the Total Intracranial Volume (TIV): [train|test]_rois.csv 284 features.

  • VBM GM 3D maps or images (vbm3d) of voxels in the MNI space: [train|test]_vbm.npz contains 3D images of shapes (121, 145, 121). This npz contains the 3D mask and the affine transformation to MNI referential. Masking the brain provide flat 331 695 input features (voxels) for each participant.

By default problem.get_[train|test]_data() return the concatenation of 284 ROIs of Grey Matter (GM) features with 331 695 features (voxels) within a brain mask. Those two blocks are higly redundant. To select only on ROIs (rois) features do:

X[:, :284]

To select only on (vbm) (voxel with the brain) features do:

X[:, 284:]

Target¶

The target can be found in [test|train]_participants.csv files, selecting the age column for regression problem.

Evaluation metrics¶

sklearn metrics

The main Evaluation metrics is the Root-mean-square deviation RMSE. We will also look at the R-squared R2.

Links¶

  • RAMP-workflow’s documentation
  • RAMP-workflow’s github
  • RAMP Kits

Installation¶

This starting kit requires Python and the following dependencies:

  • numpy
  • scipy
  • pandas
  • scikit-learn
  • matplolib
  • seaborn
  • jupyter
  • ramp-workflow

Therefore, we advise you to install Anaconda distribution which include almost all dependencies.

Only nilearn and ramp-workflow are not included by default in the Anaconda distribution. They will be installed from the execution of the notebook.

To run a submission and the notebook you will need the dependencies listed in requirements.txt. You can install the dependencies with the following command-line:

pip install -U -r requirements.txt

If you are using conda, we provide an environment.yml file for similar usage.

conda env create -f environment.yml

Then, you can activate the environment using:

conda activate brain_age

And desactivate using

conda deactivate

Getting started¶

  1. Download the data locally:
python download_data.py
  1. Execute the jupyter notebook, from the root directory using:
jupyter notebook brain_age_starting_kit.ipynb

Tune your model using the starting_kit

  1. Submission (Run locally)

The submissions need to be located in the submissions folder. For instance for linear_regression_rois, it should be located in submissions/submissions/linear_regression_rois.

Copy everything required to build your estimator in a submission file: submissions/submissions/linear_regression_rois/estimator.py. This file must contain a function get_estimator().

Run locally:

ramp-test --submission linear_regression_rois
  1. Submission on RAMP:

Using RAMP starting-kits

Descriptive statistics¶

In [1]:
%matplotlib inline
import os
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

participants_train = pd.read_csv(os.path.join("data", "train_participants.csv" ))
participants_train["set"] = 'train'
participants_test = pd.read_csv(os.path.join("data", "test_participants.csv" ))
participants_test["set"] = 'test'
participants = participants_train.append(participants_test)

sns.violinplot(x="set", y="age", data=participants)
participants[["age", "set"]].groupby("set").describe()
Out[1]:
age
count mean std min 25% 50% 75% max
set
test 90.0 47.848019 17.852717 20.071184 31.819986 44.855578 62.017112 82.187543
train 357.0 49.138846 16.095719 19.980835 35.455168 50.902122 62.061602 86.318960

Unsupervized analysis: Explore main sources of variabilities¶

Load ROIs data

In [2]:
globvol_rois_train = pd.read_csv(os.path.join("data", "train_rois.csv"))
print(globvol_rois_train.iloc[:5, :10])

   participant_id     CSF_Vol      GM_Vol      WM_Vol  l3thVen_GM_Vol  \
0             651  334.031002  603.083565  560.465707        0.046239   
1             431  300.991328  603.101354  594.425495        0.052599   
2             398  252.410342  641.965684  604.320324        0.063052   
3             419  402.533442  604.978032  490.709738        0.043962   
4             627  409.460047  579.029937  509.311974        0.042226   

   r3thVen_GM_Vol  l4thVen_GM_Vol  r4thVen_GM_Vol  lAcc_GM_Vol  rAcc_GM_Vol  
0        0.052240        0.063827        0.068097     0.419533     0.427702  
1        0.047901        0.088186        0.105241     0.373451     0.377817  
2        0.051575        0.083113        0.076834     0.447264     0.465415  
3        0.041684        0.089264        0.083791     0.391540     0.388656  
4        0.036072        0.073829        0.081725     0.339043     0.358513

train_rois.csv provides:

  • Global volumes of "tissues": CerrebroSpinal Fluid (CSF_Vol), Grey (GM_Vol) and White Matter (WM_Vol) volume of participants.
  • ROIs are starting at column l3thVen_GM_Vol. Note that rois_train.loc[:, 'l3thVen_GM_Vol':] matches problem.get_train_data()[:, :284].
In [3]:
rois_train = globvol_rois_train.loc[:, 'l3thVen_GM_Vol':]
vols_train = globvol_rois_train.loc[:, ['CSF_Vol', 'GM_Vol', 'WM_Vol']]

PCA on ROIs: explore global effect of age¶

In [4]:
from sklearn.decomposition import PCA

pca_rois = PCA(n_components=2)
PCs_rois = pca_rois.fit_transform(rois_train)
print(pca_rois.explained_variance_ratio_)

df = pd.DataFrame(dict(age=participants_train['age'], PC1_ROIs=PCs_rois[:, 0], PC2_ROIs=PCs_rois[:, 1]))

[0.54771351 0.10711543]

Main sources of variability ?

In [5]:
sns.pairplot(df)
print(df.corr())

               age      PC1_ROIs      PC2_ROIs
age       1.000000  8.177015e-01 -1.108805e-01
PC1_ROIs  0.817702  1.000000e+00 -1.133611e-17
PC2_ROIs -0.110881 -1.133611e-17  1.000000e+00

Neurobiological effect of age¶

PC1_ROIs is higlhy correlated with age Thus age seems to be the main source of variability.

Question: what is Neurobiological effect of age ? Lets explore the hypothesis of brain atrophy.

We compute brain atrophy indices:

  • GM_ratio: GM_Vol / (GM_Vol + WM_Vol + CSF_Vol). It is ratio of GM in the brain.
  • WM_ratio: WM_Vol / (GM_Vol + WM_Vol + CSF_Vol). It is ratio of WM in the brain.
In [6]:
df["GM_ratio"] = vols_train.loc[:, "GM_Vol"] / vols_train.sum(axis=1)
df["WM_ratio"] = vols_train.loc[:, "WM_Vol"] / vols_train.sum(axis=1)
df["CSF_ratio"] = vols_train.loc[:, "CSF_Vol"] / vols_train.sum(axis=1)

sns.pairplot(df[["PC1_ROIs", "age", "GM_ratio", "WM_ratio","CSF_ratio"]])
print(df[["PC1_ROIs", "age", "GM_ratio", "WM_ratio", "CSF_ratio"]].corr())

           PC1_ROIs       age  GM_ratio  WM_ratio  CSF_ratio
PC1_ROIs   1.000000  0.817702 -0.944958 -0.467386   0.960228
age        0.817702  1.000000 -0.818380 -0.345080   0.803764
GM_ratio  -0.944958 -0.818380  1.000000  0.225909  -0.890849
WM_ratio  -0.467386 -0.345080  0.225909  1.000000  -0.643805
CSF_ratio  0.960228  0.803764 -0.890849 -0.643805   1.000000

Very strong effect on GM:

  • GM_ratio correlate with PC1_ROIs
  • GM_ratio correlate with age

Effect on WM:

  • WM_ratio correlate with PC1_ROIs
  • WM_ratio correlate with age

Aging (strongly correlates with PC1) and leads to decrease of GM_ratio and increase of CSF_ratio:

Conclusion¶

Aging leads to brain atrophy with is the main source of brain variability. Therefore it should be possible to predict the age form the brain anatomy.

Machine learning¶

Import and read data

In [7]:
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Ridge
from sklearn.base import BaseEstimator
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import cross_validate
import sklearn.metrics as metrics
import problem

from sklearn.base import BaseEstimator
from sklearn.base import TransformerMixin
In [8]:
X_train, y_train = problem.get_train_data()
X_test, y_test = problem.get_test_data()

assert X_train.shape[1] == 284 + 331695

Utils functions¶

Function to compute scores

In [9]:
def cv_train_test_scores(rmse_cv_test, rmse_cv_train, r2_cv_test, r2_cv_train,
                         y_train, y_pred_train, y_test, y_pred_test):
    """Compute CV score, train and test score from a cv grid search model.

    Parameters
    ----------
    rmse_cv_test : array
        Test rmse across CV folds.
    rmse_cv_train : array
        Train rmse across CV folds.

    r2_cv_test : array
        Test R2 across CV folds.
    r2_cv_train : array
        Train R2 across CV folds.

    y_train : array
        True train values.
    y_pred_train : array
        Predicted train values.
    y_test : array
        True test values.
    y_pred_test : array
        Predicted test values.

    Returns
    -------
    info : TYPE
        DataFrame(r2_cv, r2_train, mae_train, mse_train).
    """
    # CV scores
    rmse_cv_test_mean, rmse_cv_test_sd = np.mean(rmse_cv_test), np.std(rmse_cv_test)
    rmse_cv_train_mean, rmse_cv_train_sd = np.mean(rmse_cv_train), np.std(rmse_cv_train)

    r2_cv_test_mean, r2_cv_test_sd = np.mean(r2_cv_test), np.std(r2_cv_test)
    r2_cv_train_mean, r2_cv_train_sd = np.mean(r2_cv_train), np.std(r2_cv_train)

    # Test scores
    rmse_test = np.sqrt(metrics.mean_squared_error(y_test, y_pred_test))
    r2_test = metrics.r2_score(y_test, y_pred_test)

    # Train scores
    rmse_train = np.sqrt(metrics.mean_squared_error(y_train, y_pred_train))
    r2_train = metrics.r2_score(y_train, y_pred_train)

    
    scores = pd.DataFrame([[rmse_cv_test_mean, rmse_cv_test_sd, rmse_cv_train_mean, rmse_cv_train_sd,
                            r2_cv_test_mean, rmse_cv_test_sd, r2_cv_train_mean, r2_cv_train_sd,
                            rmse_test, r2_test, rmse_train, r2_train
                           ]],
                        columns=('rmse_cv_test_mean', 'rmse_cv_test_sd', 'rmse_cv_train_mean', 'rmse_cv_train_sd',
                                 'r2_cv_test_mean', 'rmse_cv_test_sd', 'r2_cv_train_mean', 'r2_cv_train_sd',
                                 'rmse_test', 'r2_test', 'rmse_train', 'r2_train'
                                 ))

    return scores

Feature extractor of ROIs or voxels within the brain (VBM)¶

Selecting only rois or vbm images:

This can be achieved by a ROIsFeatureExtractor or VBMFeatureExtractor

In [10]:
class ROIsFeatureExtractor(BaseEstimator, TransformerMixin):
    """Select only the 284 ROIs features:"""
    def fit(self, X, y):
        return self

    def transform(self, X):
        return X[:, :284]

class VBMFeatureExtractor(BaseEstimator, TransformerMixin):
    """Select only the 284 ROIs features:"""
    def fit(self, X, y):
        return self

    def transform(self, X):
        return X[:, 284:]


fe = ROIsFeatureExtractor()
print(fe.transform(X_train).shape)

fe = VBMFeatureExtractor()
print(fe.transform(X_train).shape)

(357, 284)
(357, 331695)

Design of predictors and their evaluation using CV and test set¶

The framework is evaluated with a cross-validation approach. The metrics used are the root-mean-square error (RMSE).

First we propose a simple Regression predictor based on ROIs features only:

In [11]:
cv = problem.get_cv(X_train, y_train)

estimator = make_pipeline(ROIsFeatureExtractor(), StandardScaler(), LinearRegression())

cv_results = cross_validate(estimator, X_train, y_train, scoring=['neg_root_mean_squared_error', 'r2'], cv=cv,
                         verbose=1, return_train_score=True, n_jobs=5)

# Refit on all train
estimator.fit(X_train, y_train)
# Apply on test
y_pred_train = estimator.predict(X_train)
y_pred_test = estimator.predict(X_test)

print("Important scores are rmse_cv_test_mean and rmse_test")
cv_train_test_scores(rmse_cv_test=-cv_results['test_neg_root_mean_squared_error'],
                     rmse_cv_train=-cv_results['train_neg_root_mean_squared_error'],
                     r2_cv_test=cv_results['test_r2'],
                     r2_cv_train=cv_results['train_r2'],
                     y_train=y_train, y_pred_train=y_pred_train, y_test=y_test, y_pred_test=y_pred_test).T.round(3)

[Parallel(n_jobs=5)]: Using backend LokyBackend with 5 concurrent workers.
[Parallel(n_jobs=5)]: Done   2 out of   5 | elapsed:    2.9s remaining:    4.4s

Important scores are rmse_cv_test_mean and rmse_test

[Parallel(n_jobs=5)]: Done   5 out of   5 | elapsed:    3.3s finished
Out[11]:
0
rmse_cv_test_mean 47.793
rmse_cv_test_sd 16.807
rmse_cv_train_mean 0.639
rmse_cv_train_sd 0.288
r2_cv_test_mean -8.540
rmse_cv_test_sd 16.807
r2_cv_train_mean 0.998
r2_cv_train_sd 0.002
rmse_test 14.350
r2_test 0.347
rmse_train 2.659
r2_train 0.973

Those are really bad results:

  • Using CV we can predict the age +- 47.8 years of error !!!.
  • Using independant test set we can predict the age +- 14.3 years of error !!!.

Then we test a simple Regression predictor based on vbm features:

In [12]:
estimator = make_pipeline(VBMFeatureExtractor(), StandardScaler(), LinearRegression())

cv = problem.get_cv(X_train, y_train)

cv_results = cross_validate(estimator, X_train, y_train, scoring=['neg_root_mean_squared_error', 'r2'], cv=cv,
                         verbose=1, return_train_score=True, n_jobs=5)

cv_train_test_scores(rmse_cv_test=-cv_results['test_neg_root_mean_squared_error'],
                     rmse_cv_train=-cv_results['train_neg_root_mean_squared_error'],
                     r2_cv_test=cv_results['test_r2'],
                     r2_cv_train=cv_results['train_r2'],
                     y_train=y_train, y_pred_train=y_pred_train, y_test=y_test, y_pred_test=y_pred_test).T.round(3)

[Parallel(n_jobs=5)]: Using backend LokyBackend with 5 concurrent workers.
[Parallel(n_jobs=5)]: Done   2 out of   5 | elapsed:   54.6s remaining:  1.4min
[Parallel(n_jobs=5)]: Done   5 out of   5 | elapsed:   54.9s finished
Out[12]:
0
rmse_cv_test_mean 6.968
rmse_cv_test_sd 0.623
rmse_cv_train_mean 0.000
rmse_cv_train_sd 0.000
r2_cv_test_mean 0.808
rmse_cv_test_sd 0.623
r2_cv_train_mean 1.000
r2_cv_train_sd 0.000
rmse_test 14.350
r2_test 0.347
rmse_train 2.659
r2_train 0.973

Using CV we can predict the age +- 7 years of error. But using an independant test set we can predict the age +- 14.3 years of error !!!. Your turn to do better !

Submission¶

The submissions need to be located in the submissions folder. For instance for linear_regression_rois, it should be located in submissions/submissions/linear_regression_rois.

Copy everything required (the cell bellow) to build your estimator in a submission file: submissions/submissions/linear_regression_rois/estimator.py. This file must contain a function get_estimator():

In [13]:
import numpy as np

from sklearn.base import BaseEstimator
from sklearn.base import TransformerMixin
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Ridge
from sklearn.base import BaseEstimator
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline


class ROIsFeatureExtractor(BaseEstimator, TransformerMixin):
    """Select only the 284 ROIs features:"""
    def fit(self, X, y):
        return self

    def transform(self, X):
        return X[:, :284]


def get_estimator():
    """Build your estimator here."""
    estimator = make_pipeline(ROIsFeatureExtractor(), StandardScaler(),
                              LinearRegression())
    return estimator

Run locally:

ramp-test --submission linear_regression_rois

Submission on RAMP:

Using RAMP starting-kits

In [ ]: