Nilearn_tutorial_2018_3_Machine_learning.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Section 3: Machine learning to predict age from rs-fmri\n",
    "\n",
    "We will integrate what we've learned in the previous sections to extract data from rs-fmri images, and use that data as features in a machine learning model\n",
    "\n",
    "The dataset consists of ~150 subjects, mostly young children and some young adults. We will use rs-fmri data to predict the age of the participants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Run the cell below to ready the download script (otherwise ignore it)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from sklearn.datasets.base import Bunch\n",
    "from nilearn.datasets.utils import _get_dataset_dir, _fetch_files\n",
    "\n",
    "\n",
    "def fetch_data(n_subjects=30, data_dir=None, url=None, resume=True,\n",
    "               verbose=1):\n",
    "    \"\"\"Download and load the dataset.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    n_subjects: int, optional\n",
    "        The number of subjects to load from maximum of 40 subjects.\n",
    "        By default, 30 subjects will be loaded. If None is given,\n",
    "        all 40 subjects will be loaded.\n",
    "\n",
    "    data_dir: string, optional\n",
    "        Path of the data directory. Used to force data storage in a specified\n",
    "        location. Default: None\n",
    "\n",
    "    url: string, optional\n",
    "        Override download URL. Used for test only (or if you setup a mirror of\n",
    "        the data). Default: None\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    data: sklearn.datasets.base.Bunch\n",
    "        Dictionary-like object, the interest attributes are :\n",
    "         - 'func': Paths to functional resting-state images\n",
    "         - 'phenotypic': Explanations of preprocessing steps\n",
    "         - 'confounds': CSV files containing the nuisance variables\n",
    "\n",
    "    References\n",
    "    ----------\n",
    "    :Download:\n",
    "        https://openneuro.org/datasets/ds000228/versions/00001\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    if url is None:\n",
    "        url = 'https://openneuro.org/crn/datasets/ds000228/snapshots/00001/files/'\n",
    "\n",
    "    # Preliminary checks and declarations\n",
    "    dataset_name = 'ds000228'\n",
    "    data_dir = _get_dataset_dir(dataset_name, data_dir=data_dir,\n",
    "                                verbose=verbose)\n",
    "    max_subjects = 155\n",
    "    if n_subjects is None:\n",
    "        n_subjects = max_subjects\n",
    "    if n_subjects > max_subjects:\n",
    "        warnings.warn('Warning: there are only %d subjects' % max_subjects)\n",
    "        n_subjects = max_subjects\n",
    "    ids = range(1, n_subjects + 1)\n",
    "\n",
    "    # First, get the metadata\n",
    "    phenotypic = (\n",
    "            'participants.tsv',\n",
    "            url + 'participants.tsv', dict())\n",
    "\n",
    "    phenotypic = _fetch_files(data_dir, [phenotypic], resume=resume,\n",
    "                              verbose=verbose)[0]\n",
    "\n",
    "    # Load the csv file\n",
    "    phenotypic = np.genfromtxt(phenotypic, names=True, delimiter='\\t',\n",
    "                               dtype=None)\n",
    "\n",
    "    # Keep phenotypic information for selected subjects\n",
    "    int_ids = np.asarray(ids, dtype=int)\n",
    "    phenotypic = phenotypic[[i - 1 for i in int_ids]]\n",
    "\n",
    "    # Download dataset files\n",
    "\n",
    "    functionals = [\n",
    "        'derivatives:fmriprep:sub-pixar%03i:sub-pixar%03i_task-pixar_run-001_swrf_bold.nii.gz' % (i, i)\n",
    "        for i in ids]\n",
    "    urls = [url + name for name in functionals]\n",
    "    functionals = _fetch_files(\n",
    "        data_dir, zip(functionals, urls, (dict(),) * n_subjects),\n",
    "        resume=resume, verbose=verbose)\n",
    "\n",
    "    confounds = [\n",
    "        'derivatives:fmriprep:sub-pixar%03i:sub-pixar%03i_task-pixar_run-001_ART_and_CompCor_nuisance_regressors.mat'\n",
    "        % (i, i)\n",
    "        for i in ids]\n",
    "    confound_urls = [url + name for name in confounds]\n",
    "\n",
    "    confounds = _fetch_files(\n",
    "        data_dir, zip(confounds, confound_urls, (dict(),) * n_subjects),\n",
    "        resume=resume, verbose=verbose)\n",
    "\n",
    "    return Bunch(func=functionals, confounds=confounds,\n",
    "                 phenotypic=phenotypic, description='ds000228')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Confounds function\n",
    "\n",
    "import h5py\n",
    "import numpy as np\n",
    "\n",
    "def prepare_confounds(conf, key = 'R', transpose=True):\n",
    "    arrays = {}\n",
    "    f = h5py.File(conf)\n",
    "    for k, v in f.items():\n",
    "        arrays[k] = np.array(v)\n",
    "    \n",
    "    if transpose:\n",
    "        output = arrays[key].T\n",
    "    else:\n",
    "        output = arrays[key]\n",
    "    \n",
    "    return output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load the data\n",
    "\n",
    "If you already have downloaded the data, the next commands should load it. If you have not, the data will begin downloading. This is a lot of data (more than 14 GB) so maybe not good to download right now!\n",
    "\n",
    "You can still pop in after the \"extract features\" section"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# if the data is already stored somewhere OR you want to store it somewhere, specify this here\n",
    "#wdir = '/Users/jakevogel/Science/Nilearn_tutorial/ds000028/'\n",
    "\n",
    "wdir = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Now fetch the data\n",
    "data = fetch_data(None,data_dir=wdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many individual subjects do we have?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "len(data.func)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get Y (our target) and assess its distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Let's load the phenotype data\n",
    "data.phenotypic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pandas import read_csv\n",
    "\n",
    "pheno = read_csv(data.phenotypic)\n",
    "#pheno = read_csv(data.phenotypic, sep='\\t')\n",
    "pheno.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like there is a column for age. Let's capture it in a variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "y_age = pheno['Age']\n",
    "y_age.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe we should have a look at the distribution of our target variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.hist(y_age)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Seems pretty skewed toward younger children. Perhaps we will have better results by log-transforming age?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "log_y_age = np.log(y_age)\n",
    "plt.hist(log_y_age)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a bit better. Maybe we will want to use log-transformed age in our models?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we are going to use the same techniques we learned in section 2 to extract rs-fmri connectivity features from every subject.\n",
    "\n",
    "How are we going to do that? With a for loop.\n",
    "\n",
    "Don't worry, it's not as scary as it sounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Here is a really simple for loop\n",
    "\n",
    "for i in range(10):\n",
    "    print('the number is', i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "container = []\n",
    "for i in range(10):\n",
    "    container.append(i)\n",
    "\n",
    "container"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets construct a more complicated loop to do what we want"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we do some things we don't need to do in the loop. Let's reload our atlas, and re-iniate our masker and correlation_measure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from nilearn.input_data import NiftiLabelsMasker\n",
    "from nilearn.connectome import ConnectivityMeasure\n",
    "from nilearn import datasets\n",
    "\n",
    "# load atlas\n",
    "multiscale = datasets.fetch_atlas_basc_multiscale_2015()\n",
    "atlas_filename = multiscale.scale064\n",
    "\n",
    "# initialize masker (change verbosity)\n",
    "masker = NiftiLabelsMasker(labels_img=atlas_filename, standardize=True, \n",
    "                           memory='nilearn_cache', verbose=0)\n",
    "\n",
    "# initialize correlation measure, set to vectorize\n",
    "correlation_measure = ConnectivityMeasure(kind='correlation', vectorize=True,\n",
    "                                         discard_diagonal=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "vectorized = ConnectivityMeasure(kind='correlation', vectorize=True, discard_diagonal=True)\n",
    "not_vectorized = ConnectivityMeasure(kind='correlation', vectorize=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# example\n",
    "sub = data.func[0]\n",
    "conf = prepare_confounds(data.confounds[0])\n",
    "time_series = masker.fit_transform(sub, confounds=conf)\n",
    "v_corr = vectorized.fit_transform([time_series])[0]\n",
    "nv_corr = not_vectorized.fit_transform([time_series])[0]\n",
    "\n",
    "print('vectorized:', v_corr.shape)\n",
    "print('not vectorized:', nv_corr.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# nfeat = \n",
    "# nfeat * (nfeat-1) / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay -- now that we have that cleared up, let's run our big loop!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NOTE**: On a laptop, this might take quite a long time. It took 30 minutes on my mac. \n",
    "\n",
    "**Maybe don't run it right now!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# all_features = [] # here is where we will put the data (a container)\n",
    "\n",
    "# for i,sub in enumerate(data.func):\n",
    "#     # convert confounds to readable format\n",
    "#     conf = prepare_confounds(data.confounds[i])\n",
    "#     # extract the timeseries from the ROIs in the atlas\n",
    "#     time_series = masker.fit_transform(sub, confounds=conf)\n",
    "#     # create a region x region correlation matrix\n",
    "#     correlation_matrix = correlation_measure.fit_transform([time_series])[0]\n",
    "#     # add to our container\n",
    "#     all_features.append(correlation_matrix)\n",
    "#     # keep track of status\n",
    "#     print('finished %s of %s'%(i+1,len(data.func)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Let's save the data to disk\n",
    "\n",
    "# np.savez_compressed('BASC064_features',a = all_features, b = full_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In case you do not want to run the full loop on your computer, you can load the output of the loop here!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "feat_file = 'BASC064_features.npz'\n",
    "X_features = np.load(feat_file)['a']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X_features.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay so we've got our features."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualize our feature matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.imshow(X_features, aspect='auto')\n",
    "plt.colorbar()\n",
    "plt.title('feature matrix')\n",
    "plt.xlabel('features')\n",
    "plt.ylabel('subjects')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare data for machine learning\n",
    "\n",
    "Here, we will define a \"training sample\" where we can play around with our models. We will also set aside a \"test\" sample that we will not touch until the end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to be sure that our training and test sample are matched! We can do that with a \"stratified split\". Specifically, we will stratify by age group."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "age_groups = pheno['AgeGroup']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# Split the sample to training/test with a 60/40 ratio, and \n",
    "# stratify by age group, and also shuffle the data.\n",
    "\n",
    "X_train, X_test, y_train, y_test, ageGroup_train, ageGroup_test = train_test_split(\n",
    "                                                                X_features, \n",
    "                                                                y_age, \n",
    "                                                                age_groups,\n",
    "                                                                test_size = 0.4, \n",
    "                                                                shuffle = True,\n",
    "                                                                stratify = age_groups,\n",
    "                                                                random_state = 123\n",
    "                                                                                   )\n",
    "\n",
    "# print the size of our training and test groups\n",
    "print('training:', len(X_train),\n",
    "     'testing:', len(X_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the distributions to be sure they are matched"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.hist(y_train, label = 'train')\n",
    "plt.hist(y_test, label = 'test')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run your first model!\n",
    "\n",
    "Machine learning can get pretty fancy pretty quickly. We'll start with a very standard linear model called a Support Vector Regressor (SVR). \n",
    "\n",
    "While this may seem unambitious, simple models can be very robust. And we don't have enough data to create more complex models.\n",
    "\n",
    "For more information, see this excellent resource:\n",
    "https://hal.inria.fr/hal-01824205"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's fit our first model!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.svm import SVR\n",
    "\n",
    "l_svr = SVR(kernel='linear') # define the model\n",
    "\n",
    "l_svr.fit(X_train, y_train) # fit the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well... that was easy. Let's see how well the model learned the data!\n",
    "\n",
    "We will judge our model on two criteria:\n",
    "* R2 = r-sqaure: the variance of the test data explained by the model\n",
    "* mae = mean absolute error: how off our measurements are in absolute units (years!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_absolute_error\n",
    "\n",
    "y_pred = l_svr.predict(X_train) # predict the training data based on the model\n",
    "r2 = l_svr.score(X_train, y_train) # get the r2\n",
    "mae = mean_absolute_error(y_true = y_train, \n",
    "                          y_pred = y_pred) # get the mae"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's view our results and plot them all at once!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print('r2 = %s, mae = %s'%(r2,mae))\n",
    "\n",
    "plt.scatter(y_train, y_pred)\n",
    "plt.title('Predicted vs Observed')\n",
    "plt.xlabel('Predicted Age')\n",
    "plt.ylabel('True age')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HOLY COW! Machine learning is amazing!!! Almost a perfect fit!\n",
    "\n",
    "...which means there's something wrong. What's the problem here?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_val_predict, cross_val_score\n",
    "\n",
    "# predict\n",
    "y_pred = cross_val_predict(l_svr, X_train, y_train, \n",
    "                           groups=ageGroup_train, cv=10)\n",
    "# scores\n",
    "r2 = cross_val_score(l_svr, X_train, y_train, \n",
    "                     groups=ageGroup_train, cv=10)\n",
    "mae_score = cross_val_score(l_svr, X_train, y_train, \n",
    "                            groups=ageGroup_train, cv=10,\n",
    "                           scoring = 'neg_mean_absolute_error')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can look at the accuracy of the predictions for each fold of the cross-validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "print(r2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also look at the overall accuracy of the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import r2_score\n",
    "\n",
    "overall_r2 = r2_score(y_pred = y_pred, y_true = y_train)\n",
    "overall_mae = mean_absolute_error(y_pred = y_pred, y_true = y_train)\n",
    "print('r2 = %s, mae = %s'%(overall_r2,overall_mae))\n",
    "\n",
    "plt.scatter(y_pred, y_train)\n",
    "plt.title('Predicted vs Observed')\n",
    "plt.xlabel('Predicted Age')\n",
    "plt.ylabel('True age')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not too bad at all! But what are some things you notice about our results?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tweak your model\n",
    "\n",
    "It's very important to learn when and where its appropriate to \"tweak\" your model.\n",
    "\n",
    "Since we have done all of the previous analysis in our training data, it's find to try different models. But we **absolutely cannot** \"test\" it on our left out data. If we do, we are in great danger of overfitting.\n",
    "\n",
    "We could try other models, or tweak hyperparameters, but we are probably not powered sufficiently to do so, and would once again risk overfitting.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But we can perhaps look at the performance of the model on log-transformed data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "log_y_train = np.log(y_train) # log-transform target data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# repeat the steps above to re-fit the model \n",
    "# and assess its performance\n",
    "\n",
    "# don't forget to switch y_train to log_y_train\n",
    "y_pred = cross_val_predict(l_svr, X_train, log_y_train, groups=ageGroup_train, cv=10)\n",
    "r2 = cross_val_score(l_svr, X_train, log_y_train, groups=ageGroup_train, cv=10)\n",
    "mae_score = cross_val_score(l_svr, X_train, log_y_train, groups=ageGroup_train, cv=10,\n",
    "                           scoring = 'neg_mean_absolute_error')\n",
    "\n",
    "# don't forget to switch y_train to log_y_train\n",
    "overall_r2 = r2_score(y_pred = y_pred, y_true = log_y_train)\n",
    "overall_mae = mean_absolute_error(y_pred = y_pred, y_true = log_y_train)\n",
    "print('r2 = %s, mae = %s'%(overall_r2,overall_mae))\n",
    "\n",
    "plt.scatter(y_pred, log_y_train)\n",
    "plt.title('Predicted vs Observed')\n",
    "plt.xlabel('Predicted Age')\n",
    "plt.ylabel('True age')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What do you think about the results of this model compared to the non-transformed model?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Try fitting a new SVR model and tweak one of the many parameters. Run cross-validation and see how well it goes. Make a new cell and type SVR? to see the possible hyperparameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# new_model = SVR() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Can our model predict age in completely un-seen data?\n",
    "Now that we've fit a model we think has possibly learned how to decode age based on rs-fmri signal, let's put it to the test. We will train our model on all of the training data, and try to predict the age of the subjects we left out at the beginning of this section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we performed a transformation on our training data, we will need to transform our testing data using the *same information!* \n",
    "\n",
    "For that, we will need to create a transformer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# transform testing set based on training distribution...\n",
    "from sklearn.preprocessing import FunctionTransformer\n",
    "\n",
    "transformer = FunctionTransformer(np.log).fit(y_train.values.reshape(-1,1))\n",
    "log_y_test = transformer.transform(y_test.values.reshape(-1,1))[:,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# how did we do?\n",
    "plt.hist(log_y_train, label = 'train')\n",
    "plt.hist(log_y_test, label = 'test')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now for the moment of truth! \n",
    "\n",
    "No cross-validation needed here. We simply fit the model with the training data and use it to predict the testing data\n",
    "\n",
    "I'm so nervous. Let's just do it all in one cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "l_svr.fit(X_train, log_y_train) # fit to training data\n",
    "y_pred = l_svr.predict(X_test) # predict age using testing data\n",
    "r2 = l_svr.score(X_test, log_y_test) # get r2 score\n",
    "mae = mean_absolute_error(y_pred=y_pred, y_true=log_y_test) # get mae\n",
    "\n",
    "# print results\n",
    "print('r2 = %s, mae = %s'%(r2,mae))\n",
    "\n",
    "# plot results\n",
    "plt.scatter(y_pred, log_y_test)\n",
    "plt.title('Predicted vs actual age')\n",
    "plt.xlabel('Predicted age')\n",
    "plt.ylabel('True age')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Wow!!*** Congratulations. You just trained a machine learning model that used real rs-fmri data to predict the age of real humans.\n",
    "\n",
    "It seems like something in this data does seem to be systematically related to age ... but what?  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpreting model feature importances\n",
    "Interpreting the feature importances of a machine learning model is a real can of worms. This is an area of active research. Unfortunately, it's hard to trust the feature importance of some models. \n",
    "\n",
    "You can find a whole tutorial on this subject here:\n",
    "http://gael-varoquaux.info/interpreting_ml_tuto/index.html\n",
    "\n",
    "For now, we'll just eschew better judgement and take a look at our feature importances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can access the feature importances (weights) used my the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "l_svr.coef_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "lets plot these weights to see their distribution better"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.bar(range(l_svr.coef_.shape[-1]),l_svr.coef_[0])\n",
    "plt.title('feature importances')\n",
    "plt.xlabel('feature')\n",
    "plt.ylabel('weight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or perhaps it will be easier to visualize this information as a matrix similar to the one we started with\n",
    "\n",
    "We can use the correlation measure from before to perform an inverse transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "correlation_measure.inverse_transform(l_svr.coef_).shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from nilearn import plotting\n",
    "\n",
    "feat_exp_matrix = correlation_measure.inverse_transform(l_svr.coef_)[0]\n",
    "\n",
    "plotting.plot_matrix(feat_exp_matrix, figure=(10, 8), \n",
    "                     labels=range(feat_exp_matrix.shape[0]),\n",
    "                     reorder=False,\n",
    "                    tri='lower')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see if we can throw those features onto an actual brain.\n",
    "\n",
    "First, we'll need to gather the coordinates of each ROI of our atlas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# from nibabel import Nifti1Image\n",
    "\n",
    "# def lazy_get_coords(atlas_filename):\n",
    "#     coords = []\n",
    "#     img = image.load_img(atlas_filename)\n",
    "#     atlas_data = img.get_data()\n",
    "#     aff = img.affine\n",
    "#     values = np.unique(atlas_data)[1:]\n",
    "#     for i in values:\n",
    "#         roi = np.zeros_like(atlas_data)\n",
    "#         roi[atlas_data==i] = 1.0\n",
    "#         coords.append(plotting.find_xyz_cut_coords(Nifti1Image(roi,aff)))\n",
    "    \n",
    "#     return coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# atlas_coordinates = lazy_get_coords(atlas_filename)\n",
    "# np.savez_compressed(os.path.join(outdir,'BASC064_coordinates'),\n",
    "#                    a = atlas_coordinates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "coords = np.load('BASC064_coordinates.npz')['a']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we can use our feature matrix and the wonders of nilearn to create a connectome map where each node is an ROI, and each connection is weighted by the importance of the feature to the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plotting.plot_connectome(feat_exp_matrix, coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Appendix: Goodies and extra points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** There are is a clear distinction between adults in children in this dataset. See if you can train a classifier (for example, SVC) to predict which subjects are adults and which are children. \n",
    "\n",
    "Return your overall accuracy, but also the positive predictive value (precisions score). You can tweak your model, but remember, don't burn your test data or your results don't count! Also remember you don't have many subjects to work with in the first place. \n",
    "\n",
    "For a bonus, plot the feature importances, and correlate the feature importances with those from our linear model predicting age"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# from sklearn.svm import SVC\n",
    "# y_class_train = [1 if x==b'Adult' else 0 for x in ageGroup_train]\n",
    "# y_class_test = [1 if x==b'Adult' else 0 for x in ageGroup_test]"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}