General linear model with Xibabel¶
import numpy as np
import pandas as pd
import xarray as xr
import nibabel as nib
import xibabel as xib
# For test images.
from xibabel import testing
Make sure we have the minimal test data; we will use a test image for the General Linear Model (GLM).
# Get the data
testing.get_set('minimal')
{PosixPath('/Users/mb312/.xibabel/data/ds000009/.git/annex/objects/jv/80/MD5E-s42750881--03498e47bc89b855e1e71c8b1231e338.nii.gz/MD5E-s42750881--03498e47bc89b855e1e71c8b1231e338.nii.gz'), PosixPath('/Users/mb312/.xibabel/data/ds000009/.git/annex/objects/w9/WZ/MD5E-s12474309--4997ec21f72f4e1f7b2f9fafaf3c805e.nii.gz/MD5E-s12474309--4997ec21f72f4e1f7b2f9fafaf3c805e.nii.gz'), PosixPath('/Users/mb312/.xibabel/data/ds000009/sub-07/func/sub-07_task-balloonanalogrisktask_events.tsv'), PosixPath('/Users/mb312/.xibabel/data/ds000105/.git/annex/objects/8K/Vj/MD5E-s25669537--7e1f97ca1367d33f06e12170feeabe4a.nii.gz/MD5E-s25669537--7e1f97ca1367d33f06e12170feeabe4a.nii.gz'), PosixPath('/Users/mb312/.xibabel/data/ds000105/sub-1/func/sub-1_task-objectviewing_run-01_events.tsv')}
We start by getting the machinery to compile the design matrix.
# For constructing the design.
from nipy.modalities.fmri.design import (block_design, natural_spline)
This is the root for the images and design files:
img_path_root = (testing.DATA_PATH /
'ds000105' /
'sub-1' /
'func' /
'sub-1_task-objectviewing_run-01_')
# 4D image.
bold_path = img_path_root.with_name(img_path_root.name + 'bold.nii.gz')
# Design file.
tsv_path = img_path_root.with_name(img_path_root.name + 'events.tsv')
# Load the events ready to make a design.
event_df = pd.read_csv(tsv_path, sep='\t')
df = event_df.rename(columns={'onset': 'start', 'duration': 'end'})
df['end'] = df['start'] + df['end']
df
start | end | trial_type | |
---|---|---|---|
0 | 12.0 | 12.5 | scissors |
1 | 14.0 | 14.5 | scissors |
2 | 16.0 | 16.5 | scissors |
3 | 18.0 | 18.5 | scissors |
4 | 20.0 | 20.5 | scissors |
... | ... | ... | ... |
91 | 278.0 | 278.5 | chair |
92 | 280.0 | 280.5 | chair |
93 | 282.0 | 282.5 | chair |
94 | 284.0 | 284.5 | chair |
95 | 286.0 | 286.5 | chair |
96 rows × 3 columns
Nipy needs the design as NumPy recarrays.
block_spec = df.to_records(index=None)
GLM the Nibabel way¶
We'll use Nibabel for now to get the parameter for the design.
nib_img = nib.load(bold_path)
vol_shape = nib_img.shape[:-1]
n_vols = nib_img.shape[-1]
TR = nib_img.header.get_zooms()[-1]
# Of course this array comes naturally from xirr['time'] below.
t = np.arange(n_vols) * TR
regressors, contrasts = block_design(block_spec, t)
con_tt0 = contrasts['trial_type_0']
n_contrasts = con_tt0.shape[0]
# Add drift regressors.
drift = natural_spline(t)
n_drift = drift.shape[1]
design_matrix = np.hstack([regressors, drift])
# Contrasts for the eight event types.
con_mat = np.hstack([con_tt0, np.zeros((n_contrasts, n_drift))])
# For notation!
X = design_matrix
Nibabel's incantation to get the 4D array.
# Analysis the old way.
# The 4D array (i, j, k, time)
data = nib_img.get_fdata()
We reshape the data to time by voxels to run the estimation.
# Reshape to time by (i * j * k)
vols_by_voxels = np.reshape(data, (-1, n_vols)).T
# Run estimation with B = X+ Y
# To get B such that errors (E = Y - X @ B) are minimized in the sense
# of having the smallest np.sum(E ** 2, axis=0) for each column.
# Betas for each voxel (n_cols_in_X by n_voxels).
pX = np.linalg.pinv(X)
B = pX @ vols_by_voxels
B
is a parameters by voxels array of parameters, with one row for each column in the design matrix X
.
Run the contrasts and assemble the corresponding 3D array.
# Contrast applied. Two slopes compared
c = con_mat[4, :]
con_vec = c @ B
con_arr = np.reshape(con_vec, vol_shape)
This was all rather clunky. Xarray (inside Xibabel) makes this rather easier to write.
GLM with Xibabel¶
Load the Nifti image with Xibabel. Xibabel can also load its own zarr-based format, or from the equivalent HDF5.
xib_img = xib.load(bold_path)
xib_img
<xarray.DataArray 'sub-1_task-objectviewing_run-01_bold' (i: 40, j: 64, k: 64, time: 121)> Size: 159MB dask.array<array, shape=(40, 64, 64, 121), dtype=float64, chunksize=(40, 64, 64, 121), chunktype=numpy.ndarray> Coordinates: * i (i) int64 320B 0 1 2 3 4 5 6 7 8 9 ... 31 32 33 34 35 36 37 38 39 * j (j) int64 512B 0 1 2 3 4 5 6 7 8 9 ... 55 56 57 58 59 60 61 62 63 * k (k) int64 512B 0 1 2 3 4 5 6 7 8 9 ... 55 56 57 58 59 60 61 62 63 * time (time) float64 968B 0.0 2.5 5.0 7.5 ... 292.5 295.0 297.5 300.0 Attributes: RepetitionTime: 2.5 xib-affines: {'scanner': [[-3.5, 0.0, 0.0, 68.25], [0.0, 3.75, 0.0, -...
# Make the design
xesign = xr.DataArray(pX, dims=['p', 'time'])
We may want to specify data chunking to save memory. In this case we specify we want to run the calculation with slabs of 5 slices.
# Optional: make the data chunky.
chunked = xib_img.chunk({'k': 5})
The estimation is more obvious, because of the array axis naming.
# Do the estimation
xB = xr.dot(xesign, xib_img, dim=['time', 'time'])
This is the equivalent contrast.
xC = xr.DataArray(c, dims=['p'])
Multiplying by the contrast is safer and more obvious, because of the named axes.
x_c_arr = xr.dot(xC, xB, dim=['p', 'p'])
assert np.allclose(x_c_arr, con_arr)