notebooks/nb06_ROI_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- **Author:** [Dace Apšvalka](https://www.mrc-cbu.cam.ac.uk/people/dace.apsvalka/) \n",
    "- **Date:** February 2025  \n",
    "- **conda environment**: [mri conda environment](https://github.com/dcdace/fMRI_training/blob/main/mri_environment.yml)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# fMRI Data Analysis: Region Of Interest (ROI) Analysis\n",
    "\n",
    "Along with whole-brain analysis, it is common to perform region-of-interest (ROI) analysis on fMRI data. This involves extracting signals from specified ROIs. ROI analysis is valuable for several reasons, including **data exploration**, **statistical control**, and **functional specification**. For more information, see this brief review paper by R. Poldrack (2007): https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2555436/.\n",
    "\n",
    "ROIs can be defined based on brain structure or function. For **structural ROIs**, it is recommended to define them based on each subject's anatomy wherever possible. However, this is not always possible, or feasible. In such cases, one can use ROIs based on probabilistic or (less preferable) deterministic atlases. **Functional ROIs** can be defined from 'localiser' scans that identify specific regions exhibiting a particular response (e.g., 'face area' in the fusiform gyrus that respond more strongly to faces than to non-faces). Alternatively, functional ROIs can be defined using **independent contrasts** (make sure you don't do ['double dipping'](https://www.nature.com/articles/nn.2303)) within a factorial design or derived **from previous studies**.\n",
    "\n",
    "In this tutorial, we demonstrate three common types of ROIs: a single voxel, a sphere around a voxel, and a ROI from a probabilistic atlas.\n",
    "\n",
    "1. **A single voxel** is used when the hypothesis is very precise, aiming to identify a specific point in the brain that corresponds to a known functional or anatomical feature. When identifying peak activity within a region, reporting the single voxel with the maximum signal can be informative. There is no risk of averaging signals from different functional areas, which can occur when using larger ROIs. However, single voxels are more susceptible to noise and may represent spurious activity, leading to less reliable results. Additionally, findings based on a single voxel may be less generalisable.\n",
    "\n",
    "2. **A spherical ROI** can help to average out noise and increase the reliability of the signal, as well as generalisability across studies. However, averaging across multiple voxels can reduce anatomical and functional specificity, potentially blurring the precise location of the activity.\n",
    "\n",
    "3. **A probabilistic atlas** is commonly used when investigating brain regions that are functionally or anatomically well-characterised in the literature and show consistent properties across individuals. Using a probabilistic atlas allows for more straightforward comparisons between studies, as the ROIs are defined using the same criteria across different datasets. This enhances the reproducibility and comparability of results. Probabilistic atlases help reduce subjective decisions about where to place ROIs, thereby minimising experimenter bias. Unlike deterministic atlases (which are also commonly used), probabilistic atlases account for some variability in the exact location across subjects. However, they still represent an average across individuals, which might not fully capture the specific anatomy or function of a region in a particular subject. \n",
    "\n",
    "It is common, and often advisable, to perform **ROI analysis in a subject's native space** (using their original brain images without normalisation). In this tutorial, we will also provide an example of how to transform an MNI-space ROI to a subject's native space for native-space analysis.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Table of contents**    \n",
    "1. Import required packages and set up some stuff   \n",
    "2. Load the subject-level fMRI results    \n",
    "3. A single voxel ROI    \n",
    "3.1. Function `extract_MNI_voxel_data`    \n",
    "4. A spherical ROI  \n",
    "4.1. Function `create_voxel_sphere`   \n",
    "4.2. Function `create_MNIvoxel_sphere`   \n",
    "4.3. Function `extract_mean_ROI_values`    \n",
    "5. A ROI from a probabilistic atlas   \n",
    "6. Transforming MNI-space ROI to native-space for native-space analysis    \n",
    "7. Statistics and plotting of the ROI data    \n",
    "7.1. Paired t-test: Initial vs Repeated Faces in the Right Fussiform Gyrus    \n",
    "7.2. 3-way RM ANOVA. Initial vs Immediate Repetition vs Delayed Repetition in the Right Fusiform Gyrus    \n",
    "7.2.1. A simple plot    \n",
    "7.2.2. A 'publication-ready' plot    \n",
    "7.3. 2-by-2 RM ANOVA. Repetition effect for Familiarity x Lag in the Right Fusiform Gurys    \n",
    "7.4. Familiar vs Unfamiliar Faces in Amygdala    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----\n",
    "## Import required packages and set up some stuff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The conda environment used for this tutorial is available here: https://github.com/MRC-CBU/COGNESTIC/blob/main/mri_environment.yml \n",
    "\n",
    "import os.path as op # for file path operations\n",
    "import pandas as pd # for data manipulation\n",
    "import numpy as np # for numerical operations\n",
    "\n",
    "# Plotting libraries\n",
    "import matplotlib.pyplot as plt # basic plots\n",
    "import seaborn as sns # better looking plots\n",
    "\n",
    "# BIDS library\n",
    "from bids.layout import BIDSLayout # to query BIDS dataset\n",
    "\n",
    "# Stats libraries\n",
    "from scipy import stats \n",
    "from statsmodels.stats.anova import AnovaRM\n",
    "from statsmodels.stats.multicomp import pairwise_tukeyhsd\n",
    "\n",
    "# # Nilearn modules https://nilearn.github.io/\n",
    "from nilearn.datasets import load_mni152_template\n",
    "from nilearn.plotting import plot_roi\n",
    "from nilearn.image import resample_to_img\n",
    "from nilearn.masking import _unmask_3d, compute_brain_mask\n",
    "from nilearn.maskers import nifti_spheres_masker\n",
    "\n",
    "import nibabel as nib # NiBabel, to read and write neuroimaging data, https://nipy.org/nibabel/\n",
    "import ants # ANTsPy (Advanced Normalization Tools), https://github.com/ANTsX/ANTsPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MNI152 template will be used as a backgound for plotting MNI-space ROIs\n",
    "mni152_template = load_mni152_template() \n",
    "# Note that I could just use '= load_mni152_template()' \n",
    "# but I want to easily see where the modules are coming from\n",
    "\n",
    "# ROI MNI coordinates. Taken from the group-level results\n",
    "ROI_label = 'Right Fusiform Gyrus'\n",
    "MNI_coord = (42,\t-54,\t-18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the subject-level fMRI results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# --- Set up the paths to the data and results folders\n",
    "fmri_data_dir = 'example_data/FaceRecognition/data' # data in BIDS format\n",
    "fmri_results_dir = 'example_data/FaceRecognition/results' # results in BIDS format\n",
    "\n",
    "# --- Set up the BIDS layout\n",
    "layout = BIDSLayout(fmri_data_dir, derivatives = True)\n",
    "\n",
    "# Attach the results folder to the layout. It must complay with BIDS standards. \n",
    "# And must include dataset_description.json file!\n",
    "layout.add_derivatives(op.join(fmri_results_dir, \"first-level\"))\n",
    "\n",
    "# --- Get all subject-level effect files for the specified conditions\n",
    "conditions = ['IniFF', 'ImmFF', 'DelFF', 'IniUF', 'ImmUF', 'DelUF']\n",
    "effect_files = layout.get(desc=conditions, suffix='effect', extension='.nii.gz')\n",
    "\n",
    "# Print how many files there are to check if it is correct\n",
    "nsub = len(layout.get_subjects())\n",
    "ncond = len(conditions)\n",
    "print(f\"\\nFound {len(effect_files)} files. Should be {nsub} subjects * {ncond} conditions = {nsub * ncond}\\n\")\n",
    "# Print them for inspection\n",
    "print(*effect_files, sep=\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A single voxel ROI\n",
    "\n",
    "In this section, we demonstrate how to extract voxel-specific data from a set of fMRI result files using given MNI coordinates. It is assumed that the result files are all in MNI space. \n",
    "\n",
    "The goal of the `extract_MNI_voxel_data` function is to loop through a list of fMRI result files, identify the voxel corresponding to the provided MNI coordinates, and extract the activation value at that voxel for each subject and condition. The function will return a table (DataFrame) containing the subject ID, condition name, original MNI coordinates, corresponding voxel coordinates in voxel space, and the extracted voxel value. This data can then be used for further statistical analysis and visualisation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function `extract_MNI_voxel_data`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_MNI_voxel_data(result_files, MNI_coord):\n",
    "    \"\"\"\n",
    "    Extract voxel data from a list of results files at given MNI coordinates.\n",
    "\n",
    "    Parameters:\n",
    "    - result_files: List of BIDSImageFiles, or a list of file names as strings, containing subject-level results in MNI space.\n",
    "    - MNI_coord: Tuple of MNI coordinates (x, y, z).\n",
    "\n",
    "    Returns:\n",
    "    - DataFrame with columns: 'subject', 'condition', 'voxel', 'value'.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Initialize an empty list to store the data\n",
    "    data_list = []\n",
    "\n",
    "    # Loop through the files and extract the values to add to the list\n",
    "    for file in result_files:\n",
    "        # If file is not a BIDSImageFile, but a string, load it\n",
    "        if isinstance(file, str):\n",
    "            file = layout.get_file(file)\n",
    "      \n",
    "        # Get the subject ID\n",
    "        subject = file.get_entities()['subject']\n",
    "        \n",
    "        # Get the condition name\n",
    "        condition = file.get_entities()['desc']\n",
    "        \n",
    "        # Get the voxel coordinates from the MNI coordinates\n",
    "        img = file.get_image() # contains data as a 3D numpy array, affine matrix, header.\n",
    "        inverse_affine = np.linalg.inv(img.affine)\n",
    "        vox_coord = nib.affines.apply_affine(inverse_affine, MNI_coord)\n",
    "        \n",
    "        # Get the voxel value\n",
    "        vox_value = img.get_fdata()[int(vox_coord[0]), int(vox_coord[1]), int(vox_coord[2])]\n",
    "        \n",
    "        # Append the data to the list\n",
    "        data_list.append([subject, condition, MNI_coord, vox_coord, vox_value.item()])\n",
    "    \n",
    "    # Convert the list to a DataFrame\n",
    "    data = pd.DataFrame(data_list, columns=['subject', 'condition', 'MNI', 'voxel', 'value'])\n",
    "    \n",
    "    return data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the function to extract the contrast values of the specified voxel from the subject-level results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_from_MNI_coord = extract_MNI_voxel_data(effect_files, MNI_coord) \n",
    "\n",
    "print(data_from_MNI_coord)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A spherical ROI\n",
    "\n",
    "In this section, we provide two examples of how to create a spherical ROI around a specific voxel in a 3D MRI image and then extract the mean values from this ROI across multiple fMRI result files.\n",
    "\n",
    "* The first function, `create_voxel_sphere`, can be applied to images in any space (native, MNI, etc.). The input voxel coordinates should be in voxel-space. If your image is in MNI-space and you have the MNI coordinates around which you want to create the sphere, you can obtain the corresponding voxel coordinates either programmatically (see the function *extract_MNI_voxel_data* above – *# Get the voxel coordinates from the MNI coordinates*) or by using visualisation tools like FSLeyes. In FSLeyes, you can open the image, enter the MNI coordinates, and it will display the corresponding voxel coordinates. (You can also easily adjust the function to be more flexible, allowing it to accept either MNI or voxel-space coordinates.) The function constructs a spherical mask in voxel-space around a specified centre voxel, with a radius defined in millimetres. If a *brain_mask* is provided, the spherical ROI will be constrained within that mask.\n",
    "\n",
    "* The second function, `create_MNIvoxel_sphere`, uses built-in Nilearn functions to achieve the same result as the first function but expects the input coordinates to be in MNI-space (or another world-space) coordinates.\n",
    "\n",
    "Once the spherical ROI is created, the `extract_mean_ROI_values` function can be used to extract and compute the mean activation values within this ROI across multiple fMRI result files. The extracted data can then be used for statistical analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function `create_voxel_sphere`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_voxel_sphere(base_img, center_voxel, radius_mm, brain_mask=None):\n",
    "  \n",
    "    \"\"\"\n",
    "    Create a spherical mask around a specified voxel in a 3D MRI image.\n",
    "\n",
    "    Parameters:\n",
    "    base_img (nibabel.Nifti1Image): A 3D Nifti1Image that defines the space in which the sphere will be created..\n",
    "    center_voxel (tuple): The (x, y, z) coordinates of the center voxel, in voxel space.\n",
    "    radius_mm (float): The radius of the sphere in millimeters.\n",
    "    brain_mask (nibabel.Nifti1Image): A NIFTI image containing a brain mask. If provided,\n",
    "                                      the sphere mask will be restricted to the brain mask.\n",
    "\n",
    "    Returns:\n",
    "    nibabel.Nifti1Image: A NIfTI image containing a binary mask where voxels\n",
    "                         within the specified radius are set to 1.\n",
    "                             \n",
    "    Raises:\n",
    "    TypeError: If the input image is not a NIfTI image.\n",
    "    ValueError: If the brain mask shape does not match the base_img shape.\n",
    "    \n",
    "    \"\"\"\n",
    "      \n",
    "    # Check if the input image is a NIfTI image\n",
    "    if not isinstance(base_img, nib.Nifti1Image):\n",
    "        raise ValueError(\"Input image must be a nibabel.Nifti1Image image\")\n",
    "      \n",
    "    # Get the data array and spacing from the NIfTI image\n",
    "    data = base_img.get_fdata()\n",
    "   \n",
    "    if data.ndim == 4:  # If data is 4D, take the first volume\n",
    "        data = data[..., 0]\n",
    "        \n",
    "    spacing = base_img.header.get_zooms()\n",
    "    sx, sy, sz = spacing[:3]\n",
    "    x0, y0, z0 = center_voxel\n",
    "    \n",
    "    # Calculate the number of voxels within the maximum distance along each axis\n",
    "    rx = int(np.ceil(radius_mm / sx))\n",
    "    ry = int(np.ceil(radius_mm / sy))\n",
    "    rz = int(np.ceil(radius_mm / sz))\n",
    "    \n",
    "    # Define the search space limits\n",
    "    x_min, x_max = max(0, x0 - rx), min(data.shape[0], x0 + rx + 1)\n",
    "    y_min, y_max = max(0, y0 - ry), min(data.shape[1], y0 + ry + 1)\n",
    "    z_min, z_max = max(0, z0 - rz), min(data.shape[2], z0 + rz + 1)\n",
    "    \n",
    "    # Create a mask with the same shape as the input image\n",
    "    neighbors_mask = np.zeros(data.shape)\n",
    "    \n",
    "    for x in range(x_min, x_max):\n",
    "        for y in range(y_min, y_max):\n",
    "            for z in range(z_min, z_max):\n",
    "                # Calculate the Euclidean distance considering voxel spacing\n",
    "                distance = np.sqrt(\n",
    "                    ((x - x0) * sx) ** 2 +\n",
    "                    ((y - y0) * sy) ** 2 +\n",
    "                    ((z - z0) * sz) ** 2\n",
    "                )\n",
    "                if distance <= radius_mm:\n",
    "                    neighbors_mask[x, y, z] = 1\n",
    "    \n",
    "    # If provided, apply the brain mask to the neighbors mask\n",
    "    if brain_mask is not None:\n",
    "        brain_data = brain_mask.get_fdata()\n",
    "        if brain_data.ndim == 4:\n",
    "          brain_data = brain_data[..., 0]\n",
    "        # Check if the brain mask shape matches the base image shape       \n",
    "        if brain_data.shape != neighbors_mask.shape:\n",
    "            raise ValueError(f\"Brain mask shape must match the base image shape. Got {brain_data.shape} and {neighbors_mask.shape}\")\n",
    "        # Convert the brain mask to a binary mask   \n",
    "        brain_data = np.where(brain_data != 0, 1, 0) \n",
    "        # Apply the brain mask to the neighbors mask\n",
    "        neighbors_mask *= brain_data\n",
    "      \n",
    "    # Create a new NIfTI image using the original image's affine and header\n",
    "    sphere_mask = nib.Nifti1Image(neighbors_mask, base_img.affine, base_img.header)\n",
    "  \n",
    "    return sphere_mask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's apply the function and create a sphere around the MNI coordinates we previously used. We will use one of the effect files (the subject-level result file) as the base image. All our effect files are in the same (MNI) space, so we only need to create one ROI, which will correspond to all our result files. We have already determined the corresponding voxel-space coordinates, and we need to use those."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the NIfTI file\n",
    "base_img = nib.load(effect_files[0]) # The input here can be a BIDSImageFile or a file name string\n",
    "\n",
    "# Center voxel in voxel-space\n",
    "center_voxel = (66, 36, 27)\n",
    "\n",
    "# Radius in mm\n",
    "radius = 10\n",
    "\n",
    "# fMRIprep preprocessed data contain a brain mask which we can use to constrain the sphere mask\n",
    "brain_mask_file = layout.get(return_type='file', datatype='func', suffix='mask', desc='brain', extension='nii.gz')[0]\n",
    "brain_mask = nib.load(brain_mask_file)\n",
    "\n",
    "# Create the sphere\n",
    "sphere_mask = create_voxel_sphere(base_img, center_voxel, radius, brain_mask)\n",
    "\n",
    "# Plot the sphere mask\n",
    "plot_roi(sphere_mask, bg_img=mni152_template, title='Sphere mask', display_mode='ortho', cmap='Set1')\n",
    "\n",
    "# # Save the sphere mask to a file\n",
    "# nib.save(sphere_mask, 'path/where/to/save/sphere_mask.nii.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function `create_MNIvoxel_sphere`\n",
    "\n",
    "This function uses built-in Nilearn functions to achieve the same result as the previous function. It expects the input coordinates to be in MNI-space (or another world-space) coordinates. And the resulting sphere will be constrained to the base_image as a mask. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_MNIvoxel_sphere(base_image, center_voxel, radius_mm):\n",
    "    \"\"\"\n",
    "    Create a sphere mask based on the MNI coordinates and radius within the base_image.\n",
    "\n",
    "    Parameters:\n",
    "    - base_image (nibabel.Nifti1Image): A 3D Nifti1Image (in MNI space) that defines the space in which the sphere will be created.\n",
    "    - center_voxel (tuple): The (x, y, z) coordinates of the center voxel, in MNI space.\n",
    "    - radius_mm (float): Radius of the sphere in millimeters.\n",
    "\n",
    "    Returns:\n",
    "    - sphere_mask (nibabel.Nifti1Image): A Nifti1Image containing the sphere mask.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Compute a brain mask using the povided base image\n",
    "    space_defining_image = compute_brain_mask(base_image)\n",
    "\n",
    "    # Apply mask and get affinity matrix\n",
    "    _, A = nifti_spheres_masker._apply_mask_and_get_affinity(\n",
    "        seeds=[center_voxel],  # in MNI space; can be a list of coordinates\n",
    "        niimg=None,\n",
    "        radius=radius_mm,\n",
    "        allow_overlap=False, \n",
    "        mask_img=space_defining_image\n",
    "    )\n",
    "\n",
    "    # Unmask the sphere mask\n",
    "    sphere_mask_data = _unmask_3d(\n",
    "        X=A.toarray().flatten(), \n",
    "        mask=space_defining_image.get_fdata().astype(bool)\n",
    "    )\n",
    "\n",
    "    # Create a Nifti1Image for the sphere mask\n",
    "    sphere_mask = nib.Nifti1Image(sphere_mask_data, space_defining_image.affine)\n",
    "\n",
    "    return sphere_mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Center voxel in MNI-space\n",
    "center_voxel = MNI_coord\n",
    "\n",
    "# Radius in mm\n",
    "radius = 10\n",
    "\n",
    "# Create the sphere\n",
    "sphere_mask = create_MNIvoxel_sphere(base_img, center_voxel, radius)\n",
    "\n",
    "# Plot the sphere mask\n",
    "plot_roi(sphere_mask, bg_img=mni152_template, title='Sphere mask', display_mode='ortho', cmap='Set1')\n",
    "\n",
    "# # Save the sphere mask to a file\n",
    "# nib.save(sphere_mask, 'path/where/to/save/sphere_mask.nii.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function `extract_mean_ROI_values`\n",
    "\n",
    "Once the spherical ROI is created, we can extract and compute the mean activation values within this ROI across multiple fMRI result files. The ROI and fMRI result files must be in the same space (e.g., MNI-space)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_mean_ROI_values(result_files, ROI):\n",
    "    \"\"\"\n",
    "    Extract mean ROI values from a list of results files using a provided ROI mask.\n",
    "\n",
    "    Parameters:\n",
    "    - result_files: List of BIDSImageFiles, or a list of file names, containing subject-level results.\n",
    "    - ROI: Nifti1Image containing the mask data.\n",
    "\n",
    "    Returns:\n",
    "    - DataFrame with columns: 'subject', 'condition', 'value'.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Check if the mask and the file have the same shape (dimensions) and resample if necessary\n",
    "    first_file = result_files[0]\n",
    "    if isinstance(first_file, str):\n",
    "        first_file = layout.get_file(first_file)\n",
    "    if not ROI.shape == first_file.get_image().shape[:3]:\n",
    "      ROI = resample_to_img(ROI, first_file.get_image(), interpolation='nearest')\n",
    "      print(f\"Resampled the mask to the shape of result files\")\n",
    "        \n",
    "    # Initialize an empty list to store the data\n",
    "    data_list = []\n",
    "\n",
    "    # Loop through the files and extract the mean ROI values\n",
    "    for file in result_files:\n",
    "        # If file is not a BIDSImageFile, but a string, load it\n",
    "        if isinstance(file, str):\n",
    "          file = layout.get_file(file)\n",
    "                 \n",
    "        # Get the subject ID\n",
    "        subject = file.get_entities()['subject']\n",
    "        \n",
    "        # Get the condition name\n",
    "        condition = file.get_entities()['desc']\n",
    "        \n",
    "        # Get the voxel values within the mask\n",
    "        img_data = file.get_image()\n",
    "        mask_data = ROI.get_fdata()\n",
    "        vox_values = img_data.get_fdata()[mask_data > 0]\n",
    "        # Compute the mean value\n",
    "        mean_value = np.mean(vox_values)\n",
    "        \n",
    "        # Append the data to the list\n",
    "        data_list.append([subject, condition, mean_value])\n",
    "    \n",
    "    # Convert the list to a pandas DataFrame\n",
    "    data = pd.DataFrame(data_list, columns=['subject', 'condition', 'value'])\n",
    "    \n",
    "    return data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's apply the function to extract the mean contrast values within the ROI from each results file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_from_sphere = extract_mean_ROI_values(effect_files, sphere_mask)\n",
    "print(data_from_sphere)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A ROI from a probabilistic atlas\n",
    "\n",
    "In this section, we demonstrate how to create an ROI based on probabilistic data from the [Julich-Brain Cytoarchitectonic Atlas](https://search.kg.ebrains.eu/instances/ab191c17-8cd8-4622-aaac-eee11b2fa670). Probabilistic atlases provide voxel-wise information on the likelihood that a particular brain region corresponds to a given area, reflecting inter-individual variability in brain anatomy or function. The atlas is a single 4D image file, with each ROI stored as a separate volume in the file's fourth dimension. It is accompanied by a label file that indicates which volume corresponds to each region.\n",
    "\n",
    "We selected the Julich-Brain Cytoarchitectonic Atlas for this analysis due to its comprehensive coverage of cytoarchitectonic regions, including the amygdala, which is of particular interest in the provided example. The choice of atlas often depends on the specific requirements of the study, such as the need for detailed cytoarchitectonic information or broader anatomical coverage.\n",
    "\n",
    "Here, we focus on generating an ROI for the amygdala by extracting and combining relevant regions from the atlas, applying a probability threshold, and creating a binary mask for use in subsequent analyses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by specifying the location of the atlas's image and lable files, and setting the probability threshold. \n",
    "\n",
    "When using probabilistic atlases to define ROIs, **the choice of threshold** can significantly affect the final ROI. The threshold determines the minimum probability that a voxel belongs to a particular region, so a higher threshold will result in a smaller, more conservative ROI, while a lower threshold will include more voxels and create a larger ROI. The choice of threshold depends on your specific research question, the size and function of the region you are studying, and the need for either more inclusive or more conservative ROIs. A moderate threshold (20-30%) is frequently used as a balance between inclusivity and specificity. This threshold will exclude voxels that are less likely to belong to the region while still capturing the bulk of the structure.\n",
    "\n",
    "It’s advisable to visually inspect the resulting ROIs at different thresholds to ensure they are biologically plausible and appropriate for your analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the location of the atlas and labels files\n",
    "atlas_julich_file = 'example_data/atlases/Julich/JULICH_BRAIN_CYTOARCHITECTONIC_MAPS_2_9_MNI152_2009C_NONL_ASYM.pmaps_resized.nii.gz'\n",
    "labels_julich_file = 'example_data/atlases/Julich/JULICH_BRAIN_CYTOARCHITECTONIC_MAPS_2_9_MNI152_2009C_NONL_ASYM.txt'\n",
    "\n",
    "# Set the probability threshold\n",
    "prob_threshold = 0.25"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we identify which volumes in the atlas correspond to the amygdala, our region of interest. We do this by consulting the accompanying label text file, where each ROI in the atlas is assigned a map index. We search for all labels that contain 'amygdala' in their names and extract their map indexes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert labels txt file to dataframe\n",
    "labels_df = pd.read_csv(labels_julich_file, sep=\" \", quotechar=\"'\", names=['mapindex', 'regionname'], skiprows=1)\n",
    "print(f\"Atlas labels:\\n{labels_df}\\n----------------------------------\")\n",
    "\n",
    "# find the indexes for regions containing amygdala\n",
    "mapindex_of_interest = labels_df[labels_df['regionname'].str.contains('Amygdala')].index\n",
    "print(f\"Labels of interest:\\n{labels_df.loc[mapindex_of_interest]}\\n----------------------------------\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can create a binary mask for the amygdala from probabilistic atlas data, applying a probability threshold to ensure that only regions with sufficient likelihood of corresponding to the amygdala are included in the final mask. We do it in the following steps:\n",
    "\n",
    "1. Load the 4D atlas image.\n",
    "2. Extract the specific 3D volumes corresponding to the amygdala.\n",
    "3. Threshold the selected data to create a binary mask for each 3D volume.\n",
    "4. Combine these 3D images by summing them.\n",
    "5. Binarise the combined 3D data to create the final binary mask for the amygdala.\n",
    "6. Convert the combined data, which is a NumPy array, into a 3D NIfTI image. The affine transformation matrix from the original atlas image is used to ensure the mask aligns correctly in the same spatial orientation as the atlas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the 4D atlas image\n",
    "atlas_img = nib.load(atlas_julich_file, mmap=True)\n",
    "atlas_data = atlas_img.get_fdata()\n",
    "\n",
    "# Select the maps of interest (from the 4th dimension of the 4D image)\n",
    "selected_data = atlas_data[:, :, :, mapindex_of_interest]\n",
    "\n",
    "# Threshold the data\n",
    "thresholded_data = np.where(selected_data > prob_threshold, 1, 0)\n",
    "\n",
    "# Combine these 3D images by summing them\n",
    "new_3d_data = np.sum(thresholded_data, axis=-1)\n",
    "\n",
    "# Binarize the data\n",
    "new_3d_data = np.where(new_3d_data > 0, 1, 0)\n",
    "\n",
    "# Cast the data to int8 to avoid compatibility issues\n",
    "new_3d_data = new_3d_data.astype(np.uint8)\n",
    "\n",
    "# Create a new NIfTI image with the 3D data\n",
    "amygdala_mask = nib.Nifti1Image(new_3d_data, atlas_img.affine) # could set dtype here, but setting it above works better (no warnings)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot the created amygdala ROI and save if necessary. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the amygdala mask\n",
    "plot_roi(\n",
    "  amygdala_mask, \n",
    "  title='Amygdala mask', \n",
    "  display_mode='ortho', \n",
    "  draw_cross=False, \n",
    "  annotate=True, \n",
    "  black_bg=True, \n",
    "  bg_img = mni152_template, \n",
    "  cmap='Set1')\n",
    "\n",
    "# # Save the mask to a file\n",
    "# nib.save(amygdala_mask, '/my_path/amygdala_mask.nii.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And once we have the ROI, we can extract the contrast values using the previously created `extract_mean_ROI_values` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_from_amygdala = extract_mean_ROI_values(effect_files, amygdala_mask)\n",
    "print(data_from_amygdala)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An important aspect to note is that the dimensions of the created amygdala ROI differ from those of the results files. For data extraction to occur correctly, the dimensions of the two must match. We have addressed this in our `extract_mean_ROI_values` function. If the ROI's dimensions do not match those of the results file, the ROI is resampled to align with the results file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check the shape of the amygdala mask and the result files\n",
    "print(\"Amygdala mask shape:\", amygdala_mask.shape[:3])\n",
    "print(\"Result files shape:\", effect_files[0].get_image().shape[:3])\n",
    "\n",
    "# Check the voxel dimensions of the amygdala mask\n",
    "vox_dims = np.diag(amygdala_mask.affine)[:3]\n",
    "print(f\"Voxel dimensions of the amygdala mask: {vox_dims}\")\n",
    "\n",
    "# Check the voxel dimensions of the result files\n",
    "vox_dims = np.diag(effect_files[0].get_image().affine)[:3]\n",
    "print(f\"Voxel dimensions of the result files: {vox_dims}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transforming MNI-space ROI to native-space for native-space analysis\n",
    "\n",
    "In the final example, we demonstrate how to transform an ROI defined in MNI space (e.g., a spherical mask) into a subject's native space using a transformation matrix derived during preprocessing with fMRIPrep. Once transformed, the ROI can be used to extract activation values from the native-space data of individual subjects, which may provide more accurate results than performing the analysis on normalized data.\n",
    "\n",
    "In this example, we will use the ANTs Python package [ANTsPy](https://github.com/ANTsX/ANTsPy) to perform the transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the MNI sphere on subject's native space anatomical image\n",
    "\n",
    "sID = '15'\n",
    "sub_anat_file = layout.get(subject=sID, datatype='anat', suffix='T1w', extension='nii.gz', desc=None, return_type='file')[0]\n",
    "\n",
    "plot_roi(sphere_mask, bg_img=sub_anat_file, title=f\"{ROI_label} sphere mask in native space; incorrect\", display_mode='ortho', dim=-1, cmap='Set1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MNI_to_native_xfm_path = layout.get(subject=sID, datatype='anat', extension='h5', to='T1w', return_type='file')[0]\n",
    "print(f\"Found the following transformation file: {MNI_to_native_xfm_path}\")\n",
    "\n",
    "# save sphere mask to disk\n",
    "MNI_mask_file = 'example_data/FaceRecognition/scratch/sphere_mask.nii.gz'\n",
    "nib.save(sphere_mask, MNI_mask_file)\n",
    "\n",
    "# Apply the transformation to the sphere mask\n",
    "ROI_native_space = ants.apply_transforms(\n",
    "  fixed = ants.image_read(sub_anat_file), \n",
    "  moving = ants.image_read(MNI_mask_file),\n",
    "  transformlist = MNI_to_native_xfm_path, \n",
    "  interpolator = 'nearestNeighbor')\n",
    "\n",
    "# Convert ants image to nibabel image\n",
    "ROI_native_space = nib.Nifti1Image(ROI_native_space.numpy(), nib.load(sub_anat_file).affine)\n",
    "\n",
    "# plot the ROI in native space\n",
    "plot_roi(ROI_native_space, bg_img=sub_anat_file, title=f\"{ROI_label} sphere mask in native space; correct\", display_mode='ortho', dim=-1, cmap='Set1')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Importantly, when working in MNI space, a single ROI mask file is sufficient for the entire group since all subjects are aligned to the same space. However, in native space, each subject requires a unique ROI mask due to the individual differences in brain anatomy. As a result, the `extract_mean_ROI_values` function must be applied separately to each subject's results files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistics and plotting of the ROI data\n",
    "\n",
    "Now that we have created the ROI masks and extracted the data, we can proceed with statistical analyses as we would with any other dataset. Below are some examples of the ROI analysis results for our example dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Paired t-test: Initial vs Repeated Faces in the Right Fussiform Gyrus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data_from_MNI_coord # or data_from_sphere\n",
    "\n",
    "# Paired t-test\n",
    "\n",
    "# Separate the data into Initial and Repetition presentations\n",
    "initial = data[data['condition'].isin(['IniFF', 'IniUF'])].groupby('subject')['value'].mean()\n",
    "repeated = data[data['condition'].isin(['ImmFF', 'DelFF', 'ImmUF', 'DelUF'])].groupby('subject')['value'].mean()\n",
    "\n",
    "# Perform the paired t-test\n",
    "t_stat, p_value = stats.ttest_rel(initial, repeated)\n",
    "print(f\"Initial vs Repeated presentation of faces in the {ROI_label}. \\nPaired t-test results: t-statistic = {t_stat}, p-value = {p_value}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a long format DataFrame for the plot\n",
    "data_initial_repeated = pd.DataFrame({\n",
    "    'condition': ['initial'] * len(initial) + ['repeated'] * len(repeated),\n",
    "    'value': pd.concat([initial, repeated]).reset_index(drop=True),\n",
    "    'subject': list(initial.index) + list(repeated.index)\n",
    "})\n",
    "print(data_initial_repeated)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(4, 4))\n",
    "\n",
    "# Plot the point plot with dodge to avoid overlap\n",
    "sns.violinplot(x='condition', y='value', data=data_initial_repeated, linestyles='-', color='lightgrey')\n",
    "\n",
    "# Add lines\n",
    "for subject in data_initial_repeated['subject'].unique():\n",
    "    subject_data = data_initial_repeated[data_initial_repeated['subject'] == subject]\n",
    "    plt.plot([0, 1], subject_data['value'], color='gray', alpha=0.3)\n",
    "\n",
    "# Remove plot borders\n",
    "plt.gca().spines['top'].set_visible(False)\n",
    "plt.gca().spines['right'].set_visible(False)\n",
    "plt.gca().spines['left'].set_visible(False)\n",
    "plt.gca().spines['bottom'].set_visible(False)\n",
    "\n",
    "# Remove tick marks but keep the labels\n",
    "plt.tick_params(axis='both', which='both', length=0)  # Set tick marks length to 0\n",
    "\n",
    "plt.ylabel('Beta estimate')\n",
    "plt.title(f\"{ROI_label}\\nInitial vs Repeated presentation of faces\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-way RM ANOVA. Initial vs Immediate Repetition vs Delayed Repetition in the Right Fusiform Gyrus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Repeated measures ANOVA\n",
    "\n",
    "initial = data[data['condition'].isin(['IniFF', 'IniUF'])].groupby('subject')['value'].mean()\n",
    "immediate_rep = data[data['condition'].isin(['ImmFF', 'ImmUF'])].groupby('subject')['value'].mean()\n",
    "delayed_rep = data[data['condition'].isin(['DelFF', 'DelUF'])].groupby('subject')['value'].mean()\n",
    "\n",
    "# Create a long format DataFrame\n",
    "data_3way_presentation = pd.DataFrame({\n",
    "    'condition': ['initial'] * len(initial) + ['immediate'] * len(immediate_rep) + ['delayed'] * len(delayed_rep),\n",
    "    'value': pd.concat([initial, immediate_rep, delayed_rep]).reset_index(drop=True),\n",
    "    'subject': list(initial.index) + list(immediate_rep.index) + list(delayed_rep.index)\n",
    "})\n",
    "\n",
    "#print(data_3way_presentation)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Perform repeated measures ANOVA\n",
    "anova = AnovaRM(data_3way_presentation, 'value', 'subject', within=['condition'])\n",
    "anova_results = anova.fit()\n",
    "\n",
    "print(anova_results)\n",
    "\n",
    "tukey = pairwise_tukeyhsd(endog=data_3way_presentation['value'], groups=data_3way_presentation['condition'], alpha=0.05)\n",
    "\n",
    "print(tukey) # Tukey’s HSD: This test compares all possible pairs of means and controls the family-wise error rate.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A simple plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot data_3way_presentation\n",
    "condition_order = ['initial', 'immediate', 'delayed']\n",
    "\n",
    "plt.figure(figsize=(6, 4))\n",
    "sns.violinplot(x='condition', y='value', data=data_3way_presentation, color='lightgrey', order=condition_order)\n",
    "plt.ylabel('Beta estimate')\n",
    "plt.xlabel('Face presentation')\n",
    "plt.title(f\"{ROI_label}\\nInitial vs Immediate rep. vs Delayed rep. presentation of faces\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A 'publication-ready' plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_1by3_anova(data, subject_col, condition_col, value_col, condition_order, \n",
    "                    roi_label='', mni_coord='', title='', xlabel='', ylabel='', xticks=[]):\n",
    "    \"\"\"\n",
    "    Plots the results of a 1-by-3 ANOVA with within-subject 95% confidence intervals and significance annotations.\n",
    "    \n",
    "    Parameters:\n",
    "    - data: DataFrame containing the data (long format).\n",
    "    - subject_col: Name of the column representing subjects.\n",
    "    - condition_col: Name of the column representing the condition.\n",
    "    - value_col: Name of the column representing the dependent variable.\n",
    "    - condition_order: List specifying the order of conditions.\n",
    "    - roi_label: Label for the region of interest (ROI).\n",
    "    - mni_coord: MNI coordinates as a string.\n",
    "    - xlabel: Label for the x-axis.\n",
    "    - ylabel: Label for the y-axis.\n",
    "    - xticks: List of custom x-axis labels.\n",
    "    - title: Title for the plot.\n",
    "    \n",
    "    Returns:\n",
    "    - A plot of the ANOVA results with error bars and significance markers.\n",
    "    \"\"\"    \n",
    "    # Perform ANOVA\n",
    "    anova = AnovaRM(data_3way_presentation, 'value', 'subject', within=['condition'])\n",
    "    anova_results = anova.fit()\n",
    "    \n",
    "    # Get the ANOVA F-value and p-value for the plot title\n",
    "    F_value = anova_results.anova_table['F Value'].iloc[0]\n",
    "    \n",
    "    if anova_results.anova_table['Pr > F'].iloc[0] < 0.001:\n",
    "        p_value = \"< 0.001\"\n",
    "    else:\n",
    "        p_value = f\"= {anova_results.anova_table['Pr > F'].iloc[0]:.3f}\"\n",
    "        \n",
    "    # Perform Tukey's HSD test\n",
    "    tukey_results = pairwise_tukeyhsd(endog=data_3way_presentation['value'], groups=data_3way_presentation['condition'], alpha=0.05)\n",
    "    print(tukey_results)\n",
    "    \n",
    "    # Determine the significance markers for the pairs based on the Tukey's HSD test\n",
    "    annotations = ['***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.' for p in tukey_results.pvalues]\n",
    "    \n",
    "    # --------------------------------------------------------------\n",
    "    # Match pairs to annotations based on the Tukey's HSD test results\n",
    "    # --------------------------------------------------------------\n",
    "    # Get condition labels from Tukey results\n",
    "    tukey_conditions = np.array(tukey_results.groupsunique) \n",
    "\n",
    "    # Map conditions to their indices in the user-specified order\n",
    "    condition_indices = {cond: i for i, cond in enumerate(condition_order)}\n",
    "\n",
    "    # Extract and sort pairs according to the specified condition order\n",
    "    pair_indices = tukey_results._multicomp.pairindices \n",
    "    \n",
    "    group1_idx, group2_idx = pair_indices \n",
    "    tukey_pairs = list(zip(tukey_conditions[group1_idx], tukey_conditions[group2_idx]))\n",
    "\n",
    "    # Map Tukey pairs to specified condition indices\n",
    "    pairs = [(condition_indices[g1], condition_indices[g2]) for g1, g2 in tukey_pairs]\n",
    "\n",
    "    # --------------------------------------------------------------\n",
    "    \n",
    "    # Ensure 'subject' and 'condition' are categorical variables\n",
    "    data[subject_col] = data[subject_col].astype('category')\n",
    "    data[condition_col] = data[condition_col].astype('category')\n",
    "    \n",
    "    # Calculate the mean for each condition\n",
    "    means = data.groupby(condition_col, observed=True)[value_col].mean()\n",
    "    \n",
    "    # Subtract the participant's mean value from each score to calculate deviations\n",
    "    data['deviation'] = data[value_col] - data.groupby(subject_col, observed=True)[value_col].transform('mean')\n",
    "    \n",
    "    # Calculate the standard error of the mean for each condition\n",
    "    sem = data.groupby(condition_col, observed=True)['deviation'].std() / np.sqrt(data[subject_col].nunique())\n",
    "    \n",
    "    # Calculate the t-value for the 95% confidence interval\n",
    "    t_value = stats.t.ppf(1 - 0.025, data[subject_col].nunique() - 1)  # degrees of freedom = n-1\n",
    "    \n",
    "    # Calculate the confidence interval\n",
    "    ci = sem * t_value\n",
    "    \n",
    "    # Plot the data\n",
    "    plt.figure(figsize=(6, 8))\n",
    "    \n",
    "    # Create a point plot\n",
    "    sns.pointplot(data=data, x=condition_col, y=value_col, order=condition_order, capsize=.1, errorbar=None, color='black', markers='o', linestyles='-', estimator=np.mean)\n",
    "    \n",
    "    # Add error bars using the calculated within-subject CI\n",
    "    plt.errorbar(x=range(len(means)), y=means[condition_order], yerr=ci[condition_order], fmt='none', capsize=5, color='black')\n",
    "    \n",
    "    # Add individual data points\n",
    "    sns.stripplot(data=data, x=condition_col, y=value_col, order=condition_order, color='grey', alpha=0.6, jitter=True)\n",
    "    \n",
    "    # Annotate the plot with significance markers\n",
    "    y_max = data[value_col].max() + (data[value_col].max() * 0.04)  # Position above the highest point\n",
    "    y_range = data[value_col].max() + (data[value_col].max() * 0.07)  # Adjust the range to avoid overlapping\n",
    "    \n",
    "    for i, (pair, annotation) in enumerate(zip(pairs, annotations)):\n",
    "        x1, x2 = pair\n",
    "        y = y_max + (i * 0.05)  # Space annotations vertically\n",
    "        plt.plot([x1, x2], [y, y], color='black', lw=1.5)\n",
    "        plt.text((x1 + x2) * 0.5, y + 0.01, annotation, ha='center', va='bottom', color='black')\n",
    "    \n",
    "    # Add horizontal line at 0\n",
    "    plt.axhline(y=0, color='grey')\n",
    "    \n",
    "    # Remove plot borders\n",
    "    plt.gca().spines['top'].set_visible(False)\n",
    "    plt.gca().spines['right'].set_visible(False)\n",
    "    plt.gca().spines['left'].set_visible(False)\n",
    "    plt.gca().spines['bottom'].set_visible(False)\n",
    "    \n",
    "    # Remove tick marks but keep the labels\n",
    "    plt.tick_params(axis='both', which='both', length=0)  \n",
    "    \n",
    "    # Add labels and title\n",
    "    plt.xlabel(f'\\n{xlabel}')\n",
    "    plt.ylabel(ylabel)\n",
    "    \n",
    "    # Custom x labels\n",
    "    plt.xticks(ticks=range(len(means)), labels=xticks)\n",
    "    \n",
    "    # Set title\n",
    "    plt.title(f\"{title} \\n F = {F_value:.2f}, p {p_value}; error bars: within-subj 95% CI\")\n",
    "    \n",
    "    # Show plot\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_1by3_anova(\n",
    "    data=data_3way_presentation,\n",
    "    subject_col='subject',\n",
    "    condition_col='condition',\n",
    "    value_col='value',\n",
    "    condition_order=['initial', 'immediate', 'delayed'],\n",
    "    roi_label=ROI_label,\n",
    "    mni_coord=MNI_coord,\n",
    "    xlabel='Face presentation',\n",
    "    ylabel='BOLD signal (beta estimate)', \n",
    "    xticks=['Initial', 'Immediate\\nrepetition', 'Delayed\\nrepetition'], \n",
    "    title=f\"{ROI_label}, MNI = {MNI_coord}\"\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-by-2 RM ANOVA. Repetition effect for Familiarity x Lag in the Right Fusiform Gurys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> \n",
    "> *\"Behavioural priming has been found for both familiar and unfamiliar faces when repeated immediately, but only for familiar faces when repeated after a delay.\"* [(Henson et al., 2004)](https://www.mrc-cbu.cam.ac.uk/personal/rik.henson/personal/HensonEtAl_HBM_Abstract_04.pdf)\n",
    "> \n",
    "\n",
    "Do we find this in our fMRI results? \n",
    "\n",
    "**Do repetition effect for Familiarity (Famous, Unfamiliar) X Lag (Immediate, Delayed) ANOVA**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pivot the DataFrame and reset the (subject) index\n",
    "df_pivot = data.pivot(index='subject', columns='condition', values='value').reset_index()\n",
    "\n",
    "# Get the repetition suppression effect for each repetition condition for each subject.\n",
    "df_pivot['DelFF'] -= df_pivot['IniFF']\n",
    "df_pivot['ImmFF'] -= df_pivot['IniFF']\n",
    "df_pivot['DelUF'] -= df_pivot['IniUF']\n",
    "df_pivot['ImmUF'] -= df_pivot['IniUF']\n",
    "\n",
    "# Reshape the DataFrame to long format\n",
    "df_rs = df_pivot.melt(id_vars=['subject'], \n",
    "                      value_vars=['DelFF', 'ImmFF', 'DelUF', 'ImmUF'], \n",
    "                      var_name='condition', \n",
    "                      value_name='repetition_effect')\n",
    "\n",
    "# Extracting factors from the condition names\n",
    "df_rs['familiarity'] = df_rs['condition'].apply(lambda x: 'Famous' if 'FF' in x else 'Unfamiliar')\n",
    "df_rs['lag'] = df_rs['condition'].apply(lambda x: 'Immediate' if 'Imm' in x else 'Delayed')\n",
    "\n",
    "# Print or save the resulting DataFrame\n",
    "print(df_rs)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform the repeated measures ANOVA\n",
    "aovrm = AnovaRM(df_rs, depvar='repetition_effect', subject='subject', within=['familiarity', 'lag'])\n",
    "res = aovrm.fit()\n",
    "\n",
    "# Print the results\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_nbyn_anova(data, subject_col, condition_cols, value_col, \n",
    "                    condition_order=None, line_colors=None,  \n",
    "                    xlabel='', ylabel='', title=''):\n",
    "    \"\"\"\n",
    "    Plots the results of an n-by-n ANOVA with within-subject 95% confidence intervals.\n",
    "    \n",
    "    Parameters:\n",
    "    - data: DataFrame containing the data.\n",
    "    - subject_col: Name of the column representing subjects.\n",
    "    - condition_cols: List of columns representing the conditions (e.g., ['familiarity', 'presentation']).\n",
    "    - value_col: Name of the column representing the dependent variable.\n",
    "    - condition_order: List specifying the order of levels within each condition.\n",
    "    - line_colors: List of colors for each condition.\n",
    "    - title: Title for the plot.\n",
    "    - xlabel: Label for the x-axis.\n",
    "    - ylabel: Label for the y-axis.\n",
    "    - title: Title for the plot.\n",
    "    \n",
    "    Returns:\n",
    "    - A plot of the ANOVA results with error bars.\n",
    "    \"\"\"\n",
    "    \n",
    "    # --------------------------------------------------------------------------------------------\n",
    "    # Calculate within-subject 95% confidence intervals\n",
    "    # --------------------------------------------------------------------------------------------\n",
    "\n",
    "    # Ensure subject and condition columns are categorical\n",
    "    data[subject_col] = data[subject_col].astype('category')\n",
    "    for condition_col in condition_cols:\n",
    "        data[condition_col] = data[condition_col].astype('category')\n",
    "    \n",
    "    # Calculate the mean for each condition combination\n",
    "    means = data.groupby(condition_cols, observed=True)[value_col].mean()\n",
    "\n",
    "    # Subtract the participant's mean value from each score to calculate deviations\n",
    "    data['deviation'] = data[value_col] - data.groupby(subject_col, observed=True)[value_col].transform('mean')\n",
    "\n",
    "    # Calculate the standard error of the mean for each condition combination\n",
    "    sem = data.groupby(condition_cols, observed=True)['deviation'].std() / np.sqrt(data[subject_col].nunique())\n",
    "\n",
    "    # Calculate the t-value for the 95% confidence interval\n",
    "    t_value = stats.t.ppf(1 - 0.025, data[subject_col].nunique() - 1)\n",
    "\n",
    "    # Calculate the confidence interval\n",
    "    ci = sem * t_value\n",
    "    \n",
    "    \n",
    "    # Print the means and confidence intervals\n",
    "    print(\"\\n----------------------------------------\\n\")\n",
    "    print(\"Condition means:\" + \"\\n\" + str(means))\n",
    "    print(\"\\n----------------------------------------\\n\")\n",
    "    print(\"Condition CIs:\" + \"\\n\" + str(ci))\n",
    "    print(\"\\n----------------------------------------\\n\")\n",
    "\n",
    "    # --------------------------------------------------------------------------------------------\n",
    "    # Plot the data\n",
    "    # --------------------------------------------------------------------------------------------\n",
    "    plt.figure(figsize=(4, 4))\n",
    "    \n",
    "    ax = sns.pointplot(\n",
    "        data=data, x=condition_cols[-1], y=value_col, hue=condition_cols[0],\n",
    "        errorbar=None, palette=line_colors, order=condition_order, dodge=0.02\n",
    "    )\n",
    "\n",
    "    # Extract x positions for error bars\n",
    "    positions = [x_value for line in ax.lines for x_value in line.get_xdata()]\n",
    "\n",
    "    # Get the legend from the axis to know which line corresponds to which condition\n",
    "    legend = ax.get_legend()\n",
    "    legend_labels = [text.get_text() for text in legend.get_texts()]\n",
    "\n",
    "    # Extract means and error values\n",
    "    mean_values = []\n",
    "    error_values = []\n",
    "    for label in legend_labels:\n",
    "        for condition in condition_order:\n",
    "            mean_values.append(means[label, condition])\n",
    "            error_values.append(ci[label, condition])\n",
    "\n",
    "    # Plotting the error bars in a loop\n",
    "    for i in range(len(positions)):\n",
    "        plt.errorbar(x=positions[i], y=mean_values[i], yerr=error_values[i],\n",
    "                     fmt='none', capsize=5, color=line_colors[i // len(condition_order)])\n",
    "\n",
    "    # Remove plot borders\n",
    "    plt.gca().spines['top'].set_visible(False)\n",
    "    plt.gca().spines['right'].set_visible(False)\n",
    "    plt.gca().spines['left'].set_visible(False)\n",
    "    plt.gca().spines['bottom'].set_visible(False)\n",
    "\n",
    "    # Remove tick marks but keep the labels\n",
    "    plt.tick_params(axis='both', which='both', length=0)\n",
    "\n",
    "    # Remove the legend title\n",
    "    plt.gca().legend().set_title(None)\n",
    "\n",
    "    # Add horizontal line at 0\n",
    "    plt.axhline(y=0, color='black', linewidth=0.5)\n",
    "\n",
    "    # Add labels and title\n",
    "    plt.title(title)\n",
    "    plt.ylabel(ylabel)\n",
    "    plt.xlabel(xlabel)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_nbyn_anova(\n",
    "    data=df_rs,\n",
    "    subject_col='subject',\n",
    "    condition_cols=['familiarity', 'lag'],\n",
    "    value_col='repetition_effect',\n",
    "    condition_order=['Immediate', 'Delayed'],\n",
    "    line_colors=['#ff3d00', '#1f77b4'],\n",
    "    ylabel='Initial-Repeated (beta estimate)',\n",
    "    xlabel='Lag of repetition',\n",
    "    title=f\"{ROI_label} \\nRepetition effect: Face Familiarity x Repetition Lag\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Familiar vs Unfamiliar Faces in Amygdala"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data_from_amygdala = extract_mean_ROI_values(effect_files, amygdala_mask)\n",
    "print(data_from_amygdala)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But repetition effects can affect these results. To avoid that and only test the face familiarity effect, let's only use the initial presentations from the two conditions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "famous = data_from_amygdala[data_from_amygdala['condition'] == 'IniFF'].groupby('subject')['value'].mean()\n",
    "unfamiliar = data_from_amygdala[data_from_amygdala['condition'] == 'IniUF'].groupby('subject')['value'].mean()\n",
    "\n",
    "# Perform the paired t-test\n",
    "t_stat, p_value = stats.ttest_rel(famous, unfamiliar)\n",
    "print(f\"Paired t-test results: t-statistic = {t_stat}, p-value = {p_value}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_famous_unfamiliar = pd.DataFrame({\n",
    "    'Condition': ['Famous'] * len(famous) + ['Unfamiliar'] * len(unfamiliar),\n",
    "    'Value': pd.concat([famous, unfamiliar]),\n",
    "    'Subject': list(famous.index) + list(unfamiliar.index)\n",
    "})\n",
    "\n",
    "# Adjust the figure size to make it narrower\n",
    "plt.figure(figsize=(4, 4))  # Reduce the width (first number) to make the plot narrower\n",
    "\n",
    "# Plot the point plot with dodge to avoid overlap\n",
    "sns.violinplot(x='Condition', y='Value', data=data_famous_unfamiliar, linestyles='-', color='lightgrey')\n",
    "\n",
    "# Add lines\n",
    "for subject in data_famous_unfamiliar['Subject'].unique():\n",
    "    subject_data = data_famous_unfamiliar[data_famous_unfamiliar['Subject'] == subject]\n",
    "    plt.plot([0, 1], subject_data['Value'], color='gray', alpha=0.3)\n",
    "\n",
    "plt.ylabel('Beta estimate')\n",
    "plt.title('Paired t-test: Famous vs Unfamiliar')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Are the repetition effects indeed affecting the familiarity effect? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a data frame from amygdala_data for a 2-by-2 repeated measures ANOVA with factors 'familiarity' (Famous, Unfamiliar) and 'presentation' (Initial, Repeated). \n",
    "# Average the im and dl conditions for the repeated presentation.\n",
    "\n",
    "data_from_amygdala['familiarity'] = data_from_amygdala['condition'].apply(lambda x: 'Famous' if 'FF' in x else 'Unfamiliar')\n",
    "data_from_amygdala['presentation'] = data_from_amygdala['condition'].apply(lambda x: 'Initial' if 'Ini' in x else 'Repeated')\n",
    "\n",
    "# For repeated conditions, average the 'dl' and 'im' presentations\n",
    "def compute_average(df, condition_name):\n",
    "    return df[df['condition'].str.contains(condition_name)].groupby(['subject', 'familiarity', 'presentation'])['value'].mean().reset_index()\n",
    "\n",
    "# Get the initial presentations\n",
    "initial_df = data_from_amygdala[data_from_amygdala['presentation'] == 'Initial']\n",
    "\n",
    "# Get the repeated presentations by averaging dl and im\n",
    "repeated_df = compute_average(data_from_amygdala, 'Del|Imm')\n",
    "\n",
    "# Combine both initial and repeated into the final DataFrame for ANOVA\n",
    "data_amygdala_initial_repeated = pd.concat([initial_df, repeated_df], axis=0).sort_values(by=['subject', 'familiarity', 'presentation']).reset_index(drop=True)\n",
    "\n",
    "# Drop the now redundant 'condition' column\n",
    "data_amygdala_initial_repeated.drop(columns=['condition'], inplace=True)\n",
    "\n",
    "# Display the final DataFrame\n",
    "print(data_amygdala_initial_repeated)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform the repeated measures ANOVA\n",
    "aovrm = AnovaRM(data_amygdala_initial_repeated, depvar='value', subject='subject', within=['familiarity', 'presentation'])\n",
    "res = aovrm.fit()\n",
    "\n",
    "# Print the results\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_nbyn_anova(\n",
    "    data=data_amygdala_initial_repeated,\n",
    "    subject_col='subject',\n",
    "    condition_cols=['familiarity', 'presentation'],\n",
    "    value_col='value',\n",
    "    condition_order=['Initial', 'Repeated'],\n",
    "    line_colors=['#ff3d00', '#1f77b4'],\n",
    "    ylabel='Beta estimate',\n",
    "    xlabel='Face presentation',\n",
    "    title=\"Amygdala \\nFamiliarity x Presentation\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's explore it even further - for all three levels of presentation. This will be a 2-by-3 ANOVA. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a copy of the data\n",
    "data_amygdala_3levels = data_from_amygdala.copy()\n",
    "\n",
    "# Create new columns to identify 'familiarity' and 'presentation' levels\n",
    "data_amygdala_3levels['familiarity'] = data_amygdala_3levels['condition'].apply(lambda x: 'Famous' if 'FF' in x else 'Unfamiliar')\n",
    "\n",
    "def get_presentation_level(condition):\n",
    "    if 'Ini' in condition:\n",
    "        return 'Initial'\n",
    "    elif 'Del' in condition:\n",
    "        return 'Delayed'\n",
    "    elif 'Imm' in condition:\n",
    "        return 'Immediate'\n",
    "    return None\n",
    "\n",
    "data_amygdala_3levels['presentation'] = data_amygdala_3levels['condition'].apply(get_presentation_level)\n",
    "\n",
    "\n",
    "print(data_amygdala_3levels)  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_nbyn_anova(\n",
    "    data=data_amygdala_3levels,\n",
    "    subject_col='subject',\n",
    "    condition_cols=['familiarity', 'presentation'],\n",
    "    value_col='value',\n",
    "    condition_order=['Initial', 'Immediate', 'Delayed'],\n",
    "    line_colors=['#ff3d00', '#1f77b4'],\n",
    "    ylabel='Beta estimate',\n",
    "    xlabel='Face presentation',\n",
    "    title=\"Amygdala \\nFamiliarity x Presentation\",\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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}