01.4-encoding_miyawaki.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",
    "\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "from nilearn.input_data import MultiNiftiMasker\n",
    "from nilearn.datasets import fetch_miyawaki2008\n",
    "import pylab as plt\n",
    "import os.path as op\n",
    "from tqdm import tqdm\n",
    "\n",
    "sns.set_style('white')\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# About this tutorial\n",
    "This example is also drawn from the `nilearn` docs. You can find the original [here](https://nilearn.github.io/auto_examples/02_decoding/plot_miyawaki_encoding.html#sphx-glr-auto-examples-02-decoding-plot-miyawaki-encoding-py).\n",
    "\n",
    "This example partly reproduces the encoding model presented in\n",
    "> [Visual image reconstruction from human brain activity\n",
    "    using a combination of multiscale local image decoders](http://www.cell.com/neuron/abstract/S0896-6273%2808%2900958-6),\n",
    "    Miyawaki, Y., Uchida, H., Yamashita, O., Sato, M. A.,\n",
    "    Morito, Y., Tanabe, H. C., ... & Kamitani, Y. (2008).\n",
    "    Neuron, 60(5), 915-929.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Encoding models for visual stimuli\n",
    "\n",
    "Thus far we have focused on using predictive models for decoding. That is - predicting states of the world using brain activity. However, it is possible to do the reverse: predict brain activity using information about the world. These are known as *encoding* models.\n",
    "It is possible to fit encoding models using `nilearn` and `scikit-learn`. Here we'll use an encoding model to predict **fMRI data** from **visual stimuli**, using the dataset from\n",
    "[Miyawaki et al., 2008](http://www.cell.com/neuron/abstract/S0896-6273%2808%2900958-6>).\n",
    "\n",
    "## The Miyawaki dataset\n",
    "Participants were shown images, which consisted of random 10x10 binary\n",
    "(either black or white) pixels, and the corresponding fMRI activity was\n",
    "recorded. We will try to predict the activity in each voxel\n",
    "from the binary pixel-values of the presented images. Then we extract the\n",
    "receptive fields for a set of voxels to see which pixel location a voxel\n",
    "is most sensitive to."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preparing the data\n",
    "\n",
    "First we'll load the full dataset using `nilearn`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = fetch_miyawaki2008()\n",
    "print(dataset['description'].decode('utf-8'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This dataset is broken up into both training and testing data. Here we'll only use the training data of this study, where random binary images were shown."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# training data starts after the first 12 files\n",
    "fmri_random_runs_filenames = dataset.func[12:]\n",
    "stimuli_random_runs_filenames = dataset.label[12:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the neural data\n",
    "### Working with multiple data files\n",
    "\n",
    "The `nilearn.input_data.MultiNiftiMasker` class allows us to mask multiple Nifti files at once. We'll use this to load the fMRI data, clean and mask it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is a list of filenames, corresponding to multiple runs of the experiment\n",
    "fmri_random_runs_filenames[:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "masker = MultiNiftiMasker(mask_img=dataset.mask, detrend=True,\n",
    "                          standardize=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vectorizing our data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# This will mask and clean each dataset\n",
    "# It outputs a list of vectorized arrays\n",
    "masker.fit()\n",
    "fmri_data = masker.transform(fmri_random_runs_filenames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Type: ', type(fmri_data), '\\n---\\nLength: ', len(fmri_data))\n",
    "fmri_data[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading in stimuli\n",
    "\n",
    "These files define the stimuli that are used in this dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stimuli_random_runs_filenames[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since it's a CSV, we'll again use pandas to load it. Each CSV is one run. Within each CSV, each row is one stimulus presentation, consisting of a 10x10 grid of black and white squares:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "stimulus = pd.read_csv(stimuli_random_runs_filenames[0])\n",
    "stimulus.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reshaping so we can plot them\n",
    "\n",
    "Here we'll reshape our stimulus data so we can get a feel for what they look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# shape of the binary (i.e. black and wihte values) image in pixels\n",
    "stimulus_shape = (10, 10)\n",
    "\n",
    "# We load the visual stimuli from csv files\n",
    "stimuli = []\n",
    "for stimulus_run in stimuli_random_runs_filenames:\n",
    "    stimulus = pd.read_csv(stimulus_run, header=None).values.astype(int)\n",
    "    \n",
    "    # Reshape so it's 10 x 10\n",
    "    stimulus = stimulus.reshape((-1,) + stimulus_shape, order='F')\n",
    "    stimuli.append(stimulus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at some of these binary images:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 2, figsize=(10, 5))\n",
    "for (run, stimulus), ax in zip([(0, 124), (2, 101)], axs):\n",
    "    ax.imshow(stimuli[run][stimulus], interpolation='nearest', cmap='gray')\n",
    "    ax.set_title('Run {}, Stimulus {}'.format(run, stimulus))\n",
    "    ax.set_axis_off()\n",
    "plt.tight_layout(w_pad=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reshaping data for the model\n",
    "We now stack the fmri and stimulus data and remove an offset in the\n",
    "beginning/end. This will prep the data for our model fitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Removing the offset (2 values) so they're aligned properly\n",
    "fmri_data = np.vstack([fmri_run[2:] for fmri_run in fmri_data])\n",
    "stimuli = np.vstack([stimuli_run[:-2] for stimuli_run in stimuli]).astype(float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "fmri_data is a matrix of *samples* x *voxels*\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(fmri_data.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We flatten the last two dimensions of stimuli\n",
    "so it is a matrix of *samples* x *pixels*.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Flatten the stimuli\n",
    "stimuli = np.reshape(stimuli, (-1, stimulus_shape[0] * stimulus_shape[1]))\n",
    "\n",
    "print(stimuli.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: We don't have a fancy `ImageVectorizer` object so we're doing the reshaping by hand. If you were implementing this code in your own work, it may be a good idea to write your own `Vectorizer` object."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building the encoding models\n",
    "\n",
    "In the last tutorial we built a classifier using brain activity and labels for different stimuli that were presented. In this case, we wish to make a different kind of prediction: a *continuous* representation of the each voxel using the pixel data from the stimuli.\n",
    "\n",
    "These types of models are called **voxel-wise encoding models**, and they use *regression* to predict continuous outputs from data. Specifically, we'll use [Ridge regression](http://en.wikipedia.org/wiki/Tikhonov_regularization).\n",
    "For each voxel we fit an independent regression model,\n",
    "using the pixel-values of the visual stimuli to predict the neuronal\n",
    "activity in this voxel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import Ridge\n",
    "from sklearn.model_selection import KFold\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-validating our encoding model\n",
    "Using 10-fold cross-validation, we partition the data into 10 'folds'.\n",
    "We hold out each fold of the data for testing, then fit a ridge regression\n",
    "to the remaining 9/10 of the data.\n",
    "\n",
    "As this is an encoding model, we'll use stimuli as predictors\n",
    "and fmri_data as targets, and create predictions for the held-out 10th set.\n",
    "\n",
    "In addition, we'll **score the model** using the coefficient of determination ($R^2$). This estimates the model's ability to explain variance in the neural activity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import r2_score\n",
    "\n",
    "# Define our estimator and cross-validation objects...\n",
    "estimator = Ridge(alpha=100.)\n",
    "cv = KFold(n_splits=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a quick look at what our cross-validation is actually doing..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "for ii, (tr, tt) in enumerate(cv.split(range(10))):\n",
    "    ax.plot([tr, tr], [ii, ii+1], 'C0', lw=10)\n",
    "    ax.plot([tt, tt], [ii, ii+1], 'C2', lw=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll loop through each iteration of the CV, fit the model on a subset of data, and test it on another subset of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scores = []\n",
    "for train, test in tqdm(list(cv.split(fmri_data))):\n",
    "    # we train the Ridge estimator on the training set\n",
    "    X = stimuli.reshape(-1, 100)\n",
    "    X_train = X[train]\n",
    "    X_test = X[test]\n",
    "    y_train = fmri_data[train]\n",
    "    y_test = fmri_data[test]\n",
    "\n",
    "    # and predict the fMRI activity for the test set\n",
    "    model = Ridge(alpha=100.).fit(X_train, y_train)\n",
    "    predictions = model.predict(X_test)\n",
    "\n",
    "    # we compute how much variance our encoding model explains in each voxel\n",
    "    scores.append(r2_score(y_test, predictions, multioutput='raw_values'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mapping the encoding scores on the brain\n",
    "\n",
    "We've fit an encoding model for each voxel of the brain. Now we can inspect the model's performance at each location.\n",
    "\n",
    "To plot the scores onto the brain, we create a Nifti1Image using\n",
    "the scores output by our model. We'll again use the `NiftiMasker` to convert it back to volume space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cut_score = np.mean(scores, axis=0)\n",
    "cut_score[cut_score < 0] = 0\n",
    "\n",
    "# bring the scores into the shape of the background brain\n",
    "score_map_img = masker.inverse_transform(cut_score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.plotting import plot_stat_map\n",
    "display = plot_stat_map(score_map_img, bg_img=dataset.background,\n",
    "                        cut_coords=[-8], display_mode='z', aspect=1.25,\n",
    "                        title='Explained variance per voxel')\n",
    "plt.gcf().set_size_inches(5, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What does it mean to have a \"good\" score in regression?\n",
    "\n",
    "In regression we're predicting a continuous variable. If we do a \"good\" job of predicting, then our predictions should be correlated with the \"true\" activity of the voxel. Let's visualize a good and bad voxel below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.regplot(y_test[:, 0], predictions[:, 0])\n",
    "ax.set(title=\"A bad voxel\", xlabel=\"Voxel activity (true)\", ylabel=\"Voxel activity (predicted)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "good_scores_ixs = np.where(cut_score > .3)[0]\n",
    "this_score_ix = good_scores_ixs[10]\n",
    "\n",
    "ax = sns.regplot(y_test[:, this_score_ix], predictions[:, this_score_ix])\n",
    "ax.set(title=\"A good voxel\", xlabel=\"Voxel activity (true)\", ylabel=\"Voxel activity (predicted)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thresholding our scores and making a final plot\n",
    "We may also want to threshold our scores so that only \"valid\" scores are plotted. We'll do this thresholding below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from nilearn.image import threshold_img\n",
    "thresholded_score_map_img = threshold_img(score_map_img, threshold=.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll highlight a few voxels of particular interest. We'll inspect these more closely later on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.image.resampling import coord_transform\n",
    "\n",
    "def index_to_xy_coord(x, y, z=10):\n",
    "    '''Transforms data index to coordinates of the background + offset'''\n",
    "    coords = coord_transform(x, y, z,\n",
    "                             affine=thresholded_score_map_img.affine)\n",
    "    return np.array(coords)[np.newaxis, :] + np.array([0, 1, 0])\n",
    "\n",
    "\n",
    "xy_indices_of_special_voxels = [(30, 10), (32, 10), (31, 9), (31, 10)]\n",
    "\n",
    "display = plot_stat_map(thresholded_score_map_img, bg_img=dataset.background,\n",
    "                        cut_coords=[-8], display_mode='z', aspect=1.25,\n",
    "                        title='Explained variance per voxel')\n",
    "\n",
    "\n",
    "# creating a marker for each voxel and adding it to the statistical map\n",
    "for i, (x, y) in enumerate(xy_indices_of_special_voxels):\n",
    "    display.add_markers(index_to_xy_coord(x, y), marker_color='none',\n",
    "                        edgecolor=['b', 'r', 'magenta', 'g'][i],\n",
    "                        marker_size=140, marker='s',\n",
    "                        facecolor='none', lw=4.5)\n",
    "\n",
    "\n",
    "# re-set figure size after construction so colorbar gets rescaled too\n",
    "fig = plt.gcf()\n",
    "fig.set_size_inches(12, 12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model introspection: Estimating receptive fields\n",
    "\n",
    "Now that we can see the distribution of scores across the voxels of interest, let's take a closer look at the receptive fields for the four marked voxels. *It only makes sense to do this for voxels with a relatively high prediction score*.\n",
    "\n",
    "A voxel's [receptive field](http://en.wikipedia.org/wiki/Receptive_field)\n",
    "is the region of a stimulus (like an image) where the presence of an object (like a white instead of a black pixel), results in a change in activity\n",
    "in the voxel.\n",
    "\n",
    "In our case the receptive field is just the vector of 100\n",
    "regression  coefficients (one for each pixel) reshaped into the 10x10\n",
    "form of the original images.\n",
    "\n",
    "> The kind of model we choose to fit will affect the kind of coefficients we will uncover. For example, *Ridge regression* will tend to find a set of coefficients that are relatively normally distributed around 0. [*Lasso regression*](http://en.wikipedia.org/wiki/Lasso_(statistics)) will find a \"sparse\" solution such that only a few coefficients will be non-zero. You can play around with the effects of each below.\n",
    "\n",
    "**In practice, you should choose the model that best fits your knowledge of the underlying relationship between brain / stimuli.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LassoLarsCV, RidgeCV\n",
    "from matplotlib import gridspec\n",
    "from matplotlib.patches import Rectangle\n",
    "\n",
    "# automatically estimate the sparsity by cross-validation\n",
    "lasso = LassoLarsCV(max_iter=10)\n",
    "# Or find a smooth set of coefficients\n",
    "ridge = RidgeCV(alphas=np.logspace(-2, 2, 10))\n",
    "\n",
    "# This is the model we'll actually use\n",
    "model = lasso  # ridge\n",
    "\n",
    "# These are the indices for voxels for which we'll estimate a receptive field.\n",
    "voxels_to_fit = [1780, 1951, 2131, 1935]\n",
    "\n",
    "# This is a pixel we'll highlight in each model's receptive field.\n",
    "marked_pixel = (4, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll loop through each voxel, fit a regression model for it, and visualize the coefficients of that model.\n",
    "\n",
    "We'll plot the voxels so that the plots represent the relative positions of voxels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(2, 3, figsize=(12, 8))\n",
    "axs_to_plt = [0, 1, 2, 4]\n",
    "\n",
    "# we fit the model for each of the three voxels of the upper row\n",
    "for i, index in enumerate(voxels_to_fit):\n",
    "    # --- Fit the model ---\n",
    "    y = fmri_data[:, index]\n",
    "    model.fit(stimuli, y)\n",
    "\n",
    "    # we reshape the coefficients into the form of the original images\n",
    "    rf = model.coef_.reshape((10, 10))\n",
    "\n",
    "    # --- Visualize the RF ---\n",
    "    # add a black background\n",
    "    ax = axs.ravel()[axs_to_plt[i]]\n",
    "    ax.imshow(np.zeros_like(rf), vmin=0., vmax=1., cmap='gray')\n",
    "    ax_im = ax.imshow(np.ma.masked_less(rf, 0.1), interpolation=\"nearest\",\n",
    "                      cmap=['Blues', 'Greens', 'Reds', 'Purples'][i], vmin=0., vmax=0.75)\n",
    "    # add the marked pixel\n",
    "    ax.add_patch(Rectangle(\n",
    "        (marked_pixel[1] - .5, marked_pixel[0] - .5), 1, 1,\n",
    "        facecolor='none', edgecolor='r', lw=4))\n",
    "    plt.colorbar(ax_im, ax=ax, shrink=.7)\n",
    "\n",
    "for ax in axs.ravel()[[3, 5]]:\n",
    "    ax.set_axis_off()\n",
    "fig.suptitle('Receptive fields of the marked voxels', fontsize=25)\n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that these plots are organized according to each voxel's position relative to one another. Given this we can note two things:\n",
    "\n",
    "1. The receptive fields of the four voxels are close to each other\n",
    "2. The relative location of the \"preferred\" pixel in each RF roughly follows the voxels' position relative to one another.\n",
    "\n",
    "In this way, we've mapped out how the receptive field changes as you move across voxels in the brain. Taken as a whole, we have an idea for how spatial / visual informatino is represented in the brain, and how this representation changes as we move throughout regions of the brain.\n",
    "\n",
    "# Recap\n",
    "In this section, we've covered how we can use `sklearn` and `nilearn` in order to use information in the world to model brain activity. This is another exciting approach towards asking questions about the brain, and it only scratches the surface of what you can do with these models.\n",
    "\n"
   ]
  }
 ],
 "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
}