{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Group-level analyses\n", "This week's lab is about group-level models, multiple comparison correction (MCC), and region-of-interest (ROI) analysis. In this notebook, we will focus on group-level models, which we'll demonstrate and explain using both FSL and Python. \n", "\n", "We'll focus on the \"summary statistics\" approach again, in which we'll demonstrate how we average $c\\beta$-terms across runs (in run-level analyses) and subjects (in grouplevel analyses) using the GLM. Then, we're going to show you how to test more extensive hypotheses in grouplevel models. \n", "\n", "**What you'll learn**: after this lab, you'll ...\n", "\n", "- understand the concept of the summary statistics approach\n", "- be able to construct different grouplevel models (in FSL)\n", "\n", "**Estimated time needed to complete**: 2 hours
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Some imports\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What are group-level models?\n", "Last week, we discussed the multilevel models and the summary statistics approach. Specifically, we focused on how data from different runs are usually (in non-fMRI contexts) are analyzed in a single multilevel GLM by \"concatenating\" the data. And how, in fMRI, we usually don't take this approach due to the computational burden and use the summary statistics approach, which analyzes each run separately and subsequently aggregates the data in a second, run-level GLM. \n", "\n", "In this notebook, we will extend this idea of analyzing data from multiple levels by looking at data from multiple subjects and how to analyze this data in group-level models. We will look at two \"flavors\" of group-level analyses: parametric and non-parametric.\n", "\n", "![](https://docs.google.com/drawings/d/e/2PACX-1vQxCH3WU3nTqFlHUZb49rf9zioivGQ-flVfRpwmXQx7OF5Wm_1T6gFMYQqpqt-NPITNHUaRoVYEREgT/pub?w=965&h=745)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Parametric analyses\n", "The most often used \"flavor\" of fMRI (group-level) analyses are *parametric*: it assumes that the data can be modeled using specific probability distributions. For example, we assume that the results of statistical tests of parameters (i.e., $t$-values), given that null hypothesis is true, are distributed according to the Students $t$-distribution (with a particular degrees-of-freedom):\n", "\n", "\\begin{align}\n", "t_{c\\hat{\\beta}} \\sim \\mathcal{T}(\\mathrm{df})\n", "\\end{align}\n", "\n", "where you can read the $\\sim$ symbol as \"is distributed as\". Importantly, the validity of the computed $p$-values depends on whether the choice of distribution is appropriate. If not, you might risk inflated type 1 or type 2 errors.\n", "\n", "The first-level and run-level GLMs that we have discussed so far are examples of parametric analyses. There are also *non-parametric* versions of the GLM that do not assume any particular form of distribution; while somewhat more computationally expensive, this is become a more and more popular alternative to (group-level) parametric analyses. Importantly, the difference between parametric and non-parametric analyses is only important for the *inference* (not the *estimation*) aspect of the (group-level) analyses." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's focus on the parametric version of group-level analyses. Basically, this amounts to doing the same thing as we did last week with the run-level analyses, but this time, the results from our run-level analyses ($c\\hat{\\beta}^{*}$) across different subjects will become our target ($y^{\\dagger}$). Note that we will use the \"dagger\" ($^{\\dagger}$) superscript to denote that the mathematical terms belong to the group-level model (just like the $^{*}$ superscript in our notebooks refers to terms belonging to the run-level models).\n", "\n", "To reiterate, the results from our run-level analyses ($c\\hat{\\beta}^{*}$), or first-level analyses if we only have a single run, become our dependent variable in our group-level analysis ($y^{\\dagger}$):\n", "\n", "\\begin{align}\n", "y^{\\dagger} = c\\hat{\\beta}^{*}\n", "\\end{align}\n", "\n", "Again, the group-level represents a GLM with a particular design matrix ($\\mathbf{X}^{\\dagger}$) and parameters ($\\beta^{\\dagger}$):\n", "\n", "\\begin{align}\n", "y^{\\dagger} = \\mathbf{X}^{\\dagger}\\beta^{\\dagger} + \\epsilon^{\\dagger}\n", "\\end{align}\n", "\n", "And the group-level parameters can be estimated with OLS:\n", "\n", "\\begin{align}\n", "\\hat{\\beta}^{\\dagger} = (\\mathbf{X}^{\\dagger\\ T} \\mathbf{X}^{\\dagger})^{-1}\\mathbf{X}^{\\dagger}y^{\\dagger}\n", "\\end{align}\n", "\n", "As mentioned last week, the parameter estimation procedure (i.e., estimating $\\hat{\\beta}^{\\dagger}$) is relatively straightforward, but the inference procedure depends on the specific variance approach: fixed, random, or mixed effects. If your aim is to perform inference to the population, *you should never use a fixed-effects type of analysis across subjects*, as this will typically inflate your type 1 error greatly.\n", "\n", "That leaves us with random effects and mixed effects GLMs. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Random effects\n", "Let's first focus on a random effects type of analysis (which is commonly used in group-level analysis in the SPM software package, for example), as this is relatively easy to understand. For now, let's use simulated data. Suppose that we want to do a group-analysis of our run-level flocBLOCKED data. Specifically, we are interested in the difference between faces and places ($H_{0}: \\beta_{face}^{*} = \\beta_{house}^{*} = 0$, $H_{a}: \\beta_{face}^{*} > \\beta_{house}^{*}$). As such, we'll use the corresponding $c\\hat{\\beta}^{*}$ terms from our run-level analysis (i.e., the contrasts against baseline, COPE1 and COPE2) as our new target $y^{\\dagger}$ as follows:\n", "\n", "\\begin{align}\n", "y^{\\dagger} = \\begin{bmatrix}\n", "c\\hat{\\beta}^{*}_{1, F>0} \\\\\n", "c\\hat{\\beta}^{*}_{2, F>0} \\\\\n", "\\vdots \\\\\n", "c\\hat{\\beta}^{*}_{N, F>0} \\\\\n", "c\\hat{\\beta}^{*}_{1, P>0} \\\\\n", "c\\hat{\\beta}^{*}_{2, P>0} \\\\\n", "\\vdots \\\\\n", "c\\hat{\\beta}^{*}_{N, P>0}\n", "\\end{bmatrix}\n", "\\end{align}\n", "\n", "For our simulation, we'll assume that we have 20 subjects ($N=20$)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToThink (0 points): You might think, why compute the difference contrast at the group-level, while we could have computed this already in the first-level analysis (and subsequently average this in the run-level and group-level analyses)? In fact, this is mathematically (largely) equivalent. However, most people prefer to postpone this step to the group-level analysis. Can you think of a reason why?\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we are going to simulate the data, let's construct our design-matrix, $\\mathbf{X}^{*}$. Assuming we have 20 subjects and two types of inputs ($c\\hat{\\beta}^{*}_{F>0}$ and $c\\hat{\\beta}^{*}_{P>0}$, our design matrix will be:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 20\n", "X = np.zeros((N * 2, 2))\n", "X[:N, 0] = 1\n", "X[N:, 1] = 1\n", "\n", "print(X)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToThink (1 point): In the above design matrix, we did not specify an intercept (i.e., a vector of ones). With good reason, because our group-level GLM will crash if we did. Explain shortly why an intercept for this particular design matrix (i.e., as represented by the variable X above) is not only redundant, but will make the GLM impossible to estimate. Hint: try adding an intercept, run a GLM, and see what happens.\n", "
" ] }, { "cell_type": "markdown", "metadata": { "deletable": false, "nbgrader": { "cell_type": "markdown", "checksum": "85964998ddee0bf6e5d68d9c37ec49c2", "grade": true, "grade_id": "cell-c6671f254bc384a5", "locked": false, "points": 1, "schema_version": 3, "solution": true, "task": false } }, "source": [ "YOUR ANSWER HERE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are going to simulate data for a single voxel, which shows a slight preference (in run-level $c\\hat{\\beta}^{*}$ estimates across subjects) for faces relative to places. Specifically, we'll assume that the voxel activates on average with $\\beta^{\\dagger}_{faces} = 5$ for faces, while it activates on average with $\\beta^{\\dagger}_{places} = 4.5$ for places:\n", "\n", "\\begin{align}\n", "y^{*} = \\mathbf{X}^{\\dagger}_{faces}\\cdot 5 + \\mathbf{X}^{\\dagger}_{places} \\cdot 4.5 + \\epsilon^{\\dagger} \n", "\\end{align}\n", "\n", "where the noise, $\\epsilon^{\\dagger}$, is normally distributed with mean 0 and a standard deviation of 1:\n", "\n", "\\begin{align}\n", "\\epsilon^{\\dagger} \\sim \\mathcal{N}(0, 1)\n", "\\end{align}\n", "\n", "(The symbol $\\sim$ stands for \"is distributed as\".) We'll implement this in code below:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 20\n", "beta_f = 5\n", "beta_p = 4.5\n", "std_noise = 1\n", "np.random.seed(42)\n", "noise = np.random.normal(0, std_noise, size=N * 2)\n", "\n", "y_sim = X @ np.array([beta_f, beta_p]) + noise" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And let's plot it:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(5, 16))\n", "plt.plot(y_sim[:N], np.arange(y_sim.size // 2), marker='o', ms=10)\n", "plt.plot(y_sim[N:], np.arange(y_sim.size // 2, y_sim.size), marker='o', ms=10)\n", "plt.ylim(y_sim.size, -1)\n", "plt.legend([r\"$c\\hat{\\beta}^{*}_{face}$\", r\"$c\\hat{\\beta}^{*}_{place}$\"], fontsize=15)\n", "plt.ylabel(\"Subjects\", fontsize=20)\n", "plt.xlabel(r\"$c\\hat{\\beta}^{*}$\", fontsize=20)\n", "plt.grid()\n", "plt.title(r\"$y^{\\dagger}$ (run-level estimates)\", fontsize=25)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that these $c\\hat{\\beta}^{*}$ estimates are the only thing you need for a random effects analysis! We're simply ignoring the lower-level variance terms." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (2 points): Run a group-level random-effects GLM (using the variable X as design matrix and y_sim as dependent variable). Store the estimated parameters in a variable named gl_params. Then, compute the $t$-value for the contrast $\\beta^{\\dagger}_{face} > \\beta^{\\dagger}_{place}$ and store this in a variable named gl_tvalue.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "72ecb7721bf542a4418de552de5ebf77", "grade": false, "grade_id": "cell-4e965a2198a792ed", "locked": false, "schema_version": 3, "solution": true, "task": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Implement the ToDo here. '''\n", "from numpy.linalg import inv\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "12a16c59d4aca87aeb589e5053c4cc24", "grade": true, "grade_id": "cell-c16ea54c4c4ebc3f", "locked": true, "points": 2, "schema_version": 3, "solution": false, "task": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Tests the above ToDo. '''\n", "from niedu.tests.nii.week_6 import test_rfx_glm\n", "test_rfx_glm(X, y_sim, gl_params, gl_tvalue)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Mixed-effects (in FSL)\n", "Proper mixed-effects models are a bit too complex to implement ourselves in Python, so we'll show you how to do it in FSL! We ran the (fixed-effects) run-level models for the flocBLOCKED data for you already, which we download below:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "data_dir = os.path.join(os.path.expanduser(\"~\"), 'NI-edu-data')\n", "\n", "print(\"Downloading run-level FSL FEAT results (+- 450MB) ...\")\n", "!aws s3 sync --no-sign-request s3://openneuro.org/ds003965 {data_dir} --exclude \"*\" --include \"derivatives/fsl/sub-*/flocBLOCKED/runlevel.gfeat/*\"\n", "print(\"\\nDone!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Within the flocBLOCKED directory, the first-level (ses-1.feat and ses-2.feat) and run-level results (runlevel.gfeat) are stored for each subject." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): Start a remote desktop environment, open a terminal, navigate to the runlevel.gfeat directory of the task flocBLOCKED of subject sub-03 (so sub-03/flocBLOCKED/runlevel.gfeat), and print the directory's contents to the terminal.\n", "\n", "(Note that the following ToDos will be checked later using test-cells.)\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You should see that the runlevel.gfeat directory contains multiple cope*.feat directories. These directories correspond to the results of the run-level (intercept-only) analyses of the seven different first-level contrasts. In contrast to what you did yourself last week (in which you concatenated the COPEs from the $c\\hat{\\beta}_{\\mathrm{face}}$ and $c\\hat{\\beta}_{\\mathrm{place}}$ contrasts in a single GLM), we were lazy and used the \"Inputs are lower-level FEAT directories\" option in FEAT. This runs the exact same GLM (using an intercept-only run-level design matrix) for all first-level contrasts separately. \n", "\n", "For example, the data in the runlevel.gfeat/cope1.feat directory corresponds to the first contrast from the first-level analysis ($\\beta_{face} > 0$), the data in the runlevel.gfeat/cope2.feat directory corresponds to the second contrast from the first-level analysis ($\\beta_{place} > 0$), etc. etc.\n", "\n", "In our intercept-only run-level analysis, we only computed a run-level contrast reflecting the average across the two sessions, i.e., using the contrast vector $$. The results of this run-level contrast, for each first-level contrast, are stored in the cope1.nii.gz file." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): Open Feat, change \"First-level analysis\" to \"Higher-level analysis\" and change the input type to \"Inputs are 3D cope images from feat directories\". Then, click on the \"Stats\" tab.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Under the \"Stats\" tab, you can change the type of inference you want FEAT to do. By default, it uses a particular mixed-effects approach that FSL calls \"FLAME 1\". In general, we recommend using this type of inference. The other mixed-effects approach, \"FLAME 1+2\", is also a good choice, but may take much longer to run. The \"Mixed effects: Simple OLS\" is, unlike the same suggests, the random-effects option, which discards the lower-level variance estimates. Lastly, the \"Randomise\" option is a non-parametric (random-effects) variant, that will be discussed in more detail later.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): Set the option to \"Mixed effects: FLAME 1\". Enable the \"use automatic outlier de-weighting\" option (should be yellow). Then, go to post-stats and set the \"Thresholding\" option to \"Cluster\" with a Z threshold of 3.7 and a \"Cluster P threshold\" of 0.05. This refers to a particular type of multiple comparison correction that will be discussed later.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's define the dependent variable ($y^{\\dagger}$). Suppose, just like we did with the simulated data, that we have the following null and alternative hypothesis:\n", "\n", "$H_{0}: \\beta_{face}^{*} - \\beta_{house}^{*} = 0$
\n", "$H_{a}: \\beta_{face}^{*} - \\beta_{house}^{*} > 0$\n", "\n", "Note that the contrast evaluating the difference between these two parameters was *not* computed in the first-level or run-level analyses, so we have to create this ourselves by concatenating the COPEs corresponding to the $\\beta_{face}^{*} > 0$ and $\\beta_{house}^{*} > 0$ run-level contrasts, just like you did last week. \n", "\n", "Doing this for all subjects becomes a bit cumbersome though, so you only need to do this for sub-02, sub-03, and sub-04. Make sure the first three inputs correspond to the face-against-baseline COPE (COPE nr. 1) and the last three inputs correspond to the place-against-baseline COPE (COPE nr. 2). \n", "\n", "\\begin{align}\n", "y^{\\dagger} = \\begin{bmatrix}\n", "c\\hat{\\beta}^{*}_{\\mathrm{face>0,\\ sub-02}} \\\\\n", "c\\hat{\\beta}^{*}_{\\mathrm{face>0,\\ sub-03}} \\\\\n", "c\\hat{\\beta}^{*}_{\\mathrm{face>0,\\ sub-04}} \\\\\n", "c\\hat{\\beta}^{*}_{\\mathrm{place>0,\\ sub-02}} \\\\\n", "c\\hat{\\beta}^{*}_{\\mathrm{place>0,\\ sub-03}} \\\\\n", "c\\hat{\\beta}^{*}_{\\mathrm{place>0,\\ sub-04}} \\\\\n", "\\end{bmatrix}\n", "\\end{align}\n", "\n", "Note that you need the actual cope*.nii.gz files — *not* the .feat directory itself.\n", "\n", "(Realize, though, that having only 3 subjects for your group-analysis would be a hopelessly underpowered study.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (1 point): In the \"Data\" tab, use the option \"Inputs are 3D cope images from FEAT directories\", set the number of inputs to 6, and select the correct files. Note that within the \"Select input data\" window, you can copy/paste the paths using CTRL+C / CTRL+V, which might make the process a little easier/faster (not forget to actually modify the paths after copy/pasting).\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): In the \"Data\" tab, set the output directory to your home folder + grouplevel_task-flockBLOCKED.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (1 point): In the “Stats” tab, click on \"Full model setup\". Create two predictors (\"EVs\") that model the mean of the $c\\beta^{*}_{face} > 0$ (first EV) and the $c\\beta^{*}_{place} > 0$ (second EV). Then, add another predictor that models the effect of age on the $c\\beta^{*}_{face}$ contrast *only*. The ages of the three participants (from sub-02 to sub-04) are: $[21, 32, 25]$. Importantly, make sure you center (de-mean) the age values before entering them in FEAT, as this improves the interpretation of the other predictors (you can compute the de-meaned age values, for example, in a new code cell).\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (1 point): In the \"Contrasts & F-tests\" tab, define the following group-level contrasts:
\n", " \n", "1. $H_{0}: \\beta^{\\dagger}_{face} = 0$
$H_{a}: \\beta^{\\dagger}_{face} > 0$

\n", "2. $H_{0}: \\beta^{\\dagger}_{place} = 0$
$H_{a}: \\beta^{\\dagger}_{place} > 0$

\n", "3. $H_{0}: \\beta^{\\dagger}_{face} - \\beta^{\\dagger}_{place} = 0$
$H_{a}: \\beta^{\\dagger}_{face} - \\beta^{\\dagger}_{place} > 0$

\n", "4. $H_{0}: \\beta^{\\dagger}_{age} = 0$
$H_{a}: \\beta^{\\dagger}_{age} > 0$

\n", "5. $H_{0}: \\beta^{\\dagger}_{age} = 0$
$H_{a}: \\beta^{\\dagger}_{age} < 0$
\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): Save your design (using the \"Save\" button in the bottom of the FEAT window) in your week_6 directory and name it setup_feat_grouplevel.fsf.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below, you can test your implementation." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "66d58021c9b24d32ebf9a852e4a475e0", "grade": true, "grade_id": "cell-10a7400fbd55ef9a", "locked": true, "points": 1, "schema_version": 3, "solution": false, "task": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Tests your setup_feat.fsf file: inputs'''\n", "import os\n", "import os.path as op\n", " \n", "par_dir = op.basename(op.dirname(op.abspath('.')))\n", "if par_dir != 'solutions': # must be a student-account\n", " fsf = 'setup_feat_grouplevel.fsf'\n", " if not op.isfile(fsf):\n", " print(\"Couldn't find a 'setup_feat_group.fsf' file in your week_6 directory!\")\n", "\n", " with open(fsf, 'r') as f:\n", " fsf = f.readlines()\n", "\n", " fsf = [txt.replace('\\n', '') for txt in fsf if txt != '\\n']\n", " fsf = [txt for txt in fsf if txt != '#'] # remove commnts\n", "\n", " from niedu.tests.nii.week_6 import test_grouplevel_fsf_inputs\n", " test_grouplevel_fsf_inputs(fsf)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "2dae55e8984cb0307188fd92fff990a2", "grade": true, "grade_id": "cell-901625a5247a5817", "locked": true, "points": 1, "schema_version": 3, "solution": false, "task": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Tests your setup_feat.fsf file: Evs'''\n", "if par_dir != 'solutions':\n", " from niedu.tests.nii.week_6 import test_grouplevel_fsf_evs\n", " test_grouplevel_fsf_evs(fsf)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "ae1162cd5f462c2059323e832cc3fedc", "grade": true, "grade_id": "cell-6ce964d3df3cfae3", "locked": true, "points": 1, "schema_version": 3, "solution": false, "task": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Tests your setup_feat.fsf file: contrasts'''\n", "if par_dir != 'solutions':\n", " from niedu.tests.nii.week_6 import test_grouplevel_fsf_cons\n", " test_grouplevel_fsf_cons(fsf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Non-parametric analyses\n", "*This is an optional section*\n", "\n", "Before we go on with the topic of \"multiple comparison correction\" in the next notebook, let's focus on another type of inference in statistical models: permutation-based non-parametric analyses. Unlike parametric analyses (which we discussed in the previous section), non-parametric analyses do not assume anything about how our statistics-of-interest are distributed. Instead, it uses the data itself to estimate a distribution of your test-statistic under the null hypothesis. Note that this type of inference is also available in FSL, where it is called \"randomise\". \n", "\n", "Permutation tests work (roughly) as follows:\n", "1. The statistic (e.g., $t$-value) is computed as normal;\n", "2. The rows of design matrix ($\\mathbf{X}$) are shuffled and the statistic is (re-)computed, which is iterated a number of times (e.g., 1000);\n", "3. The (\"non-parametric\") $p$-value is computed\n", "\n", "Let's do this for our previously simulated data (the variables X and y_sim)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Computing the \"observed\" statistic\n", "First, we need to compute our \"observed\" statistic. In the framework of the GLM, this is often a $t$-statistic (or $F$-statistic). But note that permutation-based tests work for any kind of statistics. (In fact, permutation-based tests are especially useful for statistics with unknown distributions!)\n", "\n", "We are sure that you, by now, know how to compute $t$-values, so we are going to use a predefined function compute_tvalue with three arguments: X, y, and c. We'll compute the $t$-value for our data (X and y_sim) and store it in a variable called t_obs (for \"observed\")." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from niedu.utils.nii import compute_tvalue\n", "t_obs = compute_tvalue(X, y_sim, c=np.array([1, -1]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Run the permutations\n", "The second step is to, iteratively, permute the data and (re-)compute the test statistic. In most cases (assuming that the datapoints, $y$, are completely independent, i.e., there is no \"autocorrelation\"), permutation within the GLM framework is done by shuffling the rows of the design matrix ($\\mathbf{X}$). This is usually done a larger number of times (e.g., 1000 or more)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): Below, we started implementing the permutation loop. Can you finish it? You may want to use the np.random.permutation function. Make sure to save the permuted t-value each iteration (by storing it in the t_perm array).\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "3e8396f24ed1ce2a6a59932bd72c52d9", "grade": false, "grade_id": "cell-45242b07a8272101", "locked": false, "schema_version": 3, "solution": true, "task": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "\"\"\" Implement the ToDo here. \"\"\"\n", "np.random.seed(42)\n", "iters = 1000\n", "t_perm = np.zeros(iters)\n", "\n", "for i in range(iters):\n", " # YOUR CODE HERE\n", " raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "021ac6eac22ef2623aad94784dabfb0b", "grade": true, "grade_id": "cell-ff7aa4c2c4c1c543", "locked": true, "points": 0, "schema_version": 3, "solution": false, "task": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Tests the ToDo above. '''\n", "np.testing.assert_almost_equal(t_perm.mean(), 0, decimal=1)\n", "print(\"Well done!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Calculate the p-value\n", "The last step is to calculate the $p$-value. This is done by dividing the number of permuted statistics *larger* than the observed statistic plus one by the total number of permutations plus one. Or, mathematically, for $P$ permutations, the $p$-value for any statistic $s$ with observed value $s^{\\mathrm{obs}}$ and vector of permuted values $\\mathbf{s}^{\\mathrm{perm}}$ is defined as:\n", "\n", "\\begin{align}\n", "p = \\frac{\\sum_{i=1}^{P}\\mathbf{I}(\\mathbf{s}^{\\mathrm{perm}}_{i} \\geq s^{\\mathrm{obs}}) + 1}{P + 1}\n", "\\end{align}\n", "\n", "where $\\mathbf{I}$ is an indicator function, returning $1$ when the condition is true and $0$ otherwise.\n", "\n", "We can also visualize it:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(10, 5))\n", "plt.hist(t_perm)\n", "plt.axvline(t_obs, ls='--', c='red', lw=3)\n", "plt.legend(['observed', 'permuted'])\n", "plt.xlabel(r'$t$-value', fontsize=20)\n", "plt.ylabel('Frequency', fontsize=20)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the above figure, the number of values right of the dashed red line ($t^{obs}$) plus 1 divided by the number of permutations plus 1 ($P+1$) is our $p$-value! Now, try computing it yourself!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): Compute the $p$-value corresponding to our observed $t$-statistic (t_obs) and store it in a variable named p_perm.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "1e42aa3b4d6ac3ee8dc86ded23bde7d5", "grade": false, "grade_id": "cell-fe849a0f7de10faa", "locked": false, "schema_version": 3, "solution": true, "task": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "# Implement your ToDo here\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "2bf14780b29b0718f6c9914237835c01", "grade": true, "grade_id": "cell-8accc117937c5d64", "locked": true, "points": 0, "schema_version": 3, "solution": false, "task": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Tests the above ToDo. \n", "Note: this may differ slightly due to the exact implementation of the permutation procedure\n", "'''\n", "np.testing.assert_almost_equal(p_perm, 0.03, decimal=2)\n", "print(\"Well done!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): Shuffling rows works for most designs, but what if you have an intercept-only model (i.e., a one-sample $t$-test), where a design matrix consisting only of a vector of ones? Permuting the rows won't work there. An alternative for this scenario is \"sign flipping\": randomly changing ones to minus ones (with replacement).\n", "\n", "Below, we simulate some new data (y_sim2 and X2). Can you compute the non-parametric $p$-value for this one-sample $t$-test? Use a 1000 permutations and store the $p$-value in a new variable named p_perm2.\n", "\n", "Hint: to generate random \"signs\", take a look here.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "37fc77530f81857541ff047dafba31d0", "grade": false, "grade_id": "cell-1bf6a2d0c9a7d3d5", "locked": false, "schema_version": 3, "solution": true, "task": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "np.random.seed(17)\n", "y_sim2 = np.random.normal(0, 1, size=100)\n", "X2 = np.ones(100)[:, np.newaxis]\n", "\n", "iters = 1000\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "05565fb0ff874727fc82e82d0114a940", "grade": true, "grade_id": "cell-a55f54ec285a2d03", "locked": true, "points": 0, "schema_version": 3, "solution": false, "task": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Tests the above ToDo. '''\n", "# The p-value might slightly differ from 0.17 depending on the\n", "# exact method of generating random signs.\n", "np.testing.assert_almost_equal(p_perm2, 0.17, decimal=2)\n", "print(\"Well done!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Threshold-Free Cluster Enhancement (TFCE)\n", "Apart from (or perhaps because of) the fact that non-parametric analyses are blind to the distribution of effects, non-parametric analyses have an additional advantage: they can use a \"trick\" called \"Threshold-Free Cluster Enhancement\" (TFCE). This method allows for voxel-based inference while taking into account the spatial size of the effect. You could see TFCE as a method to \"boost\" high-amplitude voxels that are located within a cluster of other high-amplitude voxels. \n", "\n", "However, we don't know how these \"boosted\" values are distributed (they might or might not be $\\mathcal{T}$ distributed)! Fortunately, as said before, non-parameteric analyses don't care how the computed statistic is distributed (under the null). How TFCE works mathematically is beyond the scope of the course (but you can learn more in [this video](https://www.youtube.com/watch?v=q7cWw8WC0Ws)), but in fact recommend using this technique in combination with FSL's randomise tool, also because it provides an elegant and effective way to control for multiple comparisons (the topic of section 3)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Checking out results from group-level FEAT analyses\n", "We have run a group-level model for a couple of our lower-level contrasts. The results from the group-level analysis that we are going to take a look at will be downloaded below:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"Downloading group-level FSL FEAT results (+- 67MB) ...\")\n", "!aws s3 sync --no-sign-request s3://openneuro.org/ds003965 {data_dir} --exclude \"*\" --include \"derivatives/fsl/grouplevel_task-flocBLOCKED/contrast-faceGTother_method-FLAME1_thresh-uncorr05.gfeat/*\"\n", "print(\"\\nDone!\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "gl_dir = os.path.join(data_dir, 'derivatives', 'fsl', 'grouplevel_task-flocBLOCKED', 'contrast-faceGTother_method-FLAME1_thresh-uncorr05.gfeat')\n", "print(f\"Data is located here: {gl_dir}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This contains the results from the group-level analysis of the lower-level $4 \\cdot \\beta_{face} - \\beta_{place} - \\beta_{body} - \\beta_{character} - \\beta_{object} > 0$ contrast. (The \"faceGTother\" stands for \"face greater than other conditions\".)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): Open a terminal, navigate to the above directory, and print its contents to the terminal window.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You should see the usual files: HTML-reports, thresholded results (as nifti files), and files with information about the design (design.con, design.fsf, etc.). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): Open the report.html file (with the command firefox report.html). Click on \"Inputs\", which will show you the data that we used as inputs (i.e., $y^{\\dagger}$) for our group-level analysis.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, the inputs for our group-level analysis are not individual cope*.nii.gz files, but entire .feat directories. This is because we used the \"Inputs are lower-level FEAT directories\" option." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): Click on \"Results\". This will show you a page with the simple \"intercept-only\" design (i.e., a design for a one-sample $t$-test). Now, click on \"Lower-level contrast 1 (average)\" to view in the group-level result from this particular run-level contrast.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On this page, you can view the results, as $z$-statistic images, from the group-level analysis plotted on a set of axial slices. Below the brain plot, you can see the \"time series plot\", which does not actually show a \"time series\", but rather the input (i.e., $y^{\\dagger}$) of the voxel with the highest $z$-value. Thus, the \"time series\" here do not reflect datapoints over time (as in first-level analyses), but run-level $c\\hat{\\beta}^{*}$ estimates of different subjects!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): In your terminal, open FSLeyes. Then, click on FileAdd standard → select MNI152_T1_2mm_brain.nii.gz → click on Open.\n", "\n", "Then, click on FileAdd from directory → navigate to and select the contrast-faceGTothers_method-FLAME1_thresh-uncorr05.gfeat/cope1.feat directory → click Open.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This automatically opens the filtered_func_data.nii.gz file: a 4D nifti file with the run-level $c\\hat{\\beta}^{*}$ estimates of different subjects across the fourth dimension (i.e., data from different subjects are represented by different volumes, our new $y^{\\dagger}$). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): Set the minimum value to 100 (in the top menu) and change the colormap to \"Red-Yellow\". Note that the minimum value of 100 is rather arbitrary, but it removes some of the \"visual clutter\".\n", "\n", "Now, go to voxel $[26, 39, 25]$, which is (roughly) located in the right fusiform gyrus (the \"FFA\") and, for this participant (sub-02), seems to activate relatively strongly in response to faces (rather than to other images). \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): Increase the volume number from 0 to 1 (in the \"Location\" panel), which will show you the run-level $c\\hat{\\beta}^{*}$ estimates from a different subject (i.e., sub-03). You'll notice that this participant's \"FFA\" is in a slightly different location! In fact, when you view the other volumes (i.e., the $c\\hat{\\beta}^{*}$ maps of other subjects), you see that the anatomical location of the effects varies substantially!\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToThink (1 point): Name two reasons why an effect (i.e., a significant \"blob\") may differ in anatomical location across subjects.\n", "
" ] }, { "cell_type": "markdown", "metadata": { "deletable": false, "nbgrader": { "cell_type": "markdown", "checksum": "c0befc4e9d39c1737380da64d662cb92", "grade": true, "grade_id": "cell-f019c8631639c33a", "locked": false, "points": 1, "schema_version": 3, "solution": true, "task": false } }, "source": [ "YOUR ANSWER HERE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alright, now let's look at some real group-level results in FSLeyes!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (ungraded): In FSLeyes, delete the cope1: filtered_func_data overlay (or make it invisible by clicking on the eye icon). Then, open the thresh_zstat1.nii.gz (FileAdd from file) and set the colormap to \"Red-Yellow\".\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, you see the \"thresholded\" $z$-values plotted on the template (MNI152, 2mm) that FSL by default uses for its group-analyses. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo/ToThink (1 point): Go to the voxel at location $[37, 60, 28]$. According to the Harvard-Oxford Subcortical structural atlas, to what anatomical brain region does this voxel most likely belong? Write this down in the text cell below. Check out the notebook from last week if you forgot how to open the Atlas panel.\n", "
" ] }, { "cell_type": "markdown", "metadata": { "deletable": false, "nbgrader": { "cell_type": "markdown", "checksum": "ddfecfeb5d732aaea197cdb85b546bba", "grade": true, "grade_id": "cell-4b2caba8aa83274f", "locked": false, "points": 1, "schema_version": 3, "solution": true, "task": false } }, "source": [ "YOUR ANSWER HERE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The red-yellow voxels that you see are those that have $z$-values corresponding to $p$-values smaller than 0.05 (the other voxels with $p$ > 0.05 are set to 0). This particular thresholding procedure gives what people often refer to as \"uncorrected results\". Although they are thresholded at a particular $p$-value, the results are referred to as \"uncorrected\" because they have not been corrected for multiple comparisons — the topic of the next notebook (MCC.ipynb)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Tip!\n", " Before handing in your notebooks, we recommend restarting your kernel (KernelRestart & Clear Ouput) and running all your cells again (manually, or by CellRun all). By running all your cells one by one (from \"top\" to \"bottom\" of the notebook), you may spot potential errors that are caused by accidentally overwriting your variables or running your cells out of order (e.g., defining the variable 'x' in cell 28 which you then use in cell 15).\n", "
" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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" }, "toc": { "base_numbering": 1, "nav_menu": { "height": "296px", "width": "389px" }, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }