1. pna2015
  2. day6
In [1]:
# - 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
In [2]:
# - show figures inside the notebook
%matplotlib inline
In [3]:
# - import common modules
import numpy as np  # the Python array package
import matplotlib.pyplot as plt  # the Python plotting package
In [4]:
# - a new import for us - import numpy.linalg with a shorter name
import numpy.linalg as npl

Using the utility functions

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.

In [5]:
# import sys module
# append full path of "pna_code" directory to sys.path
import sys
sys.path.append('/Users/mb312/pna_code')

Check you can import spm_funcs and stimuli:

In [6]:
# Try importing "spm_funcs" and "stimuli" modules
import spm_funcs
import stimuli

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.

In [7]:
# Make times vector from 0 to 25 in steps of 0.1
# Sample spm_hrf function at these times
# Plot samples against time
times = np.arange(0, 25, 0.1)
hrf_samples = spm_funcs.spm_hrf(times)
plt.plot(times, hrf_samples)
Out[7]:
[<matplotlib.lines.Line2D at 0x108016290>]

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.

In [8]:
# Load 'ds114_sub009_t2r1.nii' with nibabel
# Get number of TRs from the image shape
import nibabel as nib
img = nib.load('ds114_sub009_t2r1.nii')
img.shape
Out[8]:
(64, 64, 30, 173)

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?)

In [9]:
# Use question mark to get help for stimuli.events2neural
#: stimuli.events2neural?

The stimulus file is ds114_sub009_t2r1_cond.txt.

The TR for this run was 2.5 seconds.

In [10]:
# - 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?

In [11]:
# Create neural time course prediction with "events2neural" function
# Plot the prediction
neural_prediction = stimuli.events2neural('ds114_sub009_t2r1_cond.txt', tr, img.shape[-1])
plt.plot(neural_prediction)
plt.ylim(0, 1.2)
Out[11]:
(0, 1.2)

Drop the first volume

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.

In [12]:
# Load the image data into an array
# Drop the first volume from the data array
# Drop the first element from the neural prediction vector
data = img.get_data()
fixed_data = data[..., 1:]
fixed_neural = neural_prediction[1:]

Making the hemodynamic prediction

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.

In [13]:
# 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
sample_times = np.arange(0, 25, 2.5)
hrf_samples = spm_funcs.spm_hrf(sample_times)
plt.plot(sample_times, hrf_samples)
Out[13]:
[<matplotlib.lines.Line2D at 0x10807e710>]

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.

In [14]:
# 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
convolved = np.convolve(fixed_neural, hrf_samples)
hemodynamic = convolved[:len(fixed_neural)]
plt.plot(hemodynamic)
Out[14]:
[<matplotlib.lines.Line2D at 0x108456ed0>]

Regressing the voxel time-course against the prediction

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).

In [15]:
# Get the values for voxel coordinate (45, 43, 14) and every volume
# Plot the values (voxel time course)
y = fixed_data[45, 43, 14, :]
plt.plot(y)
Out[15]:
[<matplotlib.lines.Line2D at 0x10ae972d0>]

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.

In [16]:
# Do a scatter plot of the data against the hemodynamic predictor
plt.plot(hemodynamic, y, '+')
Out[16]:
[<matplotlib.lines.Line2D at 0x10afbb5d0>]

Make X be a 2D array of 1s, with first dimension length that is the length of y, and second dimension length of 2.

In [17]:
# 2D array 'X' that will be our design matrix
# An array of 1s shape (N, 2) where N = the length of y
X = np.ones((len(y), 2))

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).

In [18]:
# Put hemodynamic prediction into first column of X
# Show the first 40 rows of the design matrix X as an image
X[:, 0] = hemodynamic
plt.imshow(X[:40], cmap='gray', interpolation='nearest')
Out[18]:
<matplotlib.image.AxesImage at 0x10b0e1390>

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:

In [19]:
# 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
B = npl.pinv(X).dot(y)
B
Out[19]:
array([   19.72782935,  2344.60385078])

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.

In [20]:
# Use X and B to make predicted data y_hat
y_hat = X.dot(B)

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.

In [21]:
# 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
plt.plot(hemodynamic, y, '+')
plt.plot(hemodynamic, y_hat)
Out[21]:
[<matplotlib.lines.Line2D at 0x10b128290>]

The residuals are just the data y minus the predicted data y_hat. Calculate the residuals, and show the standard deviation of the residuals:

In [22]:
# Calculate the residuals E by subtracing the prediction from the data
# Calculate the standard deviation of the residuals
E = y - y_hat
np.std(E)
Out[22]:
memmap(25.72736744062518)

Improving the model

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.

In [23]:
# Make a regressor that increases linearly from the first to the last value
linear = np.linspace(-1, 1, len(y))

