General linear model with Xibabel¶

In [1]:

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

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

In [2]:

# Get the data
testing.get_set('minimal')

# Get the data testing.get_set('minimal')

Out[2]:

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

In [3]:

# For constructing the design.
from nipy.modalities.fmri.design import (block_design, natural_spline)

# For constructing the design. from nipy.modalities.fmri.design import (block_design, natural_spline)

This is the root for the images and design files:

In [4]:

img_path_root = (testing.DATA_PATH /
                 'ds000105' /
                 'sub-1' /
                 'func' /
                 'sub-1_task-objectviewing_run-01_')

img_path_root = (testing.DATA_PATH / 'ds000105' / 'sub-1' / 'func' / 'sub-1_task-objectviewing_run-01_')

In [5]:

# 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')

# 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')

In [6]:

# 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

# 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

Out[6]:

	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.

In [7]:

block_spec = df.to_records(index=None)

block_spec = df.to_records(index=None)

GLM the Nibabel way¶

We'll use Nibabel for now to get the parameter for the design.

In [8]:

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

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

In [9]:

# For notation!
X = design_matrix

# For notation! X = design_matrix

Nibabel's incantation to get the 4D array.

In [10]:

# Analysis the old way.
# The 4D array (i, j, k, time)
data = nib_img.get_fdata()

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

In [11]:

# Reshape to time by (i * j * k)
vols_by_voxels = np.reshape(data, (-1, n_vols)).T

# Reshape to time by (i * j * k) vols_by_voxels = np.reshape(data, (-1, n_vols)).T

In [12]:

# 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

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

In [13]:

# Contrast applied.  Two slopes compared
c = con_mat[4, :]
con_vec = c @ B

# Contrast applied. Two slopes compared c = con_mat[4, :] con_vec = c @ B

In [14]:

con_arr = np.reshape(con_vec, vol_shape)

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.

In [15]:

xib_img = xib.load(bold_path)
xib_img

xib_img = xib.load(bold_path) xib_img

Out[15]:

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

xarray.DataArray

'sub-1_task-objectviewing_run-01_bold'

• i: 40
• j: 64
• k: 64
• time: 121

• dask.array<chunksize=(40, 64, 64, 121), meta=np.ndarray>

  Array Chunk Bytes 151.25 MiB 151.25 MiB Shape (40, 64, 64, 121) (40, 64, 64, 121) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

• Coordinates: (4)
  • i
    (i)
    int64
    0 1 2 3 4 5 6 ... 34 35 36 37 38 39

    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
           18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
           36, 37, 38, 39])

  • j
    (j)
    int64
    0 1 2 3 4 5 6 ... 58 59 60 61 62 63

    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
           18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
           36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
           54, 55, 56, 57, 58, 59, 60, 61, 62, 63])

  • k
    (k)
    int64
    0 1 2 3 4 5 6 ... 58 59 60 61 62 63

    array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
           18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
           36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
           54, 55, 56, 57, 58, 59, 60, 61, 62, 63])

  • time
    (time)
    float64
    0.0 2.5 5.0 ... 295.0 297.5 300.0

    units :

    s

    array([  0. ,   2.5,   5. ,   7.5,  10. ,  12.5,  15. ,  17.5,  20. ,  22.5,
            25. ,  27.5,  30. ,  32.5,  35. ,  37.5,  40. ,  42.5,  45. ,  47.5,
            50. ,  52.5,  55. ,  57.5,  60. ,  62.5,  65. ,  67.5,  70. ,  72.5,
            75. ,  77.5,  80. ,  82.5,  85. ,  87.5,  90. ,  92.5,  95. ,  97.5,
           100. , 102.5, 105. , 107.5, 110. , 112.5, 115. , 117.5, 120. , 122.5,
           125. , 127.5, 130. , 132.5, 135. , 137.5, 140. , 142.5, 145. , 147.5,
           150. , 152.5, 155. , 157.5, 160. , 162.5, 165. , 167.5, 170. , 172.5,
           175. , 177.5, 180. , 182.5, 185. , 187.5, 190. , 192.5, 195. , 197.5,
           200. , 202.5, 205. , 207.5, 210. , 212.5, 215. , 217.5, 220. , 222.5,
           225. , 227.5, 230. , 232.5, 235. , 237.5, 240. , 242.5, 245. , 247.5,
           250. , 252.5, 255. , 257.5, 260. , 262.5, 265. , 267.5, 270. , 272.5,
           275. , 277.5, 280. , 282.5, 285. , 287.5, 290. , 292.5, 295. , 297.5,
           300. ])

• Indexes: (4)
  • i
    PandasIndex

    PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
           18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
           36, 37, 38, 39],
          dtype='int64', name='i'))

  • j
    PandasIndex

    PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
           18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
           36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
           54, 55, 56, 57, 58, 59, 60, 61, 62, 63],
          dtype='int64', name='j'))

  • k
    PandasIndex

    PandasIndex(Index([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
           18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
           36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53,
           54, 55, 56, 57, 58, 59, 60, 61, 62, 63],
          dtype='int64', name='k'))

  • time
    PandasIndex

    PandasIndex(Index([  0.0,   2.5,   5.0,   7.5,  10.0,  12.5,  15.0,  17.5,  20.0,  22.5,
           ...
           277.5, 280.0, 282.5, 285.0, 287.5, 290.0, 292.5, 295.0, 297.5, 300.0],
          dtype='float64', name='time', length=121))

• Attributes: (2)

  RepetitionTime :

  2.5

  xib-affines :

  {'scanner': [[-3.5, 0.0, 0.0, 68.25], [0.0, 3.75, 0.0, -118.125], [0.0, 0.0, 3.75, -118.125], [0.0, 0.0, 0.0, 1.0]]}

In [16]:

# Make the design
xesign = xr.DataArray(pX, dims=['p', 'time'])

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

In [17]:

# Optional: make the data chunky.
chunked = xib_img.chunk({'k': 5})

# Optional: make the data chunky. chunked = xib_img.chunk({'k': 5})

The estimation is more obvious, because of the array axis naming.

In [18]:

# Do the estimation
xB = xr.dot(xesign, xib_img, dim=['time', 'time'])

# Do the estimation xB = xr.dot(xesign, xib_img, dim=['time', 'time'])

This is the equivalent contrast.

In [19]:

xC = xr.DataArray(c, dims=['p'])

xC = xr.DataArray(c, dims=['p'])

Multiplying by the contrast is safer and more obvious, because of the named axes.

In [20]:

x_c_arr = xr.dot(xC, xB, dim=['p', 'p'])

x_c_arr = xr.dot(xC, xB, dim=['p', 'p'])

In [21]:

assert np.allclose(x_c_arr, con_arr)

assert np.allclose(x_c_arr, con_arr)