2_glm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GLM\n",
    "\n",
    "This notebook provides brief instruction describing how to run general linear model (GLM) analysis on fMRI data with Python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Video lectures on GLM\n",
    "\n",
    "To learn more about GLM, please watch the lecture videos below.\n",
    "\n",
    "- [Principles of fMRI Part 1 Module 15: The General Linear Model - Intro - YouTube](https://www.youtube.com/watch?v=GDkLQuV4he4&list=PLfXA4opIOVrGHncHRxI3Qa5GeCSudwmxM)\n",
    "- [Principles of fMRI Part 1, Module 16: GLM applied to fMRI - YouTube](https://www.youtube.com/watch?v=OyLKMb9FNhg&list=PLfXA4opIOVrGHncHRxI3Qa5GeCSudwmxM)\n",
    "- [Principles of fMRI Part 1, Module 17: Model Building I – conditions and contrasts - YouTube](https://www.youtube.com/watch?v=7MibM1ATai4&list=PLfXA4opIOVrGHncHRxI3Qa5GeCSudwmxM)\n",
    "- [Principles of fMRI Part 1, Module 18: Model Building II – temporal basis sets - YouTube](https://www.youtube.com/watch?v=YfeMIcDWwko&list=PLfXA4opIOVrGHncHRxI3Qa5GeCSudwmxM)\n",
    "- [Principles of fMRI Part 1, Module 19: Model Building III- nuisance variables - YouTube](https://www.youtube.com/watch?v=DEtwsFdFwYc&list=PLfXA4opIOVrGHncHRxI3Qa5GeCSudwmxM)\n",
    "- [Principles of fMRI Part 1, Module 20: GLM Estimation - YouTube](https://www.youtube.com/watch?v=Ab-5AbJ8gAs&list=PLfXA4opIOVrGHncHRxI3Qa5GeCSudwmxM)\n",
    "- [Principles of fMRI Part 1, Module 21: Noise Models- AR models - YouTube](https://www.youtube.com/watch?v=Mb9LDzvhecY&list=PLfXA4opIOVrGHncHRxI3Qa5GeCSudwmxM)\n",
    "- [Principles of fMRI Part 1, Module 22: Inference- Contrasts and t-tests - YouTube](https://www.youtube.com/watch?v=NRunOo7EKD8&list=PLfXA4opIOVrGHncHRxI3Qa5GeCSudwmxM)\n",
    "- [Principles of fMRI Part 1, Module 26: Multiple Comparisons - YouTube](https://www.youtube.com/watch?v=AalIM9-5-Pk&list=PLfXA4opIOVrGHncHRxI3Qa5GeCSudwmxM)\n",
    "- [Principles of fMRI Part 1, Module 27: FWER Correction - YouTube](https://www.youtube.com/watch?v=MxQeEdVNihg&list=PLfXA4opIOVrGHncHRxI3Qa5GeCSudwmxM)\n",
    "- [Principles of fMRI Part 1, Module 28: FDR Correction - YouTube](https://www.youtube.com/watch?v=W9ogBO4GEzA&list=PLfXA4opIOVrGHncHRxI3Qa5GeCSudwmxM)\n",
    "- [Principles of fMRI Part 1, Module 29: More about multiple comparisons - YouTube](https://www.youtube.com/watch?v=N7Iittt8HrU&list=PLfXA4opIOVrGHncHRxI3Qa5GeCSudwmxM)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*For Kamitani lab members:*\n",
    "\n",
    "Please skip the following cell if you run this notebook in Kamitani lab servers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install numpy\n",
    "!pip install nipy\n",
    "!pip install nilearn\n",
    "!pip install git+https://github.com/KamitaniLab/bdpy.git"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Subject 2\n",
    "!wget https://ndownloader.figshare.com/files/28089525\\?private_link=3bd9a1c29f19649c8c0d -O data/sub-02_task-localizer_bold_preproc_native.h5\n",
    "!wget https://ndownloader.figshare.com/files/28089570\\?private_link=3bd9a1c29f19649c8c0d -O data/sub-02_anatomy_t1.nii.gz\n",
    "!wget https://ndownloader.figshare.com/files/28089528\\?private_link=3bd9a1c29f19649c8c0d -O data/sub-02_template_native.nii.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Subject 3\n",
    "!wget https://ndownloader.figshare.com/files/28089534\\?private_link=3bd9a1c29f19649c8c0d -O data/sub-03_task-localizer_bold_preproc_native.h5\n",
    "!wget https://ndownloader.figshare.com/files/28089582\\?private_link=3bd9a1c29f19649c8c0d -O data/sub-03_anatomy_t1.nii.gz\n",
    "!wget https://ndownloader.figshare.com/files/28089537\\?private_link=3bd9a1c29f19649c8c0d -O data/sub-03_template_native.nii.gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from itertools import product\n",
    "\n",
    "import bdpy\n",
    "from bdpy.mri import export_brain_image\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from nipy.modalities.fmri.glm import GeneralLinearModel\n",
    "from nipy.modalities.fmri.experimental_paradigm import BlockParadigm\n",
    "from nipy.modalities.fmri.design_matrix import make_dmtx\n",
    "import nibabel\n",
    "from nilearn import plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## fMRI data\n",
    "\n",
    "In this example, we run GLM analysis on fMRI data collected in the **higher visual areas localizer experiment**.\n",
    "The aim of the experiment is to identify visual areas related to processing of complex visual infromation such as objects, faces, or scenes.\n",
    "During the experiment, a subject was required to look at an image presented in the scanner. The image was either of them.\n",
    "\n",
    "- Object images (IntactObject)\n",
    "- Scrambled version of the object images (ScrambledObject)\n",
    "- Face images (IntaceFace)\n",
    "- Scrambled version of the face images (ScrambledFace)\n",
    "- Scene Images (IntaceScene)\n",
    "- Scrambled version of the scene images (ScrambledFace)\n",
    "\n",
    "Example images:\n",
    "\n",
    "<img src=\"https://drive.google.com/uc?export=view&id=1zyctUwtSJIWEL2xoT6ihkNMwOwgrE7mc&usp=sharing\" width=\"200px\"> <img src=\"https://drive.google.com/uc?export=view&id=1jR9-2uWJGHpJKN_CWucwYZ3ixT8nYbiT&usp=sharing\" width=\"200px\">\n",
    "\n",
    "<img src=\"https://drive.google.com/uc?export=view&id=1iXqsH8BOuTOf2bg-VSiyP5LSql7mcO9u&usp=sharing\" width=\"200px\"> <img src=\"https://drive.google.com/uc?export=view&id=1xiwe2hETdm1QaLZQYxCk_fvMWwK2PyeJ&usp=sharing\" width=\"200px\">\n",
    "\n",
    "<img src=\"https://drive.google.com/uc?export=view&id=1kCd9w-wzhI_88ahT1ubajfJefQiIcrap&usp=sharing\" width=\"200px\"> <img src=\"https://drive.google.com/uc?export=view&id=1aRW6_4DnTKAic4k4kp-sOiebdKOfJnmf&usp=sharing\" width=\"200px\">\n",
    "\n",
    "By taking contrast between intact and scrambled images of a particular domain (e.g., IntactObject vs. ScrambledObject), we can get a brain regions related to the visual processing of the domain (e.g., lateral occipital complex for object vision).\n",
    "\n",
    "Here is a list of typical brain regions identified with the experiment.\n",
    "\n",
    "- LOC (lateral occipital complex; related to object recognition; Malach et al., 1995)\n",
    "- FFA (fusiform face area; related to face perception; Kanwisher et al., 1997)\n",
    "- PPA (parahipocampal place area; related to scene perception; Epstein & Kanwisher, 1998)\n",
    "\n",
    "<img src=\"https://prosopagnosiabytmayo.files.wordpress.com/2014/11/ffa.gif\" width=\"800px\">\n",
    "\n",
    "From: Prosopagnosia | Fusiform Face Area Functions and Impact <https://prosopagnosiabytmayo.wordpress.com/>\n",
    "\n",
    "\n",
    "An experiment that identifies a brain regiton (region-of-interest; ROI) based on the function of the brain (e.g., how voxels are activated by certain stimuli) is called an *functional localizer experiment*.\n",
    "\n",
    "### Experimet information\n",
    "\n",
    "- TR: 3 sec\n",
    "- 100 volumes/run (300 sec/run)\n",
    "- 8 run\n",
    "\n",
    "```\n",
    "Blank (24 s)  Stim. (15 s)  Stim. (15 s)  Blank (15 s)                             Blank (6 s)\n",
    "+----------+  +----------+  +----------+  +----------+                             +----------+\n",
    "|          |  |  Intact  |  |Scrambled |  |          |                             |          |\n",
    "|          |  |  Object  |  |   Face   |  |          |  ... (repated 6 times) ...  |          |\n",
    "+----------+  +----------+  +----------+  +----------+                             +----------+\n",
    "```\n",
    "\n",
    "- Stimulus presentation: different images in the category were flashed in every 500 ms with blank.\n",
    "\n",
    "### fMRI data format in this example\n",
    "\n",
    "The fMRI data are saved as *BData*, machine-learning analysis oritented data format developed in Kamitani Lab. [bdpy](https://github.com/KamitaniLab/bdpy) is required to read the BData."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bdata = bdpy.BData('data/sub-02_task-localizer_bold_preproc_native.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get fMRI data\n",
    "fmri_data = bdata.select('VoxelData')\n",
    "fmri_data.shape  # volumes x voxels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the fMRI data are composed of 800 volumes and  32028 voxels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get labels\n",
    "labels = bdata.get_label('block_type')\n",
    "np.unique(labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data have nine events. Three of them ('PreRest', 'InterRest', and 'PostRest') are 'rest' event in which no visual stimulus was presented. The rest event is not included in the GLM model.\n",
    "\n",
    "In the other six events ('IntactFace', 'IntactObject', 'IntactScene', 'ScrambledFace', and 'ScrambledObject'), a visual stimulus was presented. These 'task' event should be included in the GLM model.\n",
    "\n",
    "The labels are used to create *task regressors*, which model the signal changes caused by the experimental conditions (e.g., stimuli)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get run numbers\n",
    "runs = bdata.select('Run')\n",
    "runs.shape  # samplex x 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is an array including run numbers, used for create *run regressors*, which explicitly model run-wise fluctuations of the fMRI signals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GLM setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Config\n",
    "tr = 3 # TR in sec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Construct the design matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating run regressors\n",
    "def make_runregressors(runs):\n",
    "    '''Returns run regressors and labels.'''\n",
    "    run_reg = np.zeros((runs.shape[0], len(np.unique(runs))))\n",
    "    run_reg_labels = []\n",
    "    \n",
    "    for i, r in enumerate(np.unique(runs)):\n",
    "        run_reg_labels.append('run-%02d' % r)\n",
    "        run_reg[runs.flatten() == r, i] = 1\n",
    "    \n",
    "    return run_reg, run_reg_labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating design matrix\n",
    "\n",
    "ignore_event = ['PreRest', 'InterRest', 'PostRest']\n",
    "\n",
    "condition, onset, duration = [], [], []\n",
    "last_onset = 0\n",
    "i = 0\n",
    "tr = 3\n",
    "while i < len(labels):\n",
    "    lb = labels[i]\n",
    "    nvol = 0\n",
    "    while i < len(labels) and lb == labels[i]:\n",
    "        nvol += 1\n",
    "        i += 1\n",
    "    cn = lb\n",
    "    on = last_onset\n",
    "    dr = nvol * tr\n",
    "    last_onset = on + dr\n",
    "    if lb in ignore_event:\n",
    "        continue\n",
    "    condition.append(cn)\n",
    "    onset.append(on)\n",
    "    duration.append(dr)\n",
    "\n",
    "paradigm = BlockParadigm(con_id=condition, onset=onset, duration=duration)\n",
    "\n",
    "run_regressors, run_regressors_labels = make_runregressors(runs)\n",
    "\n",
    "n_total_vols = len(labels)\n",
    "frametimes = np.linspace(0.5 * tr, n_total_vols * tr - 0.5 * tr, n_total_vols)\n",
    "\n",
    "design = make_dmtx(frametimes, paradigm,\n",
    "                   hrf_model='canonical', drift_model='cosine', hfcut=128,\n",
    "                   add_regs=run_regressors,\n",
    "                   add_reg_names=run_regressors_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Checking the design matrix\n",
    "print(design.names)\n",
    "plt.plot(design.matrix[:100, 0:6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first six regressors stand for the experimental conditons.\n",
    "The rest 45 regressors are nuisance regressors, which model signals irrelevant to the experimental conditions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run GLM\n",
    "\n",
    "Note that the fMRI data should be spatially smoothed before fitting a GLM model to make the data satisfying assumption underlying GLM.\n",
    "In this example, we omit the smoothing for simplicity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GLM\n",
    "model = GeneralLinearModel(design.matrix)\n",
    "model.fit(fmri_data, model='ar1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing contrasts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define contrast and get statistics (z score)\n",
    "contrast = [0, 1, 0, 0, -1, 0] + [0] * 45  # Intract object vs scrambled object\n",
    "z_vals = model.contrast(contrast).z_score() # z score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function converting converting z_vals to a brain image\n",
    "def convert_mri_vector_to_3d(mri_vector, voxel_xyz, brain_template):\n",
    "    mri_table = {}\n",
    "    for i in range(mri_vector.shape[0]):\n",
    "        x, y, z = voxel_xyz[0, i], voxel_xyz[1, i], voxel_xyz[2, i]\n",
    "        mri_table.update({(x, y, z): mri_vector[i]})\n",
    "\n",
    "    mri_3d = np.zeros(brain_template.shape)\n",
    "    for i, j, k in product(range(brain_template.shape[0]),\n",
    "                           range(brain_template.shape[1]),\n",
    "                           range(brain_template.shape[2])):\n",
    "        # Voxel index --> coordinates\n",
    "        x, y, z = brain_template.affine[:3, :3].dot([i, j, k]) + brain_template.affine[:3, 3]\n",
    "        if (x, y, z) in mri_table:\n",
    "            mri_3d[i, j, k] = mri_table[(x, y, z)]\n",
    "\n",
    "    return mri_3d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "voxel_xyz = np.vstack([\n",
    "    bdata.get_metadata('voxel_x', where='VoxelData'),\n",
    "    bdata.get_metadata('voxel_y', where='VoxelData'),\n",
    "    bdata.get_metadata('voxel_z', where='VoxelData'),\n",
    "])\n",
    "brain_template = nibabel.load('data/sub-02_template_native.nii.gz')\n",
    "\n",
    "z_vals_3d = convert_mri_vector_to_3d(z_vals, voxel_xyz, brain_template)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(z_vals_3d[:, :, 23])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_vals_img = nibabel.Nifti1Image(z_vals_3d, brain_template.affine)\n",
    "\n",
    "fig = plt.figure(figsize=(11.69, 8.27))\n",
    "display = plotting.plot_anat('data/sub-02_anatomy_t1.nii.gz',\n",
    "                             figure=fig, cut_coords=(-38, -78, 20))\n",
    "display.add_overlay(z_vals_img, threshold=3, colorbar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise tasks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 1**: Visualize different contrast and see which part of the brain is activated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code comes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 2**: [Optional] Try running GLM analysis on the other subject (sub-03) and visualize the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code comes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task 3**: [Optional] Why some voxels that are apparently not related to the task seem to be activated? How we can avoid to get such acitvation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your answer comes here:"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}