Fitting a GLM with Gaussian Kernels

[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_correct as motion_correct
import cedalion.sigproc.physio as physio


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

import cedalion.plots
from cedalion.plots import segmented_cmap

from pathlib import Path
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

[3]:

rec = cedalion.datasets.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

Reduce the lenght of the dataset to about 5 minutes. This keeps computing times low for presentation and maintains 3 trials for each condition. Also ignore 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"]))
]
rec["amp_cropped"] = rec["amp"].sel(time=slice(tmin,tmax))

# count 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

apply TDDR and wavelet corrections
remove two known bad channels
transform to concentrations
apply a frequency filter

[6]:

rec["od"] = cedalion.nirs.int2od(rec["amp_cropped"])
rec["od_tddr"] = motion_correct.tddr(rec["od"])
rec["od_wavelet"] = motion_correct.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.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

[8]:

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

../../_images/examples_modeling_36_glm_workshop_12_0.png

[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_modeling_36_glm_workshop_13_0.png

The montage has longer (3-3.5cm) and shorter (~1.7-2.2cm) distance channels. Define a cut-off at 22.5 mm to define 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
)

display(rec["ts_long"])
display(rec["ts_short"])

<xarray.DataArray 'concentration' (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 'concentration' (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'