S4: Model-driven (GLM) Analysis

This tutorial demonstrates how to use a General Linear Model (GLM) to model the recorded time series as a superposition of hemodynamic responses and nuisance effects.

[1]:

# This cells setups the environment when executed in Google Colab.
try:
    import google.colab
    !curl -s https://raw.githubusercontent.com/ibs-lab/cedalion/dev/scripts/colab_setup.py -o colab_setup.py
    # Select branch with --branch "branch name" (default is "dev")
    %run colab_setup.py
except ImportError:
    pass

[2]:

import cedalion
import cedalion.sigproc.quality as quality
import cedalion.sigproc.motion as motion
import cedalion.sigproc.physio as physio


from cedalion import units
import cedalion.xrutils as xrutils
import cedalion.models.glm as glm
import cedalion.data

import cedalion.vis.blocks as vbx
from cedalion.vis.anatomy import scalp_plot, plot_montage3D
from cedalion.vis.colors import p_values_cmap

import numpy as np
import xarray as xr

import matplotlib.pyplot as p
from statsmodels.stats.multitest import multipletests

xr.set_options(display_expand_data=False)
xrutils.unit_stripping_is_error()

Load Data

The dataset is loaded and events are renamed. Afterwards, the stimulus data frame is summarized.

[3]:

rec = cedalion.data.get_fingertappingDOT()

# assign better trial_type labels
rec.stim.cd.rename_events(
    {
        "1": "Control",
        "2": "FTapping/Left",
        "3": "FTapping/Right",
        "4": "BallSqueezing/Left",
        "5": "BallSqueezing/Right",
    }
)

[4]:

# count trials
rec.stim.groupby("trial_type")[["trial_type"]].count().rename(
    {"trial_type": "# trials"},
    axis=1,
)

[4]:

	# trials
trial_type
BallSqueezing/Left	17
BallSqueezing/Right	16
Control	65
FTapping/Left	16
FTapping/Right	16

Trim dataset

Trim the dataset to roughly 5 minutes to keep computation times low while still preserving three trials per condition. Also exclude the ‘BallSqueezing’ conditions.

[5]:

tmin = 5
tmax = 315

rec.stim = rec.stim[
    (tmin <= rec.stim.onset)
    & (rec.stim.onset <= tmax)
    & (rec.stim.trial_type.isin(["Control", "FTapping/Left", "FTapping/Right"]))
]
# use a slice to select a time range
rec["amp_cropped"] = rec["amp"].sel(time=slice(tmin,tmax))

# count remaining trials
rec.stim.groupby("trial_type")[["trial_type"]].count().rename(
    {"trial_type": "# trials"},
    axis=1,
)

[5]:

	# trials
trial_type
Control	10
FTapping/Left	3
FTapping/Right	3

Preprocessing

The same motion-correction methods as described in the third tutorial notebook are applied, and the bad channels identified there are removed. The cleaned optical-density time series are then converted to concentration changes. Finally, a low-pass filter with a 0.5 Hz cutoff is applied to remove the cardiac component.

[6]:

rec["od"] = cedalion.nirs.cw.int2od(rec["amp_cropped"])
rec["od_tddr"] = motion.tddr(rec["od"])
rec["od_wavelet"] = motion.wavelet(rec["od_tddr"])

bad_channels = ['S13D26', 'S14D28']

rec["od_clean"] = rec["od_wavelet"].sel(channel=~rec["od"].channel.isin(bad_channels))

# differential pathlength factors
dpf = xr.DataArray(
    [6, 6],
    dims="wavelength",
    coords={"wavelength": rec["amp"].wavelength},
)

rec["conc"] = cedalion.nirs.cw.od2conc(rec["od_clean"], rec.geo3d, dpf, spectrum="prahl")

# Here we use a lowpass-filter to remove the cardiac component.
# Drift will be modeled in the design matrix.

fmin = 0 * units.Hz
fmax = 0.5 * units.Hz

rec["conc_filtered"] = cedalion.sigproc.frequency.freq_filter(rec["conc"], fmin, fmax)

