01.3-decoding_haxby.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> credits: originally by Gael Varoquaux\n",
    "> \n",
    "> shamelessly ~~stolen~~ adapted by: Chris Holdgraf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import pylab as plt\n",
    "import os.path as op\n",
    "path_data = op.join(op.expanduser('~'), 'nilearn_data/')\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# About this tutorial\n",
    "This is a quick introduction to neuroimaging decoding with `nilearn`. It covers the following topics:\n",
    "\n",
    "* Loading a neuroimaging dataset suitable for decoding\n",
    "* Extracting an ROI and vectorizing our 3D data so it may be passed to a `sklearn` object\n",
    "* Using pandas to munge some data\n",
    "* Fitting a `sklearn` classifier on our neuroimaging data\n",
    "* Visualizing classifier weights on the subject's brain\n",
    "\n",
    "  > Note: A lot of the material in this tutorial is drawn from the `nilearn` collection of examples. Open-source packages can be a great way to learn both about a package, and about the things that package tries to do (e.g., machine learning).\n",
    "\n",
    "  > In addition, many `nilearn` developers have recently released a paper covering the topics of decoding brain states with fMRI (and other modalities).\n",
    "\n",
    "  > * [Link to original `nilearn` tutorial](https://nilearn.github.io/auto_examples/plot_decoding_tutorial.html#sphx-glr-auto-examples-plot-decoding-tutorial-py)\n",
    "  > * [Link to Varoquaux decoding paper](https://arxiv.org/pdf/1606.05201v2.pdf)\n",
    "\n",
    "# A introduction tutorial to fMRI decoding (brain -> world)\n",
    "Thus far we've covered the general topics of machine learning, but how do they apply specifically to neuroscience data? \n",
    "\n",
    "This is a short tutorial on decoding with nilearn. It reproduces the\n",
    "Haxby 2001 study on a face vs cat discrimination task in a mask of the\n",
    "ventral stream.\n",
    "\n",
    "* [Here's a link](http://www.pymvpa.org/datadb/haxby2001.html) to the Haxby 2001 Dataset (w/ a link to the paper too)\n",
    "\n",
    "In a decoding model, we ask \"what information about the world can we predict using brain activity?\" These are often called \"backward models\" as they run counter to the flow of time (for most experimental setups).\n",
    "\n",
    "Examples of **inputs** to a decoding model include:\n",
    "* Mean voxel activity in each trial\n",
    "* Full timecourses of voxel activity (e.g., with continuously-varying stimuli)\n",
    "\n",
    "Examples of **outputs** from a decoding model include:\n",
    "\n",
    "* Experimental conditions\n",
    "* Muscle movement\n",
    "* Stimulus features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What's the difference between encoding and decoding?\n",
    "\n",
    "That is...a complicated question. The simplest answer is:\n",
    "\n",
    "* Encoding == forward model == world -> model -> brain\n",
    "* Decoding == backward model == brain -> model -> world\n",
    "\n",
    "People disagree about when you should use one vs the other, and what you can interpret from each. This lecture won't try to answer any of those questoins (we'd need several beers for this). However, here is a list of useful papers covering this topic:\n",
    "\n",
    "* [Weichwald et al, 2015](https://arxiv.org/pdf/1511.04780.pdf)\n",
    "* [Kay, 2017](http://www.sciencedirect.com/science/article/pii/S1053811917306638?via%3Dihub)\n",
    "* [Naselaris, 2012](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3037423/)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparing the data\n",
    "\n",
    "## Retrieve and load the fMRI data from the Haxby study\n",
    "\n",
    "The `nilearn.datasets.fetch_haxby` function will download the\n",
    "Haxby dataset if it's not present on the disk. It'll put this in the `nilearn` data directory and only needs to be downloaded once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import datasets\n",
    "from nilearn import plotting\n",
    "\n",
    "# By default the 2nd subject will be fetched\n",
    "haxby_dataset = datasets.fetch_haxby(data_dir=path_data)\n",
    "fmri_filename = haxby_dataset.func[0]\n",
    "\n",
    "# print basic information on the dataset\n",
    "print('First subject functional nifti images (4D) are at: %s' %\n",
    "      fmri_filename)  # 4D data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_filename = haxby_dataset.mask_vt[0]\n",
    "plotting.plot_roi(mask_filename, bg_img=haxby_dataset.anat[0],\n",
    "                  cmap='Paired')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.input_data import NiftiMasker\n",
    "# Load the mask from disk\n",
    "masker = NiftiMasker(mask_img=mask_filename, standardize=True)\n",
    "\n",
    "# Fitting the transformer initializes it to operate on new data\n",
    "masker.fit(fmri_filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Now we'll transform our fMRI data\n",
    "fmri_masked = masker.transform(fmri_filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The variable \"fmri_masked\" is a numpy array. It is 2-D."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(fmri_masked)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Its shape corresponds to the number of time-points x the number of\n",
    "voxels in the mask. Note that this is much fewer than the total number of voxels in the nifty image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(fmri_masked.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the behavioral labels\n",
    "\n",
    "Now that we have our vectorized fMRI activity, we need labels for the state of the experiment in order to fit our classifier. The behavioral labels are stored in a CSV file, separated by spaces.\n",
    "\n",
    "We use `pandas` to load them in an array. This is a library that is excellent for representing and manipulating tabular data. It's got a steep learning curve but is very powerful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "# Load target information as string and give a numerical identifier to each\n",
    "labels = pd.read_csv(haxby_dataset.session_target[0], delimiter=' ')\n",
    "print(labels.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like labels has the same length as our fMRI data, meaning that they share the same time-base."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(labels.shape)\n",
    "print(fmri_masked.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll retrieve the behavioral targets from the labels. These will be the \"classes\" that we attempt to predict.\n",
    "\n",
    "Note that these labels aren't integers like before. That's fine - `sklearn` will try to be clever and convert these into integer representations when we fit the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(labels['labels'].values[:50])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Restrict the analysis to cats and faces\n",
    "\n",
    "As we can see from the targets above, the experiment contains many\n",
    "conditions. Today we'll restrict the decoding to two categories of interest: cats and faces.\n",
    "\n",
    "To do this, we'll use `pandas` to create a mask corresponding to these categories, and then use it to extract only the rows we care about."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create a mask w/ Pandas\n",
    "condition_mask = labels.eval('labels in [\"face\", \"cat\"]').values\n",
    "\n",
    "# Create a mask w/ Numpy\n",
    "# condition_mask = np.logical_or(target == b'face', target == b'cat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We apply this mask along the \"samples\" axis to restrict the\n",
    "# classification to the face vs cat discrimination\n",
    "fmri_masked = fmri_masked[condition_mask]\n",
    "targets = labels[condition_mask]['labels'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we now have fewer samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(fmri_masked.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Decoding with an SVM\n",
    "\n",
    "Now we have all the components we'll need to fit a model. We have:\n",
    "\n",
    "* Masked a subset of voxels in which we are interested.\n",
    "* Vectorized those masked voxels so that they have shape (n_samples, n_features). **This is `X`**\n",
    "* Extracted a set of labels, one for each timepoint, corresponding to the stimulus being shown at that time. **This is `y`**\n",
    "* Masked our time dimension so that we only have two classes of interest.\n",
    "\n",
    "Now we'll fit our model. As before, we'll use the [scikit-learn](http://www.scikit-learn.org>)  toolbox on the fmri_masked data.\n",
    "\n",
    "As a decoder, we'll use a Support Vector Classification, with a linear\n",
    "kernel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.svm import SVC\n",
    "svc = SVC(kernel='linear')\n",
    "print(svc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As our data is already in the shape for `sklearn`, we can quickly fit the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "svc.fit(fmri_masked, targets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predicting with / Scoring our model\n",
    "Machine learning is all about **prediction**. As such, we'll use our fit model to make a prediction about the class of some data, given the structure the model has already found:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we'll predict with the same input training data\n",
    "prediction = svc.predict(fmri_masked)\n",
    "print(prediction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`sklearn` has a number of functions for defining the \"score\" of a model. The proper one to use depends on the nature of your model and data. Here we'll use a simple scorer for classification. These functions expect two arrays:\n",
    "\n",
    "* The array of \"true\" classes for each sample\n",
    "* The array of predicted classes given our model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "print(accuracy_score(targets, prediction))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wow, a score of 100%! Ship it off to *Nature*, right? Not quite yet. This is because you should **always make predictions on data that hasn't been used to fit the model**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Validating our model\n",
    "\n",
    "The proper way to measure error rates or prediction accuracy is via\n",
    "cross-validation: leaving out some data and testing on it. There are many ways to do this.\n",
    "\n",
    "Here we'll do this by **manually leaving out datapoints** during fit. We'll set them aside and use them to score the model.\n",
    "\n",
    "Let's leave out the 30 last data points during training, and test the\n",
    "prediction on these 30 last points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "svc.fit(fmri_masked[:-30], targets[:-30])\n",
    "\n",
    "prediction = svc.predict(fmri_masked[-30:])\n",
    "print((prediction == targets[-30:]).sum() / len(targets[-30:]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, this seems unfortunate. We've now got 50% less data in order to fit the model. Ideally, we'd like to do two things:\n",
    "\n",
    "* Validate our model properly (aka, on held-out data not used in model fitting)\n",
    "* Use as much of our data as possible.\n",
    "\n",
    "It is difficult to satisfy both of these conditions properly, but *cross-validation* is one way of getting closer to this goal. \n",
    "\n",
    "For more details on this, check out the [cross validation notebook](01.5-cross_validation.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspecting the model weights\n",
    "\n",
    "\n",
    "Finally, now that our model has been fit to data, and validated on held-out data, it may be useful to inspect and display the model weights. This is often used to understand the voxels that were particularly important in discriminating these two classes.\n",
    "\n",
    "`sklearn` models that are linear store their weights in an attribute called `coef_`. We'll look at this below. Note that there is one weight per feature (in this case, voxels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coef_ = svc.coef_\n",
    "print(coef_[:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(svc.coef_.shape)\n",
    "print(fmri_masked.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Turning the weights into a nifti image\n",
    "Using our `NiftiMasker`, we can collect these weights and reshape them back into the 3-D space of our fMRI data.\n",
    "\n",
    "We need to turn it back into a Nifti image, in essence, \"inverting\"\n",
    "what the NiftiMasker has done.\n",
    "\n",
    "For this, we can call inverse_transform on the NiftiMasker:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coef_img = masker.inverse_transform(coef_)\n",
    "print(coef_img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "coef_img is now a NiftiImage. Essentially, this is like the statistical \"t-maps\" that we've visualized before. \n",
    "\n",
    "We can save the coefficients as a nii.gz file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "coef_img.to_filename('haxby_svc_weights.nii.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing our weights\n",
    "\n",
    "Finally, we can plot the weights that the model found, using the subject's anatomical as a background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.plotting import plot_stat_map, show\n",
    "\n",
    "plot_stat_map(coef_img, bg_img=haxby_dataset.anat[0],\n",
    "              title=\"SVM weights\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recap\n",
    "Above we've shown a sample pipeline for fitting a linear classifier using fMRI activity and object classes. We've covered a few of the basics in the machine learning pipeline, and shown how `nilearn` and `sklearn` complement one another.\n",
    "\n",
    "Much more information about these functions, and the options available for machine learning, can be found in the `nilearn` documentation."
   ]
  }
 ],
 "metadata": {
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}