Design Matrices¶

DesignMatrix is a dataframe that knows it describes a timeseries. It carries a sampling frequency, tracks which of its columns are convolved task regressors and which are confounds, and offers the operations a GLM needs: HRF convolution, drift terms, run-wise concatenation, and collinearity diagnostics.

The layout is the usual machine-learning one, observations by features: TRs by conditions plus nuisance regressors for a first-level analysis, or participants by conditions for a second-level one.

Build one by hand¶

A toy design for one participant: 22 TRs at 1.5 s, four stimulus conditions, each on for 2 TRs with a TR of rest between events. sampling_freq is in hertz, so it is 1 / TR.

import numpy as np

from nltools.data import DesignMatrix

TR = 1.5

dm = DesignMatrix(
    np.array(
        [
            [0, 0, 0, 0],
            [0, 0, 0, 0],
            [1, 0, 0, 0],
            [1, 0, 0, 0],
            [0, 0, 0, 0],
            [0, 1, 0, 0],
            [0, 1, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 1, 0],
            [0, 0, 1, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 1],
            [0, 0, 0, 1],
            [0, 0, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0],
            [0, 0, 0, 0],
        ]
    ),
    columns=["face_A", "face_B", "house_A", "house_B"],
    sampling_freq=1.0 / TR,
)
dm.head()

DesignMatrix(sampling_freq=0.6666666666666666, shape=(5, 4))

Printing it shows what the object knows about itself: sampling frequency, shape, and — once there are any — which columns are convolved and which are confounds.

print(dm)

DesignMatrix(sampling_freq=0.6666666666666666, shape=(22, 4))

.plot() draws the SPM/FSL-style heatmap, TRs down and regressors across:

dm.plot()

Drift and baseline regressors¶

Legendre polynomials¶

add_poly adds orthogonal Legendre polynomials on the interval -1 to 1, the same convention other packages use. Order 2 with include_lower=True gives an intercept, a linear trend and a quadratic trend.

dm_with_nuisance = dm.add_poly(2, include_lower=True)
dm_with_nuisance.plot()

Three columns appeared, and the object now lists them as confounds:

print(dm_with_nuisance)

DesignMatrix(sampling_freq=0.6666666666666666, shape=(22, 7))
  confounds (3): ['.nl_poly_0', '.nl_poly_1', '.nl_poly_2']

Discrete cosine basis¶

add_dct_basis is the other standard choice: a bank of cosine filters acting as a high-pass filter. duration sets the cutoff period in seconds and defaults to 180 s; 20 s suits this very short toy run.

dm.add_dct_basis(duration=20).plot()

Convolution¶

convolve applies a hemodynamic response function to the task regressors and skips the confound columns, so drift terms are left alone. kernel= defaults to "glover" and takes any of nilearn's HRF models ("glover_time", "glover_dispersion", "spm", "spm_time", "spm_dispersion"); pass a 1-D array for your own kernel, or a 2-D array to convolve with several at once. Convolved columns get a _c0 suffix (_c1, _c2, … for further kernels) so you can reference them by name.

convolved = dm_with_nuisance.convolve()
print(convolved)
convolved.plot()

DesignMatrix(sampling_freq=0.6666666666666666, shape=(22, 7))
  convolved (4): ['face_A_c0', 'face_B_c0', 'house_A_c0', 'house_B_c0']
  confounds (3): ['.nl_poly_0', '.nl_poly_1', '.nl_poly_2']

DesignMatrix also upsamples, downsamples and z-scores; the DesignMatrix reference has the full list.

Read a design from files¶

From an onsets file¶

Pass a path to a 2- or 3-column onsets file and DesignMatrix builds the regressors for you. run_length is the number of TRs and TR the repetition time. The constructor convolves by default; hrf_model=None gives raw boxcars instead, which is what you want when you plan to build interaction terms first. nltools ships an example file where each event lasts 10 s.

import os

from nltools.datasets import get_resource_path

RUN_TR = 2.0
onsets_file = os.path.join(get_resource_path(), "onsets_example.csv")

onsets_dm = DesignMatrix(
    onsets_file, run_length=160, TR=RUN_TR, hrf_model=None
).add_poly(1)
onsets_dm.plot()

From any table¶

Anything polars or pandas can read becomes a DesignMatrix — you only have to supply the sampling frequency. Here is a motion-realignment file from a preprocessing pipeline. .plot() forwards keyword arguments to seaborn, so vmin/vmax rescale the color range to the size of these regressors.

import pandas as pd

covariates_file = os.path.join(get_resource_path(), "covariates_example.csv")
covariates = DesignMatrix(pd.read_csv(covariates_file), sampling_freq=1.0 / RUN_TR)
covariates.plot(vmin=-1, vmax=1)