TS_NAME = "conc_filtered"

[7]:

display(rec[TS_NAME])

<xarray.DataArray 'concentration' (chromo: 2, channel: 98, time: 1352)> Size: 2MB
[µM] 0.8276 0.7548 0.7043 0.6868 0.6978 ... -0.08036 -0.07753 -0.07426 -0.07086
Coordinates:
  * chromo    (chromo) <U3 24B 'HbO' 'HbR'
  * time      (time) float64 11kB 5.046 5.276 5.505 5.734 ... 314.5 314.7 314.9
    samples   (time) int64 11kB 22 23 24 25 26 27 ... 1369 1370 1371 1372 1373
  * channel   (channel) object 784B 'S1D1' 'S1D2' 'S1D4' ... 'S14D31' 'S14D32'
    source    (channel) object 784B 'S1' 'S1' 'S1' 'S1' ... 'S14' 'S14' 'S14'
    detector  (channel) object 784B 'D1' 'D2' 'D4' 'D5' ... 'D29' 'D31' 'D32'

Montage and Channel Distances

The hexagonal montage has longer (3-3.5cm) and shorter (~1.7-2.2cm) channels. We define a cut-off at 22.5 mm to separate long and short channels.

[8]:

plot_montage3D(rec["amp"], rec.geo3d)

../../_images/examples_tutorial_4_model_driven_analysis_12_0.png

Calculate and histogram channel distances:

[9]:

distances = cedalion.nirs.channel_distances(rec["amp"], rec.geo3d)

p.figure(figsize=(8,4))
p.hist(distances, 40)
p.axvline(22.5, c="r", ls="--")
p.xlabel("channel distance / mm")
p.ylabel("channel count");

../../_images/examples_tutorial_4_model_driven_analysis_14_0.png

Split long and short channels

[10]:

rec["ts_long"], rec["ts_short"] = cedalion.nirs.split_long_short_channels(
    rec[TS_NAME], rec.geo3d, distance_threshold=22.5 * units.mm
)

# This renaming is for display purposes only.
display(rec["ts_long"].rename("long channels"))
display(rec["ts_short"].rename("short channels"))

<xarray.DataArray 'long channels' (chromo: 2, channel: 44, time: 1352)> Size: 952kB
[µM] 0.3715 0.3284 0.2975 0.2849 0.2892 ... 0.006362 0.008495 0.01122 0.01438
Coordinates:
  * chromo    (chromo) <U3 24B 'HbO' 'HbR'
  * time      (time) float64 11kB 5.046 5.276 5.505 5.734 ... 314.5 314.7 314.9
    samples   (time) int64 11kB 22 23 24 25 26 27 ... 1369 1370 1371 1372 1373
  * channel   (channel) object 352B 'S1D6' 'S1D8' 'S2D5' ... 'S14D25' 'S14D27'
    source    (channel) object 352B 'S1' 'S1' 'S2' 'S2' ... 'S13' 'S14' 'S14'
    detector  (channel) object 352B 'D6' 'D8' 'D5' 'D9' ... 'D28' 'D25' 'D27'

<xarray.DataArray 'short channels' (chromo: 2, channel: 54, time: 1352)> Size: 1MB
[µM] 0.8276 0.7548 0.7043 0.6868 0.6978 ... -0.08036 -0.07753 -0.07426 -0.07086
Coordinates:
  * chromo    (chromo) <U3 24B 'HbO' 'HbR'
  * time      (time) float64 11kB 5.046 5.276 5.505 5.734 ... 314.5 314.7 314.9
    samples   (time) int64 11kB 22 23 24 25 26 27 ... 1369 1370 1371 1372 1373
  * channel   (channel) object 432B 'S1D1' 'S1D2' 'S1D4' ... 'S14D31' 'S14D32'
    source    (channel) object 432B 'S1' 'S1' 'S1' 'S1' ... 'S14' 'S14' 'S14'
    detector  (channel) object 432B 'D1' 'D2' 'D4' 'D5' ... 'D29' 'D31' 'D32'