# Modeling a single voxel¶

A long time ago (Voxel time courses), we were looking at a single voxel time course.

Let’s get that same voxel time course back again:

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> import nibabel as nib

>>> img = nib.load('ds114_sub009_t2r1.nii')
>>> data = img.get_data()
>>> data = data[..., 4:]


The voxel coordinate (3D coordinate) that we were looking at before was at (42, 32, 19):

>>> voxel_time_course = data[42, 32, 19]
>>> plt.plot(voxel_time_course)
[...]


(png, hires.png, pdf) Now we are going to use our new convolved regressor to do a simple regression on this voxel time course.

If you don’t have it already, you will need to download ds114_sub009_t2r1_conv.txt.

>>> convolved = np.loadtxt('ds114_sub009_t2r1_conv.txt')
>>> # Knock off first 4 elements to match data
>>> convolved = convolved[4:]
>>> plt.plot(convolved)
[...]


(png, hires.png, pdf) First we make our design matrix. It has a column for the convolved regressor, and a column of ones:

>>> N = len(convolved)
>>> X = np.ones((N, 2))
>>> X[:, 0] = convolved
>>> plt.imshow(X, interpolation='nearest', cmap='gray', aspect=0.1)
<...>


(png, hires.png, pdf) $$\newcommand{\yvec}{\vec{y}}$$ $$\newcommand{\xvec}{\vec{x}}$$ $$\newcommand{\evec}{\vec{\varepsilon}}$$ $$\newcommand{Xmat}{\boldsymbol X} \newcommand{\bvec}{\vec{\beta}}$$ $$\newcommand{\bhat}{\hat{\bvec}} \newcommand{\yhat}{\hat{\yvec}}$$

As you will remember from Introduction to the general linear model, our model is:

$\yvec = \Xmat \bvec + \evec$

We can get our least squares parameter estimates for $$\bvec$$ with:

$\bhat = \Xmat^+y$

where $$\Xmat^+$$ is the pseudoinverse of $$\Xmat$$. When $$\Xmat$$ is invertible, the pseudoinverse is given by:

$\Xmat^+ = (\Xmat^T \Xmat)^{-1} \Xmat^T$

Let’s calculate the pseudoinverse for our design:

>>> import numpy.linalg as npl
>>> Xp = npl.pinv(X)
>>> Xp.shape
(2, 169)


We calculate $$\bhat$$:

>>> beta_hat = Xp.dot(voxel_time_course)
>>> beta_hat
array([   31.185514,  2029.367685])


We can then calculate $$\yhat$$ (also called the fitted data):

>>> y_hat = X.dot(beta_hat)
>>> e_vec = voxel_time_course - y_hat
>>> print(np.sum(e_vec ** 2))
41405.5727762
>>> plt.plot(voxel_time_course)
[...]
>>> plt.plot(y_hat)
[...]


(png, hires.png, pdf) 