# - compatibility with Python 3
from __future__ import print_function # print('me') instead of print 'me'
from __future__ import division # 1/2 == 0.5, not 0
# - show figures inside the notebook
%matplotlib inline
# - import common modules
import numpy as np # the Python array package
import matplotlib.pyplot as plt # the Python plotting package
# - a new import for us - import numpy.linalg with a shorter name
import numpy.linalg as npl
We need to put the utility functions on the Python search path.
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.
# import sys module
# append full path of "pna_code" directory to sys.path
Check you can import spm_funcs and stimuli:
# Try importing "spm_funcs" and "stimuli" modules
As a quick check, make a time vector from 0 up to (not including) 25, in steps of 0.1.
Use the spm_hrf function in spm_funcs to sample the SPM HRF at these times.
Plot the HRF samples against time to make sure it gives a reasonable shape.
# Make times vector from 0 to 25 in steps of 0.1
# Sample spm_hrf function at these times
# Plot samples against time
We are going to use the same FMRI image and stimulus file as last week.
Import nibabel and load the image ds114_sub009_t2r1.nii. Get the number of TRs from the image shape.
# Load 'ds114_sub009_t2r1.nii' with nibabel
# Get number of TRs from the image shape
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?)
# Use question mark to get help for stimuli.events2neural
The stimulus file is ds114_sub009_t2r1_cond.txt.
The TR for this run was 2.5 seconds.
# - TR for this run
tr = 2.5
Use the stimuli.events2neural function to load this stimulus file and create a neural time-course prediction for this run. Plot the prediction.
If you get an error here - call your friendly instructor. Maybe you didn't fix the bug from the first exercise?
# Create neural time course prediction with "events2neural" function
# Plot the prediction
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.
# Load the image data into an array
# Drop the first volume from the data array
# Drop the first element from the neural prediction vector
We need to use our new convolution skills to make a hemodynamic prediction from the neural prediction.
To convolve, we need to sample the HRF at times matching the TRs.
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.
# Make time vector to take HRF samples (0 to 25 in steps of TR)
# Sample SPM HRF at these times
# Plot the samples against time
Convolve the fixed (dropped first value) neural predictor with the HRF samples.
Drop the last (M-1) set of values which are the tail of the convolution (where M = len(hrf_samples)). Plot.
# Convolve the fixed neural predictor with the HRF samples
# Drop the tail of the convolution result so your hemodynamic
# prediction is the same length as the neural prediction.
# Plot
We are going to use the same voxel time course as last week.
That was the voxel with coordinates (45, 43, 14).
Use slicing into the (fixed, first volume dropped) image data to get the voxel value for this voxel, for every time point (volume).
Call this vector of voxel values y (our data vector).
# Get the values for voxel coordinate (45, 43, 14) and every volume
# Plot the values (voxel time course)
We are going to do simple linear regression, regressing the data y against the hemodynamic prediction.
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.
# Do a scatter plot of the data against the hemodynamic predictor
Make X be a 2D array of 1s, with first dimension length that is the length of y, and second dimension length of 2.
# 2D array 'X' that will be our design matrix
# An array of 1s shape (N, 2) where N = the length of y
Fill the first column in X with the hemodynamic prediction values.
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).
# Put hemodynamic prediction into first column of X
# Show the first 40 rows of the design matrix X as an image
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:
# Using numpy.linalg 'pinv' operator to get pseudo-inverse of the design
# Get B by taking dot product of pseudo-inverse with the data y
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).
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.
# Use X and B to make predicted data y_hat
Repeat the scatterplot of the data y against the hemodynamic prediction values.
On the same plot, plot the y_hat values against the hemodynamic prediction values.
# Repeat scatterplot of the data y against the hemodynamic predictor
# On the same plot, plot the values for my_line(x) against the hemodynamic predictor
The residuals are just the data y minus the predicted data y_hat. Calculate the residuals, and show the standard deviation of the residuals:
# Calculate the residuals E by subtracing the prediction from the data
# Calculate the standard deviation of the residuals
Looking again at the data scatterplot, it seems that there is a linear drift of the signal over time.
Maybe we can put in an extra regressor to remove that effect.
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.
# Make a regressor that increases linearly from the first to the last value
Make a new design matrix X_better of 1s, with shape N, 3, where N is the length of y:
# Make 2D array X_better with shape (length of y, 3)
Fill the first column of X_better with the hemodynamic prediction as before.
Fill the second column with the new linear regressor.
Show the first 40 rows of the design as an image.
# Set first column of X_better as hemodynamic prediction.
# Set second column of X_better as linear regressor.
# Show the first 40 rows of the design as an image.
Use pinv to fit this new design X_better to the data y, giving parameter vector B_better. Print the new B_better parameters:
# Use pinv to estimate best fit parameter vector B for this design
# Print the parameters
The parameters are (in order):
Is the slope for the linear trend (second parameter) negative or positive? Why?
Calculate the new predicted values y_hat. Calculate the new residuals. Has the standard deviation of the residuals gone up or down?
# Calculate predicted data y_hat.
# Calculate residuals
# Show standard deviation of residuals
As you remember, contrasts are vectors that combine parameters. We use them to test hypotheses about the parameters.
For example, maybe we want to test if the slope for the hemodynamic predictor is significantly different from zero.
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.
# - contrast to select first parameter (of 3)
c = np.array([1, 0, 0])
c
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]:
# t_top = dot product of c with B_better
We now have the top half of the t-statistic. The formula for the t-statistic is:
$$ t = \frac{c^T \hat\beta}{\sqrt{\mathrm{var}(c^T \hat\beta)}} $$where $\mathrm{var}(c^T \hat\beta)$ is the expected variance for $c^T \hat\beta$.
The expected variance - it turns out - works out to this (see the denominator):
$$ t = \frac{c^T \hat\beta}{\sqrt{\hat{\sigma}^2 c^T (X^T X)^+ c}} $$where $\hat{\sigma}^2$ is our estimate of variance in the residuals, and $(X^T X)^+$ is the pseudo-inverse of $X^T X$.
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.
Calculate the residual sum of squares, and divide by the degrees of freedom to give sigma_2_hat:
# Calculate the residual sum of squares = rss
# Calculate the degrees of freedom = df
# Divide the rss by df to give the estimated sigma squared
Now we need $c^T (X^T X)^+ c$.
$(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:
# Calculate pseudoinverse of (X.T X) where X is X_better
# Call the result pXtX
# Display pXtX
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.
# Calculate bcov = c (X.T X)+ c, and display it
This is the top-left value from the pXtX array. Is it obvious why, from the matrix multiplication with the c vector?
Now we can make the bottom half of the t statistic, by taking the square root of (sigma_2_hat * bcov).
# Square root of sigma_2_hat * c (X.T X)+ c
Finally we get the t statistic by dividing the top by the bottom:
# Divide t_top by t_bottom
Let's check that with the t_stat function from the GLM introduction page.
Here is the t_stat function:
# Get t distribution code from scipy library
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:
# Try getting the B vector, t-statistic, df, and p value from this function
# Check you get the same t value as you calculated above