# - 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
# - set gray colormap and nearest neighbor interpolation by default
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'
# - import numpy.linalg with a shorter name
import numpy.linalg as npl
# - import nibabel
import nibabel as nib
# - import stimuli, spm_funcs from pna_code
# Call over your friendly instructor if you get an error importing
# You might want to make sure you have the 'for-class-day9' branch checked out in pna_code
import stimuli
import spm_funcs
Estimating for the whole brain¶
Here we load the processed dataset. The dataset has had the following transformations:
- drop first 4 volumes;
- slice time correction;
- motion correction;
- coregistration of structural to functional scans;
- registration of structural scan to template;
- resampling of functional and structural scans to match template.
The processed functional image is wafds114_sub009_t2r1.nii.
See nipype_ds114_sub009_t2r1.py for details.
# - load 'wafds114_sub009_t2r1.nii'
# Get 4D data
img = nib.load('wafds114_sub009_t2r1.nii')
data = img.get_data()
data.shape
# - get number of volumes, 3D vol shape, number of voxels in one volume
n_vols = data.shape[-1]
vol_shape = data.shape[:3]
n_voxels = np.prod(vol_shape)
During processing, we dropped the first 4 volumes:
# - number of volumes dropped from the original 4D data
n_dropped = 4
We are going to run our linear model on these data. To run the model we want a two-dimensional array with one column for each voxel. The voxel time courses are therefore in the columns. The array will be shape (n_vols, n_voxels).
Be careful - remember the difference between reshape and transpose.
# Reshape data to n_voxels by n_vols, call this `data_2d`
# Then make result into n_vols by n_voxels
# Check shape is (n_vols, n_voxels)
data_2d = data.reshape((n_voxels, n_vols))
data_2d = data_2d.T
assert data_2d.shape == (n_vols, n_voxels)
Now we have our data ready for the estimation, we can compile our design matrix.
To do this, we go back over the steps we covered in the day6 regression notebok.
The TR for this run was 2.5 seconds.
# - TR in seconds
TR = 2.5
There was only one event type for this run. See the ds114 data description for details.
The experiment was an on-off block design. During "on" periods the subject saw a written noun every 3 seconds and had to think of an associated verb. During "off" periods subjects saw scrambled versions of the written nouns. The run started with 10 seconds "off" followed by 6 repeats of (30 seconds "on" followed by 30 seconds "off"). The run finshed with 30 seconds "on" and 25 seconds "off".
The event onsets in seconds and event durations are in the text file ds114_sub009_t2r1_cond.txt.
We can load this text file and return the neural prediction using our stimuli.events2neural function.
# Use IPython to check the signature of stimuli.events2neural
#: stimuli.events2neural?
The event onsets in ds114_sub009_t2r1_cond.txt refer to seconds in terms of the original run. During processing we dropped 4 scans at the beginning of the run.
So, the number of volumes we need to specify to stimuli.events2neural is the number of volumes in the original run, before we dropped these volumes.
# Create neural time course prediction with "stimuli.events2neural"
# Remember the extra volume!
# Plot your neural time just to check it is what we expect.
neural_prediction = stimuli.events2neural('ds114_sub009_t2r1_cond.txt', TR, n_vols + n_dropped)
plt.plot(neural_prediction)
We dropped the first four volumes, so we need to drop the first four elements of the predicted neural time course:
# Remove the first four elements of the neural time course
# Check the neural time course now has the same length as number of volumes
neural_prediction = neural_prediction[n_dropped:]
assert len(neural_prediction) == n_vols
Remember HRF convolution? This is revision from day5. You might want to check back there to remember how we did this last time.
We are going to use the hemodynamic response function in spm_funcs.spm_hrf:
# Use IPython to check the signature of spm_funcs.spm_hrf
#: spm_funcs.spm_hrf?
First we need a time vector giving times at which we are going to sample our HRF function. We want to sample from time 0, up to 35 seconds, sampling every 2.5 seconds (our TR):
# Make `times` vector at which to sample our HRF
times = np.arange(0, 35, TR)
Now sample the HRF at these times to get our convolution kernel:
# Sample the HRF at given times to give convolution kernel
# Plot the HRF values against times to make sure it looks OK
hrf_kernel = spm_funcs.spm_hrf(times)
plt.plot(times, hrf_kernel)
Use np.convolve to convolve the predicted neural time course with the HRF samples. This will give us a predicted hemodynamic time course.
Don't forget to drop the last length_of_kernel - 1 values that resulted from the convolution. As you remember, these values that we are trimming are the predicted values for volumes after the end of our run.
# Convolve the neural prediction with the HRF samples
# Drop the last `length_of_kernel - 1` values
# Check the result has the same length `n_vols`
# Plot the result
N = len(hrf_kernel)
hemo_prediction = np.convolve(neural_prediction, hrf_kernel)[:-(N-1)]
assert len(hemo_prediction) == n_vols
plt.plot(hemo_prediction)
We have the hemodynamic prediction for the event.
We now need something reasonable to model the drift of the signal over time.
Start with a linear drift. Use np.linspace to make a linear drift term that starts at -1 for the first volume and ends with 1 for the last volume:
# Linear drift term using `np.linspace`
# Starts at -1, ends at 1
linear_drift = np.linspace(-1, 1, n_vols)
Now compile the design matrix. We want a design matrix X with three columns:
- The hemodynamic prediction
- The linear trend regressor
- The vector of ones that will model the mean
It should be shape n_vols by 3.
# Compile design matrix X with columns:
# * The hemodynamic prediction
# * The linear trend regressor
# * The vector of ones that will model the mean
# Show with `plt.imshow`
X = np.ones((n_vols, 3))
X[:, 0] = hemo_prediction
X[:, 1] = linear_drift
plt.imshow(X, aspect='auto')
Now we can estimate the betas on our reshaped n_vols by n_voxels data, using npl.pinv:
# Estimate betas for every voxel, using pinv
# Check resulting betas array is 3 by n_voxels
betas = npl.pinv(X).dot(data_2d)
assert betas.shape == (3, n_voxels)
Define a contrast to select only the first beta. Matrix multiply (dot) your contrast by the betas array to get a 1D array with length n_voxels.
# Define contrast to select the first beta
# Take dot of contrast with betas array to get `con_data_1d` array
# (this is the top half of the t-statistic)
# Check `con_data_1d` has shape (n_voxels,)
c = np.array([1, 0, 0])
con_data_1d = c.dot(betas)
assert con_data_1d.shape == (n_voxels,)
Reshape this array to the original shape of the 3D volumes of wafds114_sub009_t2r1.nii.
Display the middle slice (over the third dimension):
# Reshape contrast result to shape `vol_shape`
# Show middle slice (over the third dimension)
con_data = con_data_1d.reshape(vol_shape)
center_slice = vol_shape[-1] // 2
plt.imshow(con_data[:, :, center_slice])
Next we want to get some idea of the error in our estimates.
Remember our t statistic formula:
$$ t = \frac{c^T \hat\beta}{\sqrt{\hat{\sigma}^2 c^T (X^T X)^+ c}} $$We get $\hat{\sigma}^2$ from the data. We will start with this estimate of error in our data.
$\hat{\sigma}^2$ can also be called the mean residual sum of squares (MRSS).
Our first task is to get the residual sum of squares (RSS):
The residual sum of squares is :
$$ \textrm{RSS} = \sum_{i=1}^{n} (Y_i - \widehat{Y_i})^2 $$$\widehat{Y}$ is the fitted data - that is - the data as predicted from the design and the estimated betas $\hat\beta$
# Calculate the fitted data from the design and the betas
# Check the fitted data is shape (n_vols, n_voxels)
fitted = X.dot(betas)
assert fitted.shape == (n_vols, n_voxels)
The residuals are the real data minus the fitted data:
# Subtract fitted data from real data to get residuals
residuals = data_2d - fitted
The RSS is sum of squared residuals. The sum is over volumes:
# Calculate the RSS
# Check this is shape (n_voxels,)
RSS = (residuals ** 2).sum(axis=0)
assert RSS.shape == (n_voxels,)
$\hat{\sigma}^2$ is the mean residual sum of squares. To get this we divide the RSS by the degrees of freedom due to error.
This is turn is given by the number of observations (n_vols) minus the number of independent columns in the design X.
Use npl.matrix_rank to calculate the degrees of freedom due to error.
# Calculate the degrees of freedom due to error
df = n_vols - npl.matrix_rank(X)
df
# Calculate the MRSS
MRSS = RSS / df
Now we have an esimate of the error in the data, let us compare this to the contrast data.
Convert the MRSS data to a 3D image, and show the middle slice over the third axis (or the same slice you showed above for the contrast data):
# Reshape MRSS result to shape `vol_shape`
# Show middle slice (over the third dimension)
MRSS_3d = MRSS.reshape(vol_shape)
plt.imshow(MRSS_3d[:, :, center_slice])
Returning to our t statistic:
$$ t = \frac{c^T \hat\beta}{\sqrt{\hat{\sigma}^2 c^T (X^T X)^+ c}} $$We have $\hat{\sigma}^2$ (the MRSS).
Now we need to calculate $c^T (X^T X)^+ c$ (the variance due to the design).
# Calculate the pseudoinverse of (the dot product of X.T and X)
# Do the left and right dot product with the c vector
# Call the result `design_variance`
iXtX = npl.pinv(X.T.dot(X))
design_variance = c.dot(iXtX).dot(c)
design_variance
Now we can calculate the t statistic.
But - before we do this - we have a small problem - which is that the data outside the brain have very small values. This will waste processing time, and lead to division by zero errors, for voxels we are not interested in.
So, first we want to make a mask to select voxels that have a reasonable signal.
To find voxels with a reasonable signal, let's first take the mean over volumes of the wafds114_sub009_t2r1.nii 4D series:
# Get mean over original 4D functional image
mean_functional = data.mean(axis=-1)
mean_functional.shape
Next plot the histogram of the (raveled, 1D version of the) mean image. Select a signal threshold that seems to separate the low signal from the high signal.
# Plot a histogram of the raveled mean image.
# Hint: remember the `bins` keyword
_ = plt.hist(mean_functional.ravel(), bins=100) # Don't output annoying bin data
Make a mask that is True for voxels that are above this threshold, False otherwise.
Show the middle slice (over the third dimension). Adjust the threshold if it doesn't look right.
# Make mask for mean image that is True for voxels over threshold
# Show middle slice (over third axis)
thresh = 700
mask = mean_functional > thresh
plt.imshow(mask[:, :, center_slice])
We are going to calculate the t-statistics within this mask.
- Find the contrast (top half of t-statistic) data within this mask;
- Find the MRSS data within this mask; (use the 3D version of MRSS)
- Use the MRSS data and
design_varianceto calculate t-statistics within the mask; - Make a 3D volume of zeros called
t_vol; - Fill the values of
t_volwithin the mask, with the t statistic values. - Show the middle slice (over the third axis).
# * Find the contrast (top half of t-statistic) data within this mask;
# * Find the MRSS data within this mask; (use the 3D version of MRSS)
# * Use the MRSS data and `design_variance` to calculate t-statistics within the mask;
# * Make a 3D volume of zeros called `t_vol`;
# * Fill the values of `t_vol` within the mask, with the t statistic values.
# * Show the middle slice (over the third axis).
t_top = con_data[mask]
MRSS_in = MRSS_3d[mask]
t_statistics = t_top / np.sqrt(MRSS_in * design_variance)
t_vol = np.zeros(vol_shape)
t_vol[mask] = t_statistics
plt.imshow(t_vol[:, :, center_slice])
Now calculate the Bonferroni p-value threshold for a family-wise error rate of 0.05 and the number of voxels within the mask:
# Calculate the Bonferroni threshold for FWE alpha 0.05, number of voxels in mask
bonf_thresh = 0.05 / np.sum(mask)
bonf_thresh
Calculate corresponding t-value using the Scipy t-distribution, and the degrees of freedom due to error. You want the inverse survival function method (isf) to get the t statistic value for a given probability value and degrees of freedom:
# - get t distribution code from scipy library
from scipy.stats import t as t_dist
# Get t statistic corresponding to Bonferroni p-threshold
t_thresh = t_dist.isf(bonf_thresh, df)
t_thresh
Make a volume that has True in voxels where the t statistic value is greater than the Bonferroni t threshold and False otherwise. Show the middle slice (over the third dimension).
# Make volume that is True where t statistic is greater than Bonferroni t thresh
# Show middle slice over third dimension
survives_bonf = t_vol > t_thresh
plt.imshow(survives_bonf[:, :, center_slice])
Save the t volume and contrast images to disk for us to play with later:
# Save the t statistic 3D volume to disk as a NIfTI image
# Save the contrast value 3D volume to disk as a NIfTI image
t_img = nib.Nifti1Image(t_vol, img.affine, img.header)
nib.save(t_img, 't_task2.nii')
con_img = nib.Nifti1Image(con_data, img.affine, img.header)
nib.save(con_img, 'con_task2.nii')
Plotting with nilearn¶
See if you can follow along using the examples at:
https://nilearn.github.io/auto_examples/manipulating_visualizing/plot_demo_plotting.html
# import plotting package from nilearn
from nilearn import plotting
# Plot the t statistic image with the glass brain display
plotting.plot_glass_brain(t_img)
# Plot the t statistic image thresholded at the Bonferroni t threshold
plotting.plot_glass_brain(t_img, threshold=t_thresh)
# Use the `plot_stat_map` command to display the t image on the default template
# Use the Bonferroni t statistic threshold
plotting.plot_stat_map(t_img, threshold = t_thresh)
