{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# The GLM, part 2: inference\n", "In this notebook, we'll continue with the GLM, focusing on statistical tests (i.e., inference) of parameters. Note that there are two notebooks this week: this one, glm_part2_inference.ipynb, and design_of_experiments.ipynb. Please do this one first.\n", "\n", "Last week, you learned how to estimate parameters of the GLM and how to interpret them. This week, we'll focus on statistical inference of those estimated parameters (and design of experiment, in another notebook). Importantly, we are going to introduce the most important formula in the context of univariate fMRI analyses: the formula for the *t-value*. Make sure you understand this formula, as we will continue to discuss it in the next weeks.\n", "\n", "**What you'll learn**: after this week's lab ... \n", "* you know how the different parts of the t-value formula and how they relate to your data and experiment;\n", "* you are able to calculate t-values and corresponding p-value of parameters from a GLM;\n", "\n", "**Estimated time needed to complete**: 1-3 hours
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# First some imports\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from numpy.linalg import inv\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "From your statistics classes, you might remember that many software packages (e.g. SPSS, R, SAS) do not only return beta-parameters of linear regression models, but also *t*-values and *p*-values associated with those parameters. Like beta-parameters, these statistics evaluate whether a beta-parameter (or combination of beta-parameters) differs significantly from 0 (or in fMRI terms: whether a voxel activates/deactivates significantly in response to one or more experimental factors).\n", "\n", "In univariate (activation-based) fMRI studies, we need statistics to evaluate the estimated parameters in context of the *uncertainty* of their estimation. As we'll discuss later in more detail, interpreting (and performing inference about) the magnitude of GLM parameters without their associated uncertainty is rarely warranted in univariate fMRI studies. To illustrate the problem with this, let's look at an example.\n", "\n", "In this example, we try to predict someone's height (in meters; $\\mathbf{y}$) using someone's weight (in kilos; $\\mathbf{X}$). (Note that the data is not necessarily representative of the true relationship between height and weight.)\n", "\n", "Anyway, let's run a linear regression using weight (in kilos) as a predictor for height (in meters)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "data = np.load('weight_height_data.npz')\n", "X, y = data['X'], data['y']\n", "\n", "plt.figure(figsize=(10, 6))\n", "plt.scatter(X, y)\n", "plt.title('Relation between weight and height (in meters)', y=1.05, fontsize=20)\n", "plt.xlabel('Weight (kg)', fontsize=20)\n", "plt.ylabel('Height (meters)', fontsize=20)\n", "\n", "Xn = np.hstack((np.ones((y.size, 1)), X))\n", "beta = inv(Xn.T @ Xn) @ Xn.T @ y\n", "y_hat = Xn @ beta\n", "mse = np.mean((y_hat - y) ** 2)\n", "plt.plot([X.min(), X.max()], [Xn.min(axis=0) @ beta, Xn.max(axis=0) @ beta], ls='-', c='r')\n", "plt.xlim((X.min(), X.max()))\n", "plt.text(70, 1.9, r'$\\hat{\\beta}_{weight} = %.5f$' % beta, fontsize=18)\n", "plt.text(70, 1.8, r'$MSE = %.5f$' % mse, fontsize=18)\n", "plt.grid()\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Well, quite a modest beta-parameter on the one hand, but on the other hand the Mean Squared Error is also quite low. \n", "Now, to illustrate the problem of interpretating 'raw' beta-weights, let's rephrase our objective of predicting height based on weight: we'll try to predict *height in centimeters* based on weight (still in kilos). So, what we'll do is just rescale the data points of $\\mathbf{y}$ (height in meters) so that they reflect height in centimeters. We can simply do this by multipling our $\\mathbf{y}$ by 100." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "y_cm = y * 100" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, you wouldn't expect our model to change, right? We only rescaled our target ... As you'll see below, this actually changes a lot!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (0 points): Run linear regression like the previous code block, but with y_cm instead of y as the target variable. You can use the same design (Xn). Calculate the beta-parameter and MSE (store them in the variables beta_cm and mse_cm).\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "7c0eae7138d8a700d8575f21ad675a3e", "grade": false, "grade_id": "cell-a67cba1915b72950", "locked": false, "schema_version": 3, "solution": true }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Implement the ToDo here. '''\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()\n", "print(beta_cm)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "1e224b06f09c0995243ea2afac3a24f9", "grade": true, "grade_id": "cell-60d8290f1add9cc2", "locked": true, "points": 0, "schema_version": 3, "solution": false, "task": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Tests the above ToDo'''\n", "np.testing.assert_almost_equal(beta_cm, beta * 100, decimal=4)\n", "np.testing.assert_almost_equal(mse_cm, mse * 10000, decimal=4)\n", "print(\"Well done!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you did it correctly, when you compare the beta-parameters between the two models (one where $y$ is in meters, and one where $y$ is in centimeters), you see a massive difference — a 100 fold difference to be exact\\*! This is a nice example where you see that the (raw) value of the beta-parameter is completely dependent on the scale of your variables. (Actually, you could either rescale $\\mathbf{X}$ or $\\mathbf{y}$); both will have a similar effect on your estimated beta-parameter.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToThink (0 points): Note that the MSE is a 10,000 times larger in the model with y_cm compared to y (in meters). From your understanding of how MSE is calculated, do you understand why?\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "ToThink (2 points): By now, you know that the scale of the data (either $X$ or $y$) influences the magnitude of the raw parameter estimates. One could argue that this is not relevant for fMRI data because all data (i.e. different voxels in the brain) all measure the same type of signal, so their scale shouldn't differ that much. This, however, is a false assumption.\n", "\n", "Think of (at least) two reasons why voxels might differ in their scale and write them down in the text cell below.\n", "
" ] }, { "cell_type": "markdown", "metadata": { "deletable": false, "nbgrader": { "cell_type": "markdown", "checksum": "d472b006c4fd2e42eba1bc902d4e02e3", "grade": true, "grade_id": "cell-2602edc5df20bc9f", "locked": false, "points": 2, "schema_version": 3, "solution": true } }, "source": [ "YOUR ANSWER HERE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## How to compute *t*-values and *p*-values\n", "So, you've seen that interpreting beta-parameters by themselves is useless because their value depends very much on the scale of your variables. But how should we, then, interpret the effects of our predictors on our target-variable? From the plots above, you probably guessed already that it has something to do with the MSE of our model (or, more generally, the model fit). That is indeed the case. As you might have noticed, not only the beta-parameters depend on the scale of your data, the errors (residuals) depend on the scale as well. In other words, not only the *effect* (beta-values) but also the *noise* (errors, MSE) depend on the scale of the variables! \n", "\n", "### *t*-values\n", "In fact, the key to getting interpretable effects of our predictors is to divide (\"normalize\") our beta-parameter(s) by some quantity that summarizes how well our model describes the data. This quantity is the **standard error of the beta-parameter**, usually denoted by $\\mathrm{SE}_{\\beta}$. The standard error of the beta-parameter can be computed by taking the square root of the **variance of the beta-parameter**. If we'd divide our beta-estimate with it's standard error, we compute a statistic you are all familiar with: the *t*-statistic! Formally:\n", "\n", "\\begin{align}\n", "t_{\\hat{\\beta}} = \\frac{\\hat{\\beta}}{\\mathrm{SE}_{\\hat{\\beta}}} = \\frac{\\hat{\\beta}}{\\sqrt{\\mathrm{var}(\\hat{\\beta})}}\n", "\\end{align}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToThink (0 points): Suppose that I know the $\\mathrm{SE}$ of a particular beta-parameter. How can I derive the variance of that parameter (i.e., how do I go from the $\\mathrm{SE}$ to the variance)? And yes, the answer is as straightforward as you'd think.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another way to think about it is that the t-value is the \"effect\" ($\\hat{\\beta}$) divided by your (un)certainty or confidence in the effect ($\\mathrm{SE}_{\\hat{\\beta}}$). In a way, you can think of t-values as \"uncertainty-normalized\" effects.\n", "\n", "So, what drives (statistical) uncertainty about \"effects\" (here: $\\hat{\\beta}$ parameters)? To find out, let's dissect the uncertainty term, $\\mathrm{SE}_{\\hat{\\beta}}$, a little more. The standard error of a parameter can interpreted conceptually as the \"unexplained variance of the model\" (or *noise*) multiplied with the \"design variance\" (or: *the variance of the parameter due to the design*). In this lab, we won't explain what *design variance* means or how to compute this, as this is the topic of the second notebook of this week (design_of_experiments).\n", "\n", "For now, we treat \"design variance\", here, as some known (constant) value given the design matrix ($\\mathbf{X}$). So, with this information, we can construct a conceptual formula for the standard error of our parameter(s):\n", "\n", "\\begin{align}\n", "\\mathbf{SE}_{\\hat{\\beta}} = \\sqrt{\\mathrm{noise} \\cdot \\mathrm{design\\ variance}}\n", "\\end{align}\n", "\n", "Now we also create a \"conceptual formula\" for the *t*-statistic:\n", "\n", "\\begin{align}\n", "t_{\\hat{\\beta}} = \\frac{\\hat{\\beta}}{\\mathrm{SE}_{\\hat{\\beta}}} = \\frac{\\mathrm{effect}}{\\sqrt{\\mathrm{noise} \\cdot \\mathrm{design\\ variance}}}\n", "\\end{align}\n", "\n", "**This (conceptual) formula involving effects, noise, and design variance is probably the most important concept of this course**. The effects (*t*-values) we measure in GLM analyses of fMRI data depend on two things: the effect measured ($\\hat{\\beta}$) and the (un)certainty of the effect ($SE_{\\hat{\\beta}}$), of which the latter term can be divided into the unexplained variance (\"noise\") and the design variance (uncertainty of the parameter due to the design).\n", "\n", "These two terms (noise and design variance) will be central to the next couple of weeks of this course. In this week's second notebook (topic: design of experiments), we'll focus on how to optimize our *t*-values by minimizing the \"design variance\" term. Next week (topic: preprocessing), we'll focus on how to (further) optimize our *t*-values by minimizing the error/noise.\n", "\n", "While we're going to ignore the design variance for now, we are, however, going to learn how to calculate the \"noise\" term.\n", "\n", "In fact, the noise term is *very* similar to the MSE, but instead of taking the *mean* of the squared residuals, we sum the squared residuals (\"sums of squared erros\", SSE) and divide it by the model's degrees of freedom (DF). People usually use the $\\hat{\\sigma}^{2}$ symbol for this noise term:\n", "\n", "\\begin{align}\n", "\\mathrm{noise} = \\hat{\\sigma}^{2} = \\frac{\\sum_{i=1}^{N}(\\hat{y_{i}} - y_{i})^2}{\\mathrm{df}} \n", "\\end{align}\n", "\n", "where the degrees of freedom (df) are defined as the number of samples ($N$) minus the number of predictors *including the intercept* ($P$):\n", "\n", "\\begin{align}\n", "\\mathrm{df} = N - P\n", "\\end{align}\n", "\n", "So, the formula of the *t*-statistic becomes:\n", "\n", "\\begin{align}\n", "t_{\\hat{\\beta}} = \\frac{\\hat{\\beta}}{\\sqrt{\\frac{\\sum_{i=1}^{N}(\\hat{y_{i}} - y_{i})^2}{N - P} \\cdot \\mathrm{design\\ variance}}}\n", "\\end{align}\n", "\n", "Alright, enough formulas. Let's see how we can compute these terms in Python. We're going to calculate the *t*-statistic of the weight-predictor for both models (the meter and the centimeter model) to see whether we can show that essentially the (normalized) effect of weight on height in meters is the same as the effect on height in centimeters; in other words, we are going to investigate whether the conversion to *t*-values \"normalizes\" the beta-parameters.\n", "\n", "First, we'll create a function for you to calculate the design-variance. You *don't* have to understand how this works; we're going to explain this to you in detail next week." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def design_variance(X, which_predictor=1):\n", " ''' Returns the design variance of a predictor (or contrast) in X.\n", " \n", " Parameters\n", " ----------\n", " X : numpy array\n", " Array of shape (N, P)\n", " which_predictor : int or list/array\n", " The index of the predictor you want the design var from.\n", " Note that 0 refers to the intercept!\n", " Alternatively, \"which_predictor\" can be a contrast-vector\n", " (which will be discussed later this lab).\n", " \n", " Returns\n", " -------\n", " des_var : float\n", " Design variance of the specified predictor/contrast from X.\n", " '''\n", " \n", " is_single = isinstance(which_predictor, int)\n", " if is_single:\n", " idx = which_predictor\n", " else:\n", " idx = np.array(which_predictor) != 0\n", " \n", " c = np.zeros(X.shape)\n", " c[idx] = 1 if is_single == 1 else which_predictor[idx]\n", " des_var = c.dot(np.linalg.inv(X.T.dot(X))).dot(c.T)\n", " return des_var" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, if you want the design variance of the 'weight' parameter in the varianble Xn from before, you do:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# use which_predictor=1, because the weight-column in Xn is at index 1 (index 0 = intercept)\n", "design_variance_weight_predictor = design_variance(Xn, which_predictor=1)\n", "print(\"Design variance of weight predictor is: %.6f \" % design_variance_weight_predictor)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alright, now we only need to calculate our noise-term ($\\hat{\\sigma}^2$):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Let's just redo the linear regression (for clarity)\n", "beta_meter = inv(Xn.T @ Xn) @ Xn.T @ y\n", "y_hat_meter = Xn @ beta_meter\n", "\n", "N = y.size\n", "P = Xn.shape\n", "df = (N - P)\n", "print(\"Degrees of freedom: %i\" % df)\n", "sigma_hat = np.sum((y - y_hat_meter) ** 2) / df\n", "print(\"Sigma-hat (noise) is: %.3f\" % sigma_hat)\n", "design_variance_weight = design_variance(Xn, 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can calculate the *t*-value:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t_meter = beta_meter / np.sqrt(sigma_hat * design_variance_weight)\n", "print(\"The t-value for the weight-parameter (beta = %.3f) is: %.3f\" % (beta_meter, t_meter))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's it! There's not much more to calculating *t*-values in linear regression. Now it's up to you to do the same thing and calculate the *t*-value for the model of height in centimeters, and check if it is the same as the *t*-value for the weight parameter in the model with height in meters." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (1 point): Calculate the t-statistic for the beta from the centimeter-model you calculated earlier. Store the value in a new variable named t_centimeter. Note: you don't have to calculate the design variance again (because X hasn't changed!) — you can reuse the variable design_variance_weight.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "8506e6c8ddee3462108fa37fd270437b", "grade": false, "grade_id": "cell-1b502342df415d39", "locked": false, "schema_version": 3, "solution": true }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Implement your ToDo here. '''\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "55245eef9d935d1159eeb27fbf104307", "grade": true, "grade_id": "cell-722437956591ffd0", "locked": true, "points": 1, "schema_version": 3, "solution": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Tests the above ToDo. '''\n", "try:\n", " np.testing.assert_almost_equal(t_centimeter, t_meter)\n", "except AssertionError as e:\n", " print(\"The t-value using height in centimeters is not the same as when using height in meters!\")\n", " raise(e)\n", " \n", "print(\"Well done!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### P-values\n", "As you can see, calculating *t*-values solves the \"problem\" of uninterpretable beta-parameters!\n", "\n", "Now, the last thing you need to know is how to calculate the statistical significance of your *t*-value, or in other words, how you calculate the corresponding *p*-value. You probably remember that the *p*-value corresponds to the area under the curve of a *t*-distribution associated with your observed *t*-value *and more extreme values*: \n", "![test](http://www.nku.edu/~statistics/Test_o12.gif)\n", "*Image credits: Frank Dietrich and Mike Collins, Northern Kentucky University*\n", "\n", "The function stats.t.sf(t_value, df) from the scipy package does exactly this. Importantly, this function *always* returns the right-tailed p-value. For negative t-values, however, you'd want the left-tailed *p*-value. One way to remedy this, is to always pass the absolute value of your *t*-value - np.abs(t_value) to the stats.t.sf() function. Also, the stats.t.sf() function by default returns the one-sided *p*-value. If you'd want the two-sided *p*-value, you can simply multiply the returned *p*-value by two to get the corresponding two-sided *p*-value. \n", "\n", "Let's see how we'd do that in practice:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy import stats\n", "\n", "# take the absolute by np.abs(t)\n", "p_value = stats.t.sf(np.abs(t_meter), df) * 2 # multiply by two to create a two-tailed p-value\n", "print('The p-value corresponding to t(%i) = %.3f is: %.8f' % (df, t_meter, p_value))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Contrasts\n", "We're almost done! We're really at 99% of what you should know about the GLM and fMRI analysis (except for some important caveats that have to do with GLM assumptions, that we'll discuss next week). The only major concept that we need to discuss is **contrasts**. Contrasts are basically follow-up statistical tests of GLM parameters, with which you can implement any (linear) statistical test that you are familiar with. *t*-tests, *F*-tests, ANCOVAs — they can all be realized with the GLM and the right contrast(s). (Again, if you want to know more about this equivalence between the GLM and common statistical tests, check out this [blog post](https://lindeloev.github.io/tests-as-linear/).) Importantly, the choice of contrast should reflect the hypothesis that you want to test.\n", "\n", "### *t*-tests\n", "T-tests in the GLM can be implemented in two general ways:\n", "\n", "**1. Using a contrast of a parameters \"against baseline\"**\n", "\n", "This type of contrast basically tests the hypothesis: \"Does my predictor(s) have *any* effect on my dependent variable?\" In other words, it tests the following hypothesis:\n", "* $H_{0}: \\beta = 0$ (our null-hypothesis, i.e. no effect)\n", "* $H_{a}: \\beta \\neq 0$ (our two-sided alternative hypothesis, i.e. *some* effect)\n", "\n", "Note that a directional alternative hypothesis is also possible, i.e., $H_{a}: \\beta > 0$ or $H_{a}: \\beta < 0$.\n", "\n", "**2. Using a contrast between parameters**\n", "\n", "This type of contrast basically tests hypotheses such as \"Does predictor 1 have a larger effect on my dependent variable than predictor 2?\". In other words, it tests the following hypothesis:\n", "* $H_{0}: \\beta_{1} - \\beta_{2} = 0$ (our null-hypothesis, i.e. there is no difference)\n", "* $H_{a}: \\beta_{1} - \\beta_{2} \\neq 0$ (our alternative hypotehsis, i.e. there is some difference)\n", "\n", "Let's look at an example of how we would evaluate a simple hypothesis that a beta has an *some* effect on the dependent variable. Say we'd have an experimental design with 6 conditions:\n", "\n", "* condition 1: images of **male** faces with a **happy** expression\n", "* condition 2: images of **male** faces with a **sad** expression\n", "* condition 3: images of **male** faces with a **neutral** expression\n", "* condition 4: images of **female** faces with a **happy** expression\n", "* condition 5: images of **female** faces with a **sad** expression\n", "* condition 6: images of **female** faces with a **neutral** expression\n", "\n", "Let's assume we have fMRI data from a run with 100 volumes. We then have a target-signal of shape ($100,$) and a design-matrix (after convolution with a canonical HRF) of shape ($100 \\times 7$) (the first predictor is the intercept!). We load in this data below:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "data = np.load('data_contrast_example.npz')\n", "X, y = data['X'], data['y']\n", "\n", "print(\"Shape of X: %s\" % (X.shape,))\n", "print(\"Shape of y: %s\" % (y.shape,))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After performing linear regression with these 6 predictors (after convolving the stimulus-onset times with an HRF, etc. etc.), you end up with 7 beta values:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "betas = inv(X.T @ X) @ X.T @ y\n", "betas = betas.squeeze() # remove singleton dimension; this is important for later\n", "print(\"Betas corresonding to our 6 conditions (and intercept):\\n%r\" % betas.T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first beta corresponds to the intercept, the second beta to the male/happy predictor, the third beta to the male/sad predictor, etc. etc. Now, suppose that we'd like to test whether images of male faces with a sad expression have an influence on voxel activity (our dependent variable). \n", "\n", "The first thing you need to do is extract this particular beta value from the array with beta values (I know this sounds really trivial, but bear with me):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "beta_male_sad = betas\n", "print(\"The 'extracted' beta is %.3f\" % beta_male_sad)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In neuroimaging analyses, however, this is usually done slightly differently: using **contrast-vectors**. Basically, it specifies your specific hypothesis about your beta(s) of interest in a vector. Before explaining it in more detail, let's look at it in a code example:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Again, we'd want to test whether the beta of \"male_sad\" is different from 0\n", "contrast_vector = np.array([0, 0, 1, 0, 0, 0, 0])\n", "contrast = (betas * contrast_vector).sum() # we simply elementwise multiply the contrast-vector with the betas and sum it!\n", "print('The beta-contrast is: %.3f' % contrast)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\"Wow, what a tedious way to just select the third value of the beta-array\", you might think. And, in a way, this is indeed somewhat tedious for a contrast against baseline. But let's look at a case where you would want to investigate whether two betas are different - let's say whether male sad faces have a larger effect on our voxel than male happy faces. Again, you *could* do this:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "beta_difference = betas - betas\n", "print(\"Difference between betas: %.3f\" % beta_difference)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... but you could also use a contrast-vector:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "contrast_vector = np.array([0, -1, 1, 0, 0, 0, 0])\n", "contrast = (betas * contrast_vector).sum()\n", "print('The contrast between beta 2 and beta 1 is: %.3f' % contrast)\n", "print('This is exactly the same as beta - beta: %.3f' % (betas-betas))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\"Alright, so using contrast-vectors is just a fancy way of extracting and subtracting betas from each other ...\", you might think. In a way, that's true. But you have to realize that once the hypotheses you want to test become more complicated, using contrast-vectors actually starts to make sense.\n", "\n", "Let's look at some more elaborate hypotheses. First, let's test whether male faces lead to higher voxel activity than female faces, *regardless of emotion*:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# male faces > female faces\n", "contrast_vector = [0, 1, 1, 1, -1, -1, -1]\n", "male_female_contrast = (contrast_vector * betas).sum()\n", "print(\"Male - female contrast (regardless of expression): %.2f\" % male_female_contrast)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... or whether emotional faces (regardless of *which* exact emotion) lead to higher activity than neutral faces:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Emotion (regardless of which emotion, i.e., regardless of sad/happy) - neutral\n", "contrast_vector = np.array([0, 1, 1, -2, 1, 1, -2])\n", "emo_neutral_contrast = (contrast_vector * betas).sum()\n", "print(\"Emotion - neutral contrast (regardless of which emotion): %.2f\" % emo_neutral_contrast)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "See how contrast-vectors come in handy when calculating (more intricate) comparisons? In the male-female contrast, for example, instead 'manually' picking out the betas of 'sad_male' and 'happy_male', averaging them, and subtracting their average beta from the average 'female' betas ('happy_female', 'sad_female'), you can \n", "specify a contrast-vector, multiply it with your betas, and sum them. That's it." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "ToThink (1 point): In the last contrast (emo_neural_contrast), we set all the \"emotional\" predictors (sad/happy) to 1, but the neutral predictors to minus 2 ... Why are these set to -2 and not -1? Write your answer below.\n", "
" ] }, { "cell_type": "markdown", "metadata": { "deletable": false, "nbgrader": { "cell_type": "markdown", "checksum": "5e8a98d537439a52de1fc1912f951a1b", "grade": true, "grade_id": "cell-c9b2ee94e3e03078", "locked": false, "points": 1, "schema_version": 3, "solution": true } }, "source": [ "YOUR ANSWER HERE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "ToDo (1 point): Create a contrast vector for the hypothesis: sad faces (regardless whether it's male or female) activate this voxel more than neutral faces (regardless of whether it's male/female). Multiply this contrast vector with the betas and store the result in a variable named contrast_todo.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "a89f18ad5a9781174ae700515e70ecf6", "grade": false, "grade_id": "cell-49f8094366dfb9fa", "locked": false, "schema_version": 3, "solution": true }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "# Implement the sad - neutral contrast here:\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "c6a4c18673134bd033625428dc653115", "grade": true, "grade_id": "cell-8a31e9963406e314", "locked": true, "points": 1, "schema_version": 3, "solution": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Tests the above ToDo. '''\n", "from niedu.tests.nii.week_3 import test_contrast_todo_1\n", "test_contrast_todo_1(betas, contrast_todo)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We're not only telling you about contrasts because we think it's an elegant way of computing beta-comparisons, but also because virtually every major neuroimaging software package uses them, so that you can specify what hypotheses you exactly want to test! You'll also see this when we're going to work with FSL (in week 5) to perform automated whole-brain linear regression analyses.\n", "\n", "Knowing how contrast-vectors work, we now can extend our formula for *t*-tests of beta-parameters such that they can describe **every possible test** (not only *t*-tests, but also ANOVAs, *F*-tests, etc.) of betas \"against baseline\" or between betas that you can think of: \n", "\n", "Our \"old\" formula of the *t*-test of a beta-parameter:\n", "\\begin{align}\n", "t_{\\hat{\\beta}} = \\frac{\\hat{\\beta}_{j}}{\\mathrm{SE}_{\\hat{\\beta}}}\n", "\\end{align}\n", "\n", "And now our \"generalized\" version of the *t*-test of *any* contrast/hypothesis:\n", "\n", "\\begin{align}\n", "t_{\\mathbf{c}\\hat{\\beta}} = \\frac{\\sum_{j=1}^{P}{c_{j}\\hat{\\beta}_{j}}}{\\mathrm{SE}_{\\mathbf{c}\\hat{\\beta}}} \n", "\\end{align}\n", "\n", "in which $\\mathbf{c}$ represents the entire contrast-vector, and $c_{j}$ represents the $j^{\\mathrm{th}}$ value in our contrast vector. By the way, we can simplify the (notation of the) numerator a little bit using some matrix algebra. Remember that multiplying two (equal length) vectors with each other and then summing the values together is the same thing as the (inner) \"dot product\" between the two vectors? \n", "\n", "This means that you can also evaluate this elementwise multiplication and sum of the contrast-vector and the betas using the dot-product:\n", "\n", "\\begin{align}\n", "t_{\\mathbf{c}\\hat{\\beta}} = \\frac{\\mathbf{c}\\hat{\\beta}}{\\mathrm{SE}_{\\mathbf{c}\\hat{\\beta}}} \n", "\\end{align}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (0 points): Convince yourself that the elementwise multiplication and sum is mathematically exactly the same as the dot product! Below, we initialized a hypothetical vector with beta-values (some_betas) and a hypothetical contrast-vector (some_cvec). First, implement the \"multiply and sum\" approach and then implement the \"dot product\" approach. You should find that it gives you exactly the same value: -3.34\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "some_betas = np.array([1.23, 2.95, 3.33, 4.19])\n", "some_cvec = np.array([1, 1, -1, -1])\n", "\n", "# Try to implement both approaches and convince yourself that it's\n", "# mathematically the same!\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, you need the contrast vector in the *numerator* of the *t*-value formula (i.e., $\\mathbf{c}\\hat{\\beta}$), but it turns out that you actually also need the contrast-vector in the denominator, because it's part of the calculation of design variance. Again, we will discuss how this works exactly in the next notebook. In the function design_variance, it is also possible to calculate design variance for a particular contrast (not just a single predictor) by passing a contrast vector to the which_predictor argument.\n", "\n", "We'll show this below:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# E.g., get design-variance of happy/male - sad/male\n", "c_vec = np.array([0, 1, -1, 0, 0, 0, 0]) # our contrast vector!\n", "dvar = design_variance(X, which_predictor=c_vec) # pass c_vec to which_predictor\n", "print(\"Design variance of happy/male - sad/male: %.3f\" % dvar)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the rest of ToDos this lab, make sure to pass your contrast-vector to the design_variance function in order to calculate it correctly.\n", "\n", "Now you know enough to do it yourself!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "ToDo (2 points):\n", " \n", "Calculate the *t*-value and *p*-value for the hypothesis \"sad faces have a larger effect than happy faces (regardless of gender) on our dependent variable\" (i.e. voxel activity). In other words, test the hypothesis: $\\beta_{sad} - \\beta_{happy} \\neq 0$ (note that this is a two-sided test!).\n", "\n", "Store the *t*-value and *p*-value in the variables tval_todo and pval_todo respectively. We reload the variables below (we'll call them X_new and y_new) to make sure you're working with the correct data. Note that the X_new variable already contains an intercept; the other six columns correspond to the different predictors (male/hapy, male/sad, etc.). In summary, you have to do the following:\n", "\n", "- (you don't have to calculate the betas; this has already been done (stored in the variable betas)\n", "- calculate \"sigma-hat\" ($\\mathrm{SSE} / \\mathrm{df}$)\n", "- calculate design-variance (use the design_variance function with a proper contrast-vector)\n", "- calculate the contrast ($\\mathbf{c}\\hat{\\beta}$)\n", "- calculate the t-value and p-value\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "b1c00f8bf7e944964c52b5ff63ef73e3", "grade": false, "grade_id": "cell-55833a9a2174215c", "locked": false, "schema_version": 3, "solution": true }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "data = np.load('data_contrast_example.npz')\n", "X_new, y_new = data['X'], data['y']\n", "\n", "print(\"Shape of X: %s\" % (X_new.shape,))\n", "print(\"Shape of y: %s\" % (y_new.shape,))\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "d1335cf1eee1eac6e83cfb94d6a65872", "grade": true, "grade_id": "cell-411fa3ab43d50400", "locked": true, "points": 2, "schema_version": 3, "solution": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Part 1 of testing the above ToDo. ''' \n", "from niedu.tests.nii.week_3 import test_contrast_todo_2\n", "test_contrast_todo_2(X_new, y_new, betas, tval_todo, pval_todo)\n", "print(\"Well done!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### *F*-tests on contrasts\n", "In the previous section we discussed how to calculate *t*-values for single contrasts. However, sometimes you might have a hypothesis about multiple contrasts at the same time. This may sound weird, but let's consider an experiment.\n", "\n", "Suppose you have data from an experiment in which you showed images circles which were either blue, red, or green. In that case, you have three predictors. Then, you could have very specific question, like \"Do blue circles activate a voxel significantly compared to baseline\", which corresponds to the following null and alternative hypothesis:\n", "\n", "* $H_{0}: \\beta_{blue} = 0$ (our null-hypothesis, i.e. there is no activation compared to baseline)\n", "* $H_{a}: \\beta_{blue} > 0$ (our alternative hypothesis, i.e. blue activates relative to baseline)\n", "\n", "However, you can also have a more general question, like \"Does the presentation of *any* circle (regardless of color) activate a voxel compared to baseline?\". This question represents the following null and alternative hypothesis:\n", "\n", "* $H_{0}: \\beta_{blue} = \\beta_{red} = \\beta_{green} = 0$\n", "* $H_{a}: (\\beta_{blue} > 0) \\vee (\\beta_{red} > 0) \\vee (\\beta_{green} > 0)$\n", "\n", "The $\\vee$ symbol in the alternative hypothesis means \"or\". So the alternative hypothesis nicely illustrates our question: is there *any* condition (circle) that activates a voxel more than baseline? This hypothesis-test might sound familiar, because it encompasses the *F*-test! In other words, an *F*-test tests *a collection of contrasts* together. In the example here, the *F*-test tests the following contrasts together (ignoring the intercept) of our beta-parameters:\n", "\n", "* [1, 0, 0] ($\\mathrm{red} > 0$)\n", "* [0, 1, 0] ($\\mathrm{blue} > 0$)\n", "* [0, 0, 1] ($\\mathrm{green} > 0$)\n", "\n", "Thus, a *F*-test basically tests this contrast-*matrix* all at once! Therefore, the *F*-tests is a type of \"omnibus test\"! \n", "\n", "Now, let's look at the math behind the *F*-statistic. The *F*-statistic for set of $K$ contrasts (i.e., the number of rows in the contrast-matrix) is defined as follows:\n", "\n", "\\begin{align}\n", "F = (\\mathbf{c}\\hat{\\beta})^{T}[K\\mathbf{c}((X^{T}X)^{-1}\\hat{\\sigma}^{2})\\mathbf{c}^{T}]^{-1}(\\mathbf{c}\\hat{\\beta})\n", "\\end{align}\n", "\n", "With a little imagination, you can see how the *F*-test is an extension of the *t*-test of a single contrast to accomodate testing a set of contrasts together. Don't worry, you don't have to understand how the formula for the *F*-statistic works mathematically and you don't have to implement this in Python. But you *do* need to understand what type of hypothesis an *F*-test tests! \n", "\n", "Let's practice this in a ToDo!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "ToDo (1 point)\n", " \n", "Remember the temporal basis sets from before? Suppose we have an experiment with two conditions (\"A\" and \"B\") and suppose we've created a design matrix based on convolution with a single-gamma basis set (with a canonical HRF, its temporal derivative, and its dispersion derivative). Together with the intercept, the design matrix thus has 7 columns (2 conditions * 3 HRF + intercept).\n", "\n", "The order of the columns is as follows:\n", "* column 1: intercept\n", "* column 2: canonical HRF \"A\"\n", "* column 3: temporal deriv \"A\"\n", "* column 4: dispersion deriv \"A\"\n", "* column 5: canonical HRF \"B\"\n", "* column 6: temporal deriv \"B\"\n", "* column 7: dispersion deriv \"B\"\n", "\n", "Suppose I want to test whether there is *any* difference in response to condition \"A\" ($A > 0$) compared to baseline, and *I don't care what element of the HRF caused it*. I can use an F-test for this. What would the corresponding contrast-*matrix* (in which each row represents a different contrast) look like? \n", "\n", "We've created an 'empty' (all-zeros) 2D matrix below with three rows. It's up to you to fill in the matrix such that it can be used to test the above question/hypothesis.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "b4ec5445e719a70ea9483902eb83327a", "grade": false, "grade_id": "cell-82c295ab029883fe", "locked": false, "schema_version": 3, "solution": true }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "# Fill in the correct values!\n", "contrast_matrix = np.array([\n", "\n", " [0, 0, 0, 0, 0, 0, 0],\n", " [0, 0, 0, 0, 0, 0, 0],\n", " [0, 0, 0, 0, 0, 0, 0]\n", "\n", "])\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "8dff3dd4e20a31d27684ad772be6e827", "grade": true, "grade_id": "cell-7f33290dbbfa631d", "locked": true, "points": 1, "schema_version": 3, "solution": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "''' Tests the above ToDo. '''\n", "from niedu.tests.nii.week_3 import test_definition_ftest\n", "\n", "test_definition_ftest(contrast_matrix)\n", "print(\"Well done!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Summary\n", "Alright, now you know basically everything about how to perform a univariate fMRI analysis! \n", "\n", "\"Wait, that's it?\", you might ask (or not). Well, yeah, regular univariate analyses as you might read about in scientific journals do basically what you've just learned, but then not on a single voxel, but on each voxel in the brain separately. Basically just a gigantic for-loop across voxels in which everytime the same design ($\\mathbf{X}$) is used to predict a new voxel-signal ($\\mathbf{y}$). Afterwards, the *t*-values of the contrast (hypothesis) you're interested in are plotted back onto a brain, color-code it (high t-values yellow, low t-values red), and voilà, you have your pretty brain plot." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToThink (1 point): More explained variance (i.e., a smaller \"sums of squared error\" term) does not always mean that your t-value is higher. Explain how this might happen.\n", "
" ] }, { "cell_type": "markdown", "metadata": { "deletable": false, "nbgrader": { "cell_type": "markdown", "checksum": "f10eabffe20ce417ef8507044d2155af", "grade": true, "grade_id": "cell-50d9ec0c4060b7ff", "locked": false, "points": 1, "schema_version": 3, "solution": true, "task": false } }, "source": [ "YOUR ANSWER HERE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " ToDo (2 points): Suppose that, within the hypothesized face-experiment explained earlier, you want to know which parts of the brain show (significantly) more activity during periods without stimuli (i.e., no faces were shown, i.e., \"rest\") than during periods with stimuli. Define a contrast vector which would test this hypothesis and store it in a variable cvec_rest. Remember: the original face experiment had 7 predictors (the first one being the intercept, followed by 6 face predictors).\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "nbgrader": { "cell_type": "code", "checksum": "642023e320be2c5fc4473b0581bf3f12", "grade": false, "grade_id": "cell-fc092da92c99eae3", "locked": false, "schema_version": 3, "solution": true, "task": false }, "tags": [ "raises-exception", "remove-output" ] }, "outputs": [], "source": [ "# Implement the assignment here\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "cell_type": "code", "checksum": "78529115950f515d3fd84b33cde34190", "grade": true, "grade_id": "cell-7ae9b02813f2102f", "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_3 import test_rest_vs_stim_contrast\n", "test_rest_vs_stim_contrast(cvec_rest)" ] }, { "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": {}, "number_sections": true, "sideBar": true, "skip_h1_title": true, "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 }