01.2-fmri_to_sklearn.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nilearn: From fMRI to scikit-learn\n",
    "\n",
    "## Vectorizing is really useful...\n",
    "\n",
    "Scikit-learn has a recurring pattern in its API. Whether you're fitting a regressor, a classifier, or a feature transformation object, they almost all take data as an `X` and (sometimes) `y` array.\n",
    "\n",
    "However, in neuroimaging (and many other fields) our data doesn't naturally come in this form.\n",
    "\n",
    "* sklearn shapes: (samples x features)\n",
    "* nipy shapes: (x, y, z, time/conditions/etc)\n",
    "\n",
    "The process of converting your Nifti-shaped files into a scikit-learn format is often called **vectorizing**. Effectively, this means going from:\n",
    "\n",
    "**your data shape -> (samples, features)**\n",
    "\n",
    "This notebook shows you how to do this in the context of `nilearn`."
   ]
  },
  {
   "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",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn import datasets\n",
    "from nilearn import plotting\n",
    "from nilearn import image\n",
    "\n",
    "# By default the 2nd subject will be fetched\n",
    "haxby_dataset = datasets.fetch_haxby()\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": "markdown",
   "metadata": {},
   "source": [
    "## Preparing the fMRI volume for model fitting\n",
    "Machine learning generally follows this pattern:\n",
    "\n",
    "1. Data preprocessing\n",
    "2. feature extraction/engineering\n",
    "3. model fitting\n",
    "4. model validation\n",
    "5. model inspection\n",
    "\n",
    "This covers a part of point 2.\n",
    "\n",
    "### Extracting a subset of voxels\n",
    "We can easily load volume data into python using `nibabel`. However, we generally don't want the full volume, but instead prefer to use a subset of voxels based on some preferred region of the brain (e.g. the cortical surface).\n",
    "\n",
    "Effectively, we wish to do two things:\n",
    "\n",
    "1. Extract a subset of voxels from the nifty file (masking)\n",
    "2. Reshape these voxels so that they can be used to fit a model (vectorizing)\n",
    "\n",
    "`nilearn` allows us to do this easily, turning the 4-D matrix of into an array of shape (n_samples, n_features). This allows us to use scikit-learn to do machine learning.\n",
    "\n",
    "### Visualizing the mask\n",
    "First, remember that masks are simply a 3D array in the same space as our MRI data. They are basically boolean values that say whether or not to keep each voxel. \n",
    "\n",
    "We'll load and plot it below. The mask is a mask of the Ventral Temporal streaming area coming from the Haxby study:"
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "Notice that only a few of the voxels are being shown here. These are the voxels that the mask has extracted. Now we use the NiftiMasker.\n",
    "\n",
    "### Masking / vectorizing our data\n",
    "We will use the `nilearn.input_data.NiftiMasker` to extract the\n",
    "fMRI data on a mask and convert it to data series.\n",
    "\n",
    "We first create a masker, giving it the options that we care\n",
    "about. Here we also use standardizing of the data, as it is often important\n",
    "for decoding. This scales the data so that its mean / variance is more consistent across voxels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.input_data import NiftiMasker\n",
    "masker = NiftiMasker(mask_img=mask_filename, standardize=True)\n",
    "masker"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**A note on transformers**\n",
    "\n",
    "The `masker` is an object that can perform data masking for us. It follows a very similar API to `scikit-learn`. Objects that transform other objects are called `transformers`, and generally have a `transform` method on top of the `fit` method that we covered earlier.\n",
    "\n",
    "We'll `fit` the masker on our mask, then use it to vectorize our data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fitting the transformer initializes it to operate on new data\n",
    "masker.fit(fmri_filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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": [
    "### De-vectorizing back to a nifti image\n",
    "After we've done our analysis using the vectorized data, we then want to reshape it back to `i,j,k` space. \n",
    "\n",
    "Using our `NiftiMasker`, we can turn the dta 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": [
    "nifti_img = masker.inverse_transform(fmri_masked)\n",
    "print(nifti_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting.plot_stat_map(image.index_img(nifti_img, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that maskers will work for different kinds of inputs to `inverse_transform`..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that this will work fine\n",
    "nifti_img = masker.inverse_transform(fmri_masked[:10])\n",
    "print(nifti_img.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Or even this...\n",
    "nifti_img = masker.inverse_transform(fmri_masked[0])\n",
    "print(nifti_img.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# But changing the number of features will break it...\n",
    "nifti_img = masker.inverse_transform(fmri_masked[:, :-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vectorizing (by itself) removes information\n",
    "\n",
    "Vectorizing is often the final step of the processing chain before we actually fit a model on our data. It often leaves the data in a state with _minimal_ knowledge about the relationships between features. Make sure that you keep track of what each feature means (or use a tool like the Nifti masker) to know how to interpret your model results. "
   ]
  }
 ],
 "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
}