{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Temporal preprocessing\n", "Preprocessing of (f)MRI data is quite a complex topic, involving techniques from signal processing (filtering), linear algebra/statistics (prewhitening, autocorrelation correction), and numerical optimization (image registration). Additionally, there is ongoing discussion about which preprocessing steps are necessary/sufficient and what preprocessing parameters are optimal. As always, this depends on your specific research question, the type and quality of your data, and the intended analysis.\n", "\n", "In this week, we'll discuss a couple (not all!) of preprocessing steps that are common in univariate fMRI analyses across two notebooks, one discussing temporal processing (the current notebook) and one discussing spatial preprocessing (the next notebook, `spatial_preprocessing.ipynb`).\n", "\n", "**What you'll learn**: after this lab, you'll ...\n", "\n", "- be able to explain the influence of preprocessing on the measured effects using the t-value formula\n", "- understand (the advantage of) temporal filtering from both the time-domain and frequency-domain\n", "- understand the necessity of prewhitening given the assumptions of the GLM\n", "- understand the advantage of spatial filtering (smoothing)\n", "- understand how to handle outliers\n", "- know how to implement the concepts above in Python\n", "\n", "**Estimated time needed to complete**: 8-12 hours" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "As we said before, preprocessing is a topic that almost warrants its own course. Nonetheless, we'll try to show you (and let you practice with) some of the most common and important preprocessing operations. Additionally, we'll introduce the concept of the fast fourier transform, which allows us to analyze our signal in the frequency domain, which helps to understand several preprocessing steps, such as temporal filtering." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The *t*-value formula — yet again\n", "The previous two weeks, you have learned that, essentially, we want to find large effects (calculated as *t*-values) of our contrasts by optimizing various parts of the *t*-value formula. Conceptually, the *t*-value formula for a particular contrast ($\\mathbf{c}$) can be written as:\n", "\n", "\\begin{align}\n", "t\\mathrm{-value} = \\frac{\\mathrm{effect}}{\\mathrm{uncertainty}} = \\frac{\\mathrm{effect}}{\\sqrt{\\mathrm{noise} \\cdot \\mathrm{design\\ variance}}} = \\frac{\\mathbf{c}\\hat{\\beta}}{\\sqrt{\\hat{\\sigma}^{2}\\mathbf{c}(\\mathbf{X}^{T}\\mathbf{X})^{-1}\\mathbf{c}^{T}}}\n", "\\end{align}\n", "\n", "Last week you've learned that by ensuring low design variance (through *high* predictor variance and *low* predictor correlations) leads to larger normalized effects (higher *t*-values). This week, we will discuss the other term of the denominator of this formula: the noise (also called the residual variance or unexplained variance), which is defined in the *t*-value formula as follows:\n", "\n", "\\begin{align}\n", "\\mathrm{noise} = \\frac{\\mathrm{SSE}}{\\mathrm{DF}} = \\frac{\\sum_{i=1}^{N}(y_{i} - \\hat{y}_{i})^{2}}{N - P}\n", "\\end{align}\n", "\n", "Through preprocessing, we aim to reduce the difference between our prediction ($\\hat{\\mathbf{y}}$) and our true signal ($\\mathbf{y}$), thus reducing the noise-term of the formula and thereby optimizing our (\"normalized\") effects.\n", "\n", "### On temporal vs. spatial preprocessing\n", "Roughly speaking, you can divide fMRI preprocessing into two types of operations: temporal and spatial. Temporal preprocessing involves operations that filter or otherwise affect the properties of your data across the time dimension (i.e., the time series). This is, of course, only applicable to *functional* MRI data (not structural or diffusion MRI). Examples of temporal operations are high-pass filtering and slice-time correction. Spatial preprocessing involves operations that filter or otherwise affect the spatial properties of your data (such as spatial orienatation, resolution, and shape). Examples of these spatial operations are spatial filtering (\"smoothing\"), distortion correction, motion correction (realignment), and spatial normalization/resampling.\n", "\n", "In this lab, we'll discuss these operations per theme: we'll start with temporal operations and then move to spatial operations. Note that this *is not* the order in which data is usually preprocessed (see page 35 of the \"Handbook of functional MRI Data Anslysis\" textbook).\n", "\n", "Let's start with temporal preprocessing." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# First, some imports!\n", "import numpy as np\n", "from numpy.linalg import inv\n", "from nilearn.glm.first_level.hemodynamic_models import glover_hrf\n", "from scipy.interpolate import interp1d\n", "\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Temporal preprocessing\n", "In this section, we will discuss how temporal preprocessing may greatly reduce the error term of our statistics. We'll also delve into a variant of the GLM that deals with autocorrelation in our signals appropriately (\"generalized least squares\", or GLS)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Slice-time correction\n", "Slice-time correction is a temporal resampling technique that corrects for the fact that, in (most) BOLD-MRI scan sequences, volumes are acquired slice by slice. This means that each slices is acquired at a slightly different time. For example, suppose that we are acquiring an fMRI run with a TR of 2 seconds and our volumes consist of 40 slices (acquired axially, inferior → superior). In this case, the acquisition of each slice takes $\\frac{TR}{N_{slice}} = \\frac{2}{40} = 0.05$ seconds (let's call this `dt`). This means that, for the very first volume, the \"onset\" of the first slice is at 0 seconds (and lasts until 0.05 seconds), while the onset of the last slice is at 1.95 (and lasts until 2.00 seconds). In code, we can calculate these slice onsets **within a volume** using the `np.linspace(start, stop, n_steps)` function:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "TR = 2\n", "n_slices = 40\n", "slice_onsets_within_volume = np.linspace(0, TR, n_slices, endpoint=False)\n", "# Note, this would have worked as well: np.arange(0, TR, TR / 40)\n", "\n", "print(\"Length of slice onsets: %i\" % slice_onsets_within_volume.size, end='\\n\\n')\n", "print(\"Onsets of slices for first volume: %s\" % (slice_onsets_within_volume,), end='\\n\\n')\n", "print(\"Acquisition of first slice (of first volume) started at %.2f sec.\" % slice_onsets_within_volume[0])\n", "print(\"Acquisition of last slice (of first volume) started at %.2f sec.\" % slice_onsets_within_volume[-1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's assume we do a very simple experiment, lasting 32 seconds (i.e., 16 volumes with a TR of 2), in which we show a single stimulus at $t=0$. Furthermore, let's assume that the *entire brain* responds to this stimulus with an idealized (noiseless response). In other words, the response of each voxel will look exactly like the HRF:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "oversampling = 100\n", "exp_length = 32\n", "n_vols = 32 // TR\n", "\n", "ideal_response = glover_hrf(tr=TR, time_length=exp_length, oversampling=TR * oversampling)\n", "ideal_response /= ideal_response.max()\n", "t = np.linspace(0, exp_length, ideal_response.size, endpoint=False)\n", "\n", "plt.figure(figsize=(15, 4))\n", "plt.plot(t, ideal_response)\n", "plt.grid()\n", "plt.title(\"Idealized response of all voxels\", fontsize=25)\n", "plt.xlabel(\"Time (seconds)\", fontsize=20)\n", "plt.ylabel(\"Activation (A.U.)\", fontsize=20)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We must, however, take into accout that, even when our data would represent an idealized response, it would be sampled at different timepoints due to a difference in onsets for different slices. The first slice would for example be sampled at $[0, 2, 4, 6, ..., 30]$ seconds, while, for example, the middle slice would be sampled at $[1, 3, 5, ... 31]$ seconds. We can calculate the onset per slice **across volumes** using `np.linspace` again:\n", "\n", "```python\n", "# Note: this is the onset of slice number \"slice_num\" **across volumes**\n", "onsets_slice_across_volumes = np.linspace(dt * slice_num, length_exp + dt * slice_num, n_vols, endpoint=False)\n", "```\n", "\n", "where `dt` is the time it takes per slice (i.e., $dt = \\frac{TR}{N_{slice}}$). So, to get the slice times for, e.g. slice number 23, for the experiment we have so far (i.e., with a TR of 2, 40 slices, and 32 second duration), we would compute the onsets as follows:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "dt = TR / n_slices\n", "slice_num = 23\n", "\n", "# note the slice_num - 1, as Python uses 0-based indexing\n", "start = dt * (slice_num - 1)\n", "end = exp_length + dt * (slice_num - 1)\n", "onsets_slice23_across_volumes = np.linspace(start, end, n_vols, endpoint=False)\n", "print(onsets_slice23_across_volumes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, even with a noiseless response, the data will look different per slice. We'll plot this below for the first slice and the middle slice (on top of the idealized response):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "onsets_s0_across_volumes = np.linspace(0, exp_length, n_vols, endpoint=False)\n", "onsets_s20_across_volumes = np.linspace(dt * 19, exp_length + dt * 19, n_vols, endpoint=False)\n", "\n", "print(\"Onset acquisition of slice 0 (first 10 volumes): %s\" % (onsets_s0_across_volumes[:10],))\n", "print(\"Onset acquisition of slice 20 (first 10 volumes): %s\" % (onsets_s20_across_volumes[:10],))\n", "\n", "# \"Fit\" the interpolation \n", "resampler = interp1d(t, ideal_response)\n", "\n", "plt.figure(figsize=(15, 4))\n", "plt.plot(t, ideal_response)\n", "plt.grid()\n", "plt.title(\"Resampled slicewise responses\", fontsize=25)\n", "plt.xlabel(\"Time (seconds)\", fontsize=20)\n", "plt.ylabel(\"Activation (A.U.)\", fontsize=20)\n", "\n", "for onsets in [onsets_s0_across_volumes, onsets_s20_across_volumes]:\n", " # Do not use this code as an example to solve the next ToDo!\n", " # Because it does it the other way around.\n", " slice_sig = resampler(onsets) # resample to slice times! kind of like inverse slice-time correction\n", " plt.plot(onsets, slice_sig, marker='*', ls='None', ms=10)\n", "\n", "plt.legend(['Idealized response', 'Actual response slice 0', 'Actual response slice 20'],\n", " fontsize=15)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, if we would ignore the fact that these signals were acquired at different times (by assuming they were all acquired at $t=0, t=2, t=4$ etc.), our model will be suboptimal for all slices (except for slice 0). We'll show this below by plotting the resampled data from slice 0, 20, and 39 on top of the idealized response (which represents our predictor):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(15, 4))\n", "\n", "slice_data = np.zeros((slice_sig.size, 3))\n", "slice_onsets_across_volumes = np.zeros((slice_sig.size, 3))\n", "for i, slc in enumerate([0, 19, 39]):\n", " t_slice = np.linspace(dt * slc, exp_length + dt * slc, n_vols, endpoint=False)\n", " slice_resp = resampler(t_slice)\n", " slice_data[:, i] = slice_resp # save for later\n", " slice_onsets_across_volumes[:, i] = t_slice\n", " plt.subplot(1, 3, i+1)\n", " plt.plot(t, ideal_response)\n", " plt.plot(onsets_s0_across_volumes, slice_resp)\n", " plt.title(\"Slice %i\" % (slc + 1), fontsize=25)\n", " plt.grid()\n", " \n", " if i == 0:\n", " plt.ylabel(\"Activation (A.U.)\", fontsize=20)\n", " \n", " if i == 1:\n", " plt.xlabel(\"Time (seconds)\", fontsize=20)\n", " plt.legend(['Model', 'Actual response'], fontsize=12)\n", "\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, \"later\" slices are quite misaligned relative to the idealized response/model (which assume that each slice is acquired at the start of the volume). We can fix this by *resampling the data from different slices to the onset times of a reference slice*. Basically, this amounts to saying: \"what would the signal from slice x look like if it was acquired at the onsets of a reference slice\". This process is called *slice-time correction*. \n", "\n", "Remember the temporal resampling of our \"high precision\" predictors that we did in week 2? We're going to do something similar to our data, here. Basically, we are going to resample our slice-wise timesieres corresponding to a particular set of onsets (let's call these the `original_onsets`) to the set of onsets of a particular reference slice (let's call these the `desired_onsets`). Using the same `interp1d` function from week 2, we can do the following:\n", "\n", "```python\n", "resampler = interp1d(original_onsets, original_signal)\n", "std_signal = resampler(desired_onsets)\n", "```\n", "\n", "where `original_onsets` is a numpy array with onsets of each slice **across volumes** and the `original_signal` refers to the actual signal of that voxel (here: slice). They should be the same size. The `desired_onsets` is another numpy array with, *in the case of slice-time correction*, the same size as `original_onsets`.\n", "\n", "Below, we'll pick the first slice as the reference slice." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(15, 8))\n", "\n", "plt.subplot(2, 1, 1)\n", "for i in range(slice_data.shape[1]):\n", " plt.plot(onsets_s0_across_volumes, slice_data[:, i])\n", " plt.grid()\n", " plt.ylabel(\"Activation (A.U.)\", fontsize=20)\n", " plt.title(\"Original slicewise responses\", fontsize=25)\n", "\n", "# Define reference slice as the first one (with onsets [0, 2, 4, 6, etc])\n", "desired_onsets = onsets_s0_across_volumes\n", "\n", "plt.subplot(2, 1, 2)\n", "for i in range(slice_data.shape[1]):\n", " \n", " # Resample slice-specific onsets to desired onsets\n", " resampler = interp1d(slice_onsets_across_volumes[:, i], slice_data[:, i], fill_value='extrapolate')\n", " stc_slice = resampler(desired_onsets)\n", " \n", " # Plot relative to the reference onsets, i.e., onsets_s0 (first slice)\n", " plt.plot(onsets_s0_across_volumes, stc_slice)\n", " plt.grid()\n", " plt.title(\"Slice-time corrected data\", fontsize=25)\n", " plt.xlabel(\"Time (seconds)\", fontsize=20)\n", " plt.ylabel(\"Activation (A.U.)\", fontsize=20)\n", " \n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, the resampled time series resemble each other much more after slice-time correction! It's showing some \"edge artifacts\" which are caused by extrapolation of the time series. In actual neuroimaging software packages, higher-order resampling (such as \"cubic\" or \"spline\" resampling) are used to somewhat mitigate this. Moreover, instead of the first slice, usually the middle slice is chosen as the reference slice, which reduces the amount of resampling (and extrapolation) that has to be done, further mitigating artifacts. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "