Here we try the simplest possible general linear model on some example data
%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
You are used to seeing these now:
import numpy as np
import matplotlib.pyplot as plt
Set defaults for the figures
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'
Our faithful image library
import nibabel as nib
Let's get the data:
import os
HOME = os.path.expanduser('~')
DATA_PATH = os.path.join(HOME, 'data')
img_fname = os.path.join(DATA_PATH, 'ds105', 'sub001', 'BOLD', 'task001_run001', 'bold.nii.gz')
img = nib.load(img_fname)
n_vols = img.shape[-1]
n_vols
121
Now let's get the task. There were blocks of different types of objects in an on-off pattern. We'll pool the objects to get just 'on-off':
# from sub001/model/model001/onsets/task001_run001
block_starts_seconds = np.array([12, 48, 84, 120, 156, 192, 228, 264])
block_length_seconds = 22
TR = 2.5
Same in scans:
block_starts_scans = block_starts_seconds / TR
block_length_scans = block_length_seconds / TR
block_starts_scans
array([ 4.8, 19.2, 33.6, 48. , 62.4, 76.8, 91.2, 105.6])
Make a regressor. Let's give the HRF a delay of 2 scans or so.
delay = 2
regressor = np.zeros((n_vols,))
for start in block_starts_scans:
start_tr = round(start + delay)
end_tr = round(start + delay + block_length_scans)
regressor[start_tr:end_tr+1] = 1
print start_tr, end_tr
plt.plot(regressor)
7.0 16.0 21.0 30.0 36.0 44.0 50.0 59.0 64.0 73.0 79.0 88.0 93.0 102.0 108.0 116.0
[<matplotlib.lines.Line2D at 0x4022b10>]
Make the design:
X = np.column_stack((regressor, np.ones((n_vols,))))
data = img.get_data()
slice_size = np.prod(data.shape[:-1])
print(slice_size)
Y = np.reshape(data, (slice_size, n_vols)).T
Y.shape
163840
(121, 163840)
X
array([[ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 1., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.], [ 0., 1.]])
Fit the model!
B = np.linalg.pinv(X).dot(Y)
B.shape
(2, 163840)
Make B into an image again:
B = np.reshape(B, (2,) + img.shape[:-1])
B.shape
(2, 40, 64, 64)
Let's have a look at the image:
plt.imshow(B[0, :, :, 40])
<matplotlib.image.AxesImage at 0x4534c10>
B_selected = B[:, 19, 12, 40]
B_selected
array([ 14.78175313, 540.37209302])
Yselected = data[19, 12, 40]
Yselected.shape
(121,)
Bselected = np.linalg.pinv(X).dot(Yselected)
Bselected
array([ 14.78175313, 540.37209302])
plt.plot(Yselected - np.mean(Yselected))
plt.plot(regressor * 40, 'r')
[<matplotlib.lines.Line2D at 0x40408d0>]
What is the 'fitted' time course?
Yhat = X.dot(Bselected)
plt.plot(Yhat, 'r')
plt.plot(Yselected)
[<matplotlib.lines.Line2D at 0x486a310>]
For convenience, let's package the model fit up into a function:
def fit_model(data, X):
# Shape of design
n_vols, n_params = X.shape
assert n_vols == data.shape[-1]
# Reshape data to be time by voxels
Y = np.reshape(data, (-1, n_vols)).T
B = np.linalg.pinv(X).dot(Y)
# Reshape betas to image
B = np.reshape(B, (n_params,) + data.shape[:-1])
return B
How would we we remove these drifts? We need to add something to $X$.
# What should we add to X?
B_better = fit_model(data, X_better)
B_better_selected = B_better[:, 19, 12, 40]
Yhat_better = X_better.dot(B_better_selected)
plt.plot(Yhat_better, 'r')
plt.plot(Yselected)
[<matplotlib.lines.Line2D at 0x50e7c50>]
B_better.shape
(5, 40, 64, 64)
plt.imshow(B_better[0, :, :, 40])
<matplotlib.image.AxesImage at 0x530e990>
Can we do better than this? Time for some hemodynamic modeling!
import nipy.modalities.fmri.hrf as hrfs
Tab completion right?
What does an HRF look like?
t = np.linspace(0, 30, 100)
plt.plot(t, hrfs.glovert(t))
plt.xlabel('time (seconds)')
<matplotlib.text.Text at 0x6dc7510>
How do we apply that to our data? Convolution!