Make a new design matrix X_better of 1s, with shape N, 3, where N is the length of y:

In [24]:
# Make 2D array X_better with shape (length of y, 3)
X_better = np.ones((len(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.

In [25]:
# 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.
X_better[:, 0] = hemodynamic
X_better[:, 1] = linear
plt.imshow(X_better[:40], cmap='gray', interpolation='nearest')
Out[25]:
<matplotlib.image.AxesImage at 0x10b2b8e50>

Use pinv to fit this new design X_better to the data y, giving parameter vector B_better. Print the new B_better parameters:

In [26]:
# Use pinv to estimate best fit parameter vector B for this design
# Print the parameters
B_better = npl.pinv(X_better).dot(y)
B_better
Out[26]:
array([   17.20561032,   -34.16623799,  2345.83563217])

The parameters are (in order):

  • slope for hemodynamic regressor;
  • slope for linear trend regressor;
  • intercept.

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?

In [27]:
# Calculate predicted data y_hat.
# Calculate residuals
# Show standard deviation of residuals
y_hat = X_better.dot(B_better)
E = y - y_hat
np.std(E)
Out[27]:
memmap(16.42808837342474)

Contrasts and the t-statistic

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.

In [28]:
# - contrast to select first parameter (of 3)
c = np.array([1, 0, 0])
c
Out[28]:
array([1, 0, 0])

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]:

In [29]:
# t_top = dot product of c with B_better
t_top = c.dot(B_better)
t_top
Out[29]:
17.205610322002929

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:

In [30]:
# Calculate the residual sum of squares = rss
# Calculate the degrees of freedom = df
# Divide the rss by df to give the estimated sigma squared
rss = np.sum(E ** 2)
df = len(y) - 3
sigma_2_hat = rss / df
sigma_2_hat
Out[30]:
memmap(274.67289389390027)

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:

In [31]:
# Calculate pseudoinverse of (X.T X) where X is X_better
# Call the result pXtX
# Display pXtX
pXtX = npl.pinv(X_better.T.dot(X_better))
pXtX
Out[31]:
array([[ 0.02256102,  0.00127805, -0.01101817],
       [ 0.00127805,  0.01731262, -0.00062416],
       [-0.01101817, -0.00062416,  0.01119492]])

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.

In [32]:
# Calculate bcov = c (X.T X)+ c, and display it
bcov = c.dot(pXtX).dot(c)
bcov
Out[32]:
0.022561023805170197

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).

In [33]:
# Square root of sigma_2_hat * c (X.T X)+ c
t_bottom = np.sqrt(sigma_2_hat * bcov)
t_bottom
Out[33]:
memmap(2.489357687793233)

Finally we get the t statistic by dividing the top by the bottom:

In [34]:
# Divide t_top by t_bottom
t_top / t_bottom
Out[34]:
6.9116665742219494

Let's check that with the t_stat function from the GLM introduction page.

Here is the t_stat function:

In [35]:
# Get t distribution code from scipy library
from scipy.stats import t as t_dist

def t_stat(y, X, c):
    """ betas, t statistic and significance test given data, design matrix, contrast
    
    This is OLS estimation; we assume the errors to have independent
    and identical normal distributions around zero for each $i$ in 
    $\e_i$ (i.i.d).
    """
    # Make sure y, X, c are all arrays
    y = np.asarray(y)
    X = np.asarray(X)
    c = np.atleast_2d(c).T  # As column vector
    # Calculate the parameters - b hat
    beta = npl.pinv(X).dot(y)
    # The fitted values - y hat
    fitted = X.dot(beta)
    # Residual error
    errors = y - fitted
    # Residual sum of squares
    RSS = (errors**2).sum(axis=0)
    # Degrees of freedom is the number of observations n minus the number
    # of independent regressors we have used.  If all the regressor 
    # columns in X are independent then the (matrix rank of X) == p
    # (where p the number of columns in X). If there is one column that 
    # can be expressed as a linear sum of the other columns then 
    # (matrix rank of X) will be p - 1 - and so on.
    df = X.shape[0] - npl.matrix_rank(X)
    # Mean residual sum of squares
    MRSS = RSS / df
    # calculate bottom half of t statistic
    SE = np.sqrt(MRSS * c.T.dot(npl.pinv(X.T.dot(X)).dot(c)))
    t = c.T.dot(beta) / SE
    # Get p value for t value using cumulative density dunction
    # (CDF) of t distribution 
    ltp = t_dist.cdf(t, df) # lower tail p
    p = 1 - ltp # upper tail p
    return beta, t, df, p

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:

In [36]:
# 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
b, t, df, p = t_stat(y, X_better, c)
t
Out[36]:
array([[ 6.91166657]])