{ "metadata": { "name": "", "signature": "sha256:24163e36ecf41cc409dc9bbc13ebfd7d2d0ec7a8ca6a033408219d57475fbe0c" }, "name": "regression", "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "code", "collapsed": false, "input": [ "# - compatibility with Python 3\n", "from __future__ import print_function # print('me') instead of print 'me'\n", "from __future__ import division # 1/2 == 0.5, not 0" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - show figures inside the notebook\n", "%matplotlib inline" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - import common modules\n", "import numpy as np # the Python array package\n", "import matplotlib.pyplot as plt # the Python plotting package" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "code", "collapsed": false, "input": [ "# - a new import for us - import numpy.linalg with a shorter name\n", "import numpy.linalg as npl" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Using the utility functions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need to put the utility functions on the Python search path.\n", "\n", "Import the `sys` module, and append the full path to the `pna_code` directory to `sys.path`. If you need to jog your memory on how to do this, have a look at the top of `day4/first_statistics_solutions.ipynb` notebook." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# import sys module\n", "# append full path of \"pna_code\" directory to sys.path" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check you can import `spm_funcs` and `stimuli`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Try importing \"spm_funcs\" and \"stimuli\" modules" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a quick check, make a time vector from 0 up to (not including) 25, in steps of 0.1.\n", "\n", "Use the `spm_hrf` function in `spm_funcs` to sample the SPM HRF at these times.\n", "\n", "Plot the HRF samples against time to make sure it gives a reasonable shape." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Make times vector from 0 to 25 in steps of 0.1\n", "# Sample spm_hrf function at these times\n", "# Plot samples against time" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are going to use the same FMRI image and stimulus file as last week." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Import nibabel and load the image `ds114_sub009_t2r1.nii`. Get the number of TRs from the image shape." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Load 'ds114_sub009_t2r1.nii' with nibabel\n", "# Get number of TRs from the image shape" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use the IPython `?` at the end of the function name to check the arguments needed for `stimuli.events2neural`. (Where does the help text come from?)" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Use question mark to get help for stimuli.events2neural" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The stimulus file is `ds114_sub009_t2r1_cond.txt`.\n", "\n", "The TR for this run was 2.5 seconds." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# - TR for this run\n", "tr = 2.5" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use the `stimuli.events2neural` function to load this stimulus file and create a neural time-course prediction for this run. Plot the prediction.\n", "\n", "If you get an error here - call your friendly instructor. Maybe you didn't fix the bug from the first exercise?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Create neural time course prediction with \"events2neural\" function\n", "# Plot the prediction" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Drop the first volume" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you remember, this run has a very strange first volume. Load the image data into an array. Drop the first volume from the image data. Drop the first element from the neural prediction." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Load the image data into an array\n", "# Drop the first volume from the data array\n", "# Drop the first element from the neural prediction vector" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Making the hemodynamic prediction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need to use our new convolution skills to make a hemodynamic prediction from the neural prediction.\n", "\n", "To convolve, we need to sample the HRF at times matching the TRs.\n", "\n", "Make a time vector from 0 to 25 in steps of the TR (2.5). Sample the SPM HRF at these times. Plot to check it looks OK." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Make time vector to take HRF samples (0 to 25 in steps of TR)\n", "# Sample SPM HRF at these times\n", "# Plot the samples against time" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Convolve the fixed (dropped first value) neural predictor with the HRF samples. \n", "\n", "Drop the last (M-1) set of values which are the tail of the convolution (where `M = len(hrf_samples)`). Plot." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Convolve the fixed neural predictor with the HRF samples\n", "# Drop the tail of the convolution result so your hemodynamic\n", "# prediction is the same length as the neural prediction.\n", "# Plot" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Regressing the voxel time-course against the prediction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are going to use the same voxel time course as last week.\n", "\n", "That was the voxel with coordinates (45, 43, 14).\n", "\n", "Use slicing into the (fixed, first volume dropped) image data to get the voxel value for this voxel, for every time point (volume).\n", "\n", "Call this vector of voxel values `y` (our data vector)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Get the values for voxel coordinate (45, 43, 14) and every volume\n", "# Plot the values (voxel time course)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are going to do simple linear regression, regressing the data `y` against the hemodynamic prediction." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As experienced regressers, we first do a scatterplot of the data `y` on the y-axis, against the hemodynamic prediction on the x-axis. Remember you will need to use not-default marker to get a scatter plot, instead of points joined by lines." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Do a scatter plot of the data against the hemodynamic predictor" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make `X` be a 2D array of 1s, with first dimension length that is the length of `y`, and second dimension length of 2." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# 2D array 'X' that will be our design matrix\n", "# An array of 1s shape (N, 2) where N = the length of y" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fill the first column in `X` with the hemodynamic prediction values.\n", "\n", "Use `plt.imshow` to show the first 40 rows of the design as an image (if you show the whole thing, it will be too long and narrow to see what's going on)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Put hemodynamic prediction into first column of X\n", "# Show the first 40 rows of the design matrix X as an image" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use the `npl.pinv` function to get the pseudo-inverse of the design `X`. Take the matrix multiplication (dot product) of this pseudo-inverse with the data `y` to get the estimated parameters `B`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Using numpy.linalg 'pinv' operator to get pseudo-inverse of the design\n", "# Get B by taking dot product of pseudo-inverse with the data y" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first value of `B` is the slope (parameter for the hemodynamic regressor), and the second is the intercept (parameter for the column of ones).\n", "\n", "The *predicted data* $\\hat{y}$ is given by $X B$ - the matrix multiplication of the design `X` on the beta estimates `B`. Calculate the predicted data `y_hat`." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Use X and B to make predicted data y_hat" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Repeat the scatterplot of the data `y` against the hemodynamic prediction values.\n", "\n", "On the same plot, plot the `y_hat` values against the hemodynamic prediction values. " ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Repeat scatterplot of the data y against the hemodynamic predictor\n", "# On the same plot, plot the values for my_line(x) against the hemodynamic predictor" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The residuals are just the data `y` minus the predicted data `y_hat`. Calculate the residuals, and show the standard deviation of the residuals:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Calculate the residuals E by subtracing the prediction from the data\n", "# Calculate the standard deviation of the residuals" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Improving the model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Looking again at the data scatterplot, it seems that there is a linear drift of the signal over time.\n", "\n", "Maybe we can put in an extra regressor to remove that effect." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make a new regressor, the same length as `y`, with values increasing linearly from the first to the last value. For example, you could do this by using `np.linspace` to get N = (length of y) samples between -1 and 1." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Make a regressor that increases linearly from the first to the last value" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Make a new design matrix `X_better` of 1s, with shape N, 3, where N is the length of `y`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Make 2D array X_better with shape (length of y, 3)" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fill the first column of `X_better` with the hemodynamic prediction as before.\n", "\n", "Fill the second column with the new linear regressor.\n", "\n", "Show the first 40 rows of the design as an image." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Set first column of X_better as hemodynamic prediction.\n", "# Set second column of X_better as linear regressor.\n", "# Show the first 40 rows of the design as an image." ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use `pinv` to fit this new design `X_better` to the data `y`, giving parameter vector `B_better`. Print the new `B_better` parameters:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Use pinv to estimate best fit parameter vector B for this design\n", "# Print the parameters" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The parameters are (in order):\n", "\n", "* slope for hemodynamic regressor;\n", "* slope for linear trend regressor;\n", "* intercept.\n", "\n", "Is the slope for the linear trend (second parameter) negative or positive? Why?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate the new predicted values `y_hat`. Calculate the new residuals. Has the standard deviation of the residuals gone up or down?" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Calculate predicted data y_hat.\n", "# Calculate residuals\n", "# Show standard deviation of residuals" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Contrasts and the t-statistic" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you remember, contrasts are vectors that combine parameters. We use them to test hypotheses about the parameters.\n", "\n", "For example, maybe we want to test if the slope for the hemodynamic predictor is significantly different from zero.\n", "\n", "The contrast $[1, 0, 0]$ will select the first parameter, when you do the dot product of the contrast with the parameter vector `B_better`." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# - contrast to select first parameter (of 3)\n", "c = np.array([1, 0, 0])\n", "c" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Take the dot product of `c` with the parameter vector `B_better`. This is `1 * B_better[0] + 0 * B_better[1] + 0 * B_better[2]` == `B_better[0]`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# t_top = dot product of c with B_better" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now have the top half of the t-statistic. The formula for the t-statistic is:\n", "\n", "$$\n", "t = \\frac{c^T \\hat\\beta}{\\sqrt{\\mathrm{var}(c^T \\hat\\beta)}}\n", "$$\n", "\n", "where $\\mathrm{var}(c^T \\hat\\beta)$ is the expected variance for $c^T \\hat\\beta$.\n", "\n", "The expected variance - it turns out - works out to this (see the denominator):\n", "\n", "$$\n", "t = \\frac{c^T \\hat\\beta}{\\sqrt{\\hat{\\sigma}^2 c^T (X^T X)^+ c}}\n", "$$\n", "\n", "where $\\hat{\\sigma}^2$ is our estimate of variance in the residuals, and $(X^T X)^+$ is the [pseudo-inverse](https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse) of $X^T X$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So, our next job is calculate $\\hat{\\sigma}^2$. This is our estimate of the variance in the residuals. The best estimate is given by the *residual sum of squares* ($\\sum_{i=1}^n{e_i^2}$) divided by the *degrees of freedom*. The *degrees of freedom* in our case are the number of observations (== number of rows == `len(y)`) minus the number of independent columns in the design matrix `X_better`. In fact all three columns in the design are independent from each other, so the number of independent columns is 3.\n", "\n", "Calculate the residual sum of squares, and divide by the degrees of freedom to give `sigma_2_hat`:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Calculate the residual sum of squares = rss\n", "# Calculate the degrees of freedom = df\n", "# Divide the rss by df to give the estimated sigma squared" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we need $c^T (X^T X)^+ c$.\n", "\n", "$(X^T X)^+$ is just the `pinv` of the dot product of `X_better.T` with `X_better`. This gives a measure of the covariation of the design columns. Calculate that:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Calculate pseudoinverse of (X.T X) where X is X_better\n", "# Call the result pXtX\n", "# Display pXtX" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now take the dot product of `c`, onto `pXtX`, and the dot product of the result, onto `c` again. This might be called `bcov` - referring to the covariance of the beta (parameter) estimates." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Calculate bcov = c (X.T X)+ c, and display it" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is the top-left value from the `pXtX` array. Is it obvious why, from the matrix multiplication with the `c` vector?\n", "\n", "Now we can make the bottom half of the t statistic, by taking the square root of (`sigma_2_hat * bcov`)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Square root of sigma_2_hat * c (X.T X)+ c" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we get the t statistic by dividing the top by the bottom:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Divide t_top by t_bottom" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's check that with the `t_stat` function from the [GLM introduction page](http://perrin.dynevor.org/glm_intro.html).\n", "\n", "Here is the `t_stat` function:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Get t distribution code from scipy library" ], "language": "python", "metadata": {}, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Call the `t_stat` function, passing the data, design and contrast, and see if you get back the same t-statistic that you calculated before:" ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Try getting the B vector, t-statistic, df, and p value from this function\n", "# Check you get the same t value as you calculated above" ], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }