{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Nistats: Functional MRI in Python\n",
    "=========================================================\n",
    "\n",
    "[Nistats](https://nistats.github.io/) is a Python module for fast and easy functional MRI statistical analysis.\n",
    "\n",
    "It leverages [Nilearn](), [Nibabel]() and other Python libraries from the Python scientific stack like [Scipy](), [Numpy]() and [Pandas]().\n",
    "\n",
    "In this tutorial, we're going to explore `nistats` functionality by analyzing a single subject single run example using a General Linear Model (GLM). We're gonna use the same example dataset (ds000114) as from the `nibabel`\n",
    "and `nilearn` tutorials. As this is a multi run multi task dataset, we've to decide on a run and a task we want to analyze. Let's go with `ses-test` and `task-fingerfootlips` of `sub-01`.\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setting and inspecting the data\n",
    "=========================\n",
    "\n",
    "At first, we have to set and indicate the data we want to analyze. As stated above, we're going to use the anatomical image and the preprocessed functional image of `sub-01` from `ses-test`. The preprocessing was conducted through [fmriprep](https://fmriprep.readthedocs.io/en/stable/index.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fmri_img = '/data/ds000114/derivatives/fmriprep/sub-01/ses-test/func/sub-01_ses-test_task-fingerfootlips_bold_space-mni152nlin2009casym_preproc.nii.gz'\n",
    "anat_img = '/data/ds000114/sub-01/ses-test/anat/sub-01_ses-test_T1w.nii.gz'\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can display the mean functional image and the subject's anatomy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.image import mean_img\n",
    "mean_img = mean_img(fmri_img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nilearn.plotting import plot_stat_map, plot_anat, plot_img, show\n",
    "plot_img(mean_img)\n",
    "plot_anat(anat_img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specifying the experimental paradigm\n",
    "------------------------------------\n",
    "\n",
    "We must now provide a description of the experiment, that is, define the\n",
    "timing of the task and rest periods. This is typically\n",
    "provided in an events.tsv file.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "events = pd.read_table('/data/ds000114/task-fingerfootlips_events.tsv')\n",
    "print(events)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Performing the GLM analysis\n",
    "---------------------------\n",
    "\n",
    "It is now time to create and estimate a ``FirstLevelModel`` object, that will generate the *design matrix* using the  information provided by the ``events`` object.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nistats.first_level_model import FirstLevelModel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameters of the first-level model\n",
    "\n",
    "* t_r=?? is the time of repetition of acquisitions\n",
    "* noise_model='ar1' specifies the noise covariance model: a lag-1 dependence\n",
    "* standardize=False means that we do not want to rescale the time\n",
    "series to mean 0, variance 1\n",
    "* hrf_model='spm' means that we rely on the SPM \"canonical hrf\"\n",
    "model (without time or dispersion derivatives)\n",
    "* drift_model='cosine' means that we model the signal drifts as slow oscillating time functions\n",
    "* period_cut=160(s) defines the cutoff frequency (its inverse actually).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need the TR of the functional images, luckily we can extract that information using `nibabel`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!nib-ls /data/ds000114/derivatives/fmriprep/sub-01/ses-test/func/sub-01_ses-test_task-fingerfootlips_bold_space-mni152nlin2009casym_preproc.nii.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see the `TR` is 2.5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fmri_glm = FirstLevelModel(t_r=2.5,\n",
    "                           noise_model='ar1',\n",
    "                           standardize=False,\n",
    "                           hrf_model='spm',\n",
    "                           drift_model='cosine',\n",
    "                           period_cut=160)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usually, we also want to include confounds computed during preprocessing (e.g., motion, global signal, etc.) as regressors of no interest. In our example, these were computed by `fmriprep` and can be found in `derivatives/fmriprep/sub-01/func/`. We can use `pandas` to inspect that file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "confounds = pd.read_csv('/data/ds000114/derivatives/fmriprep/sub-01/ses-test/func/sub-01_ses-test_task-fingerfootlips_bold_confounds.tsv', delimiter='\\t')\n",
    "confounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparable to other neuroimaging softwards, we have a timepoint x confound dataframe. However, `fmriprep` computes way more confounds than most of you are used to and that require a bit of reading to understand and therefore utilize properly. We therefore and for the sake of simplicity stick to the \"classic\" ones: `WhiteMatter`, `GlobalSignal`, `FramewiseDisplacement` and the `motion correction parameters` in translation and rotation: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "confounds_glm = confounds[['WhiteMatter', 'GlobalSignal', 'FramewiseDisplacement', 'X', 'Y', 'Z', 'RotX', 'RotY', 'RotZ']].replace(np.nan, 0)\n",
    "confounds_glm\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have specified the model, we can run it on the fMRI image\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fmri_glm = fmri_glm.fit(fmri_img, events, confounds_glm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can inspect the design matrix (rows represent time, and\n",
    "columns contain the predictors).\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "design_matrix = fmri_glm.design_matrices_[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Formally, we have taken the first design matrix, because the model is\n",
    "implictily meant to for multiple runs.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nistats.reporting import plot_design_matrix\n",
    "plot_design_matrix(design_matrix)\n",
    "import matplotlib.pyplot as plt\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the design matrix image to disk, first creating a directory where you want to write the images:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "outdir = 'results'\n",
    "if not os.path.exists(outdir):\n",
    "    os.mkdir(outdir)\n",
    "\n",
    "from os.path import join\n",
    "plot_design_matrix(design_matrix, output_file=join(outdir, 'design_matrix.png'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first column contains the expected reponse profile of regions which are\n",
    "sensitive to the \"Finger\" task. Let's plot this first column:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(design_matrix['Finger'])\n",
    "plt.xlabel('scan')\n",
    "plt.title('Expected Response for condition \"Finger\"')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Detecting voxels with significant effects\n",
    "-----------------------------------------\n",
    "\n",
    "To access the estimated coefficients (Betas of the GLM model), we\n",
    "created constrast with a single '1' in each of the columns: The role\n",
    "of the contrast is to select some columns of the model --and\n",
    "potentially weight them-- to study the associated statistics. So in\n",
    "a nutshell, a contrast is a weigted combination of the estimated\n",
    "effects.  Here we can define canonical contrasts that just consider\n",
    "the two condition in isolation ---let's call them \"conditions\"---\n",
    "then a contrast that makes the difference between these conditions.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import array\n",
    "conditions = {\n",
    "    'active - Finger': array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),\n",
    "    'active - Foot':   array([0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),\n",
    "    'active - Lips':   array([0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at it: plot the coefficients of the contrast, indexed by\n",
    "the names of the columns of the design matrix.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nistats.reporting import plot_contrast_matrix\n",
    "plot_contrast_matrix(conditions['active - Finger'], design_matrix=design_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we compute the estimated effect. It is in BOLD signal unit,\n",
    "but has no statistical guarantees, because it does not take into\n",
    "account the associated variance.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eff_map = fmri_glm.compute_contrast(conditions['active - Finger'],\n",
    "                                    output_type='effect_size')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to get statistical significance, we form a t-statistic, and\n",
    "directly convert is into z-scale. The z-scale means that the values\n",
    "are scaled to match a standard Gaussian distribution (mean=0,\n",
    "variance=1), across voxels, if there were now effects in the data.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_map = fmri_glm.compute_contrast(conditions['active - Finger'],\n",
    "                                  output_type='z_score')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot thresholded z scores map.\n",
    "\n",
    "We display it on top of the average\n",
    "functional image of the series (could be the anatomical image of the\n",
    "subject).  We use arbitrarily a threshold of 3.0 in z-scale. We'll\n",
    "see later how to use corrected thresholds.  we show to display 3\n",
    "axial views: display_mode='z', cut_coords=3\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_stat_map(z_map, bg_img=mean_img, threshold=3.0,\n",
    "              display_mode='z', cut_coords=3, black_bg=True,\n",
    "              title='active - Finger (Z>3)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Statistical signifiance testing. One should worry about the\n",
    "statistical validity of the procedure: here we used an arbitrary\n",
    "threshold of 3.0 but the threshold should provide some guarantees on\n",
    "the risk of false detections (aka type-1 errors in statistics). One\n",
    "first suggestion is to control the false positive rate (fpr) at a\n",
    "certain level, e.g. 0.001: this means that there is.1% chance of\n",
    "declaring active an inactive voxel.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nistats.thresholding import map_threshold\n",
    "_, threshold = map_threshold(z_map, level=.001, height_control='fpr')\n",
    "print('Uncorrected p<0.001 threshold: %.3f' % threshold)\n",
    "plot_stat_map(z_map, bg_img=mean_img, threshold=threshold,\n",
    "              display_mode='z', cut_coords=3, black_bg=True,\n",
    "              title='active - Finger (p<0.001)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem is that with this you expect 0.001 * n_voxels to show up\n",
    "while they're not active --- tens to hundreds of voxels. A more\n",
    "conservative solution is to control the family wise errro rate,\n",
    "i.e. the probability of making ony one false detection, say at\n",
    "5%. For that we use the so-called Bonferroni correction:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, threshold = map_threshold(z_map, level=.05, height_control='bonferroni')\n",
    "print('Bonferroni-corrected, p<0.05 threshold: %.3f' % threshold)\n",
    "plot_stat_map(z_map, bg_img=mean_img, threshold=threshold,\n",
    "              display_mode='z', cut_coords=3, black_bg=True,\n",
    "              title='active - Finger (p<0.05, corrected)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is quite conservative indeed !  A popular alternative is to\n",
    "control the false discovery rate, i.e. the expected proportion of\n",
    "false discoveries among detections. This is called the false\n",
    "disovery rate.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, threshold = map_threshold(z_map, level=.05, height_control='fdr')\n",
    "print('False Discovery rate = 0.05 threshold: %.3f' % threshold)\n",
    "plot_stat_map(z_map, bg_img=mean_img, threshold=threshold,\n",
    "              display_mode='z', cut_coords=3, black_bg=True,\n",
    "              title='active - Finger (fdr=0.05)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally people like to discard isolated voxels (aka \"small\n",
    "clusters\") from these images. It is possible to generate a\n",
    "thresholded map with small clusters removed by providing a\n",
    "cluster_threshold argument. here clusters smaller than 10 voxels\n",
    "will be discarded.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clean_map, threshold = map_threshold(\n",
    "    z_map, level=.05, height_control='fdr', cluster_threshold=10)\n",
    "plot_stat_map(clean_map, bg_img=mean_img, threshold=threshold,\n",
    "              display_mode='z', cut_coords=3, black_bg=True,\n",
    "              title='active - Finger (fdr=0.05), clusters > 10 voxels')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can save the effect and zscore maps to the disk\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_map.to_filename(join(outdir, 'active_finger_z_map.nii.gz'))\n",
    "eff_map.to_filename(join(outdir, 'active_finger_eff_map.nii.gz'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Report the found positions in a table\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nistats.reporting import get_clusters_table\n",
    "table = get_clusters_table(z_map, stat_threshold=threshold,\n",
    "                           cluster_threshold=20)\n",
    "print(table)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This table can be saved for future use:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "table.to_csv(join(outdir, 'table.csv'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or use [atlasreader](https://github.com/miykael/atlasreader) to get even more information and informative figures:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from atlasreader import create_output\n",
    "create_output(join(outdir, 'active_finger_z_map.nii.gz'), cluster_extent=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at the csv file containing relevant information about the peak of each cluster. This table contains the cluster association and location of each peak, its signal value at this location, the cluster extent (in mm, not in number of voxels), as well as the membership of each peak, given a particular atlas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "peak_info = pd.read_csv('results/active_finger_z_map_peaks.csv')\n",
    "peak_info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the clusters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "cluster_info = pd.read_csv('results/active_finger_z_map_clusters.csv')\n",
    "cluster_info"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each cluster, we also get a corresponding visualization, saved as `.png`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "Image(\"results/active_finger_z_map.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image(\"results/active_finger_z_map_cluster01.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image(\"results/active_finger_z_map_cluster02.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Performing an F-test\n",
    "\n",
    "\"active vs rest\" is a typical t test: condition versus\n",
    "baseline. Another popular type of test is an F test in which one\n",
    "seeks whether a certain combination of conditions (possibly two-,\n",
    "three- or higher-dimensional) explains a significant proportion of\n",
    "the signal.  Here one might for instance test which voxels are well\n",
    "explained by combination of the active and rest condition.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "effects_of_interest = np.vstack((conditions['active - Finger'], conditions['active - Lips']))\n",
    "plot_contrast_matrix(effects_of_interest, design_matrix)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specify the contrast and compute the correspoding map. Actually, the\n",
    "contrast specification is done exactly the same way as for t\n",
    "contrasts.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_map = fmri_glm.compute_contrast(effects_of_interest,\n",
    "                                  output_type='z_score')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the statistic has been converted to a z-variable, which\n",
    "makes it easier to represent it.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clean_map, threshold = map_threshold(\n",
    "    z_map, level=.05, height_control='fdr', cluster_threshold=10)\n",
    "plot_stat_map(clean_map, bg_img=mean_img, threshold=threshold,\n",
    "              display_mode='z', cut_coords=3, black_bg=True,\n",
    "              title='Effects of interest (fdr=0.05), clusters > 10 voxels')\n",
    "plt.show()"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}