GLM Fingertapping Example

[1]:

import cedalion
import cedalion.datasets
import cedalion.io
import cedalion.models.glm as glm
import cedalion.nirs
import cedalion.plots as plots
import cedalion.sigproc.frequency
import matplotlib.pyplot as p
import pandas as pd
import xarray as xr
from cedalion import units

xr.set_options(display_expand_data=False);

Loading and preprocessing the dataset

This notebook uses a finger-tapping dataset in BIDS layout provided by Rob Luke. It can can be downloaded via cedalion.datasets.

[2]:

rec = cedalion.datasets.get_fingertapping()

# rename trials
rec.stim.cd.rename_events(
    {
        "1.0": "control",
        "2.0": "Tapping/Left",
        "3.0": "Tapping/Right",
        "15.0": "sentinel",
    }
)
rec.stim = rec.stim[rec.stim.trial_type != "sentinel"]

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

# calculate optical density and concentrations
rec["od"] = cedalion.nirs.int2od(rec["amp"])
rec["conc"] = cedalion.nirs.od2conc(rec["od"], rec.geo3d, dpf, spectrum="prahl")

# Bandpass filter remove cardiac component and slow drifts.
# Here we use a highpass to remove drift. Another possible option would be to
# use drift regressors in the design matrix.
fmin = 0.02 * units.Hz
fmax = 0.3 * units.Hz

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

display(rec)

<Recording |  timeseries: ['amp', 'od', 'conc', 'conc_filtered'],  masks: [],  stim: ['control', 'Tapping/Left', 'Tapping/Right'],  aux_ts: [],  aux_obj: []>

Plot freq. filtered concentration data for two channels on the left (S1D1, S1D3) and right (S5D5, S5D7) hemispheres.

[3]:

ts = rec["conc_filtered"]

f, ax = p.subplots(4, 1, sharex=True, figsize=(12, 6))
for i, ch in enumerate(["S1D1", "S1D3", "S5D5", "S5D7"]):
    ax[i].plot(ts.time, ts.sel(channel=ch, chromo="HbO"), "r-", label="HbO")
    ax[i].plot(ts.time, ts.sel(channel=ch, chromo="HbR"), "b-", label="HbR")
    ax[i].set_title(f"Ch. {ch}")
    cedalion.plots.plot_stim_markers(ax[i], rec.stim, y=1)
    ax[i].set_ylabel(r"$\Delta$ c / uM")

ax[0].legend(ncol=6)
ax[3].set_label("time / s")
ax[3].set_xlim(0,300)
p.tight_layout()

../../_images/examples_modeling_32_glm_fingertapping_example_5_0.png

Build design matrix

use the glm.make_design_matrix method to build regressors
to account for signal components from superficial layers use short-distance channel regression: for each long channel the closest short channel is selected. From these the channel-wise regressor’short’ is derived.

[4]:

# split time series into two based on channel distance
ts_long, ts_short = cedalion.nirs.split_long_short_channels(
    rec["conc_filtered"], rec.geo3d, distance_threshold=1.5 * units.cm
)

# build regressors
dm, channel_wise_regressors = glm.make_design_matrix(
    ts_long,
    ts_short,
    rec.stim,
    rec.geo3d,
    basis_function=glm.Gamma(tau=0 * units.s, sigma=3 * units.s, T=3 * units.s),
    drift_order=1,
    short_channel_method="closest",
)

The design matrix dm holds all regressors that apply to all channels. It has dimensions ‘time’, ‘chromo’ and ‘regressor’. Regressors have string labels.

[5]:

display(dm)

channel_wise_regressors is list of additional xr.DataArrays that contain regressors which differ between channels. Each such array may contain only one regressor (i.e. the size of the regressor dimension must be 1). The regressors for each channel are arranged in the additional ‘channel’ dimension.

[6]:

display(channel_wise_regressors[0])

<xarray.DataArray 'concentration' (regressor: 1, chromo: 2, channel: 20,
                                   time: 23239)> Size: 7MB
[µM] 0.1092 0.1416 0.1736 0.2054 0.2366 ... 0.03124 0.02521 0.01864 0.01182
Coordinates:
  * chromo         (chromo) <U3 24B 'HbO' 'HbR'
  * time           (time) float64 186kB 0.0 0.128 0.256 ... 2.974e+03 2.974e+03
    samples        (time) int64 186kB 0 1 2 3 4 ... 23235 23236 23237 23238
  * channel        (channel) object 160B 'S1D1' 'S1D2' 'S1D3' ... 'S8D7' 'S8D8'
    short_channel  (channel) object 160B 'S1D9' 'S1D9' ... 'S8D16' 'S8D16'
    comp_group     (channel) int64 160B 0 0 0 1 1 1 2 2 3 3 4 4 4 5 5 5 6 6 7 7
  * regressor      (regressor) <U5 20B 'short'

Fitting the model

[8]:

betas = glm.fit(ts_long, dm, channel_wise_regressors, noise_model="ols")

display(betas)

#pd.set_option('display.max_rows', None)
display(betas.rename("beta").to_dataframe())

			beta
regressor	channel	chromo
HRF control	S1D1	HbO	0.033086
	S1D1	HbR	-0.017289
	S1D2	HbO	0.007820
	S1D2	HbR	-0.000802
	S1D3	HbO	0.026215
...	...	...	...
short	S7D7	HbR	0.113651
	S8D7	HbO	0.615008
	S8D7	HbR	0.226749
	S8D8	HbO	0.407106
	S8D8	HbR	-0.212301

240 rows × 1 columns

Model Predictions

using glm.predict one can scale the regressors in dm and channel_wise_regressors with the estimated coefficients to obtain a model prediction
by giving only a subset of betas to glm.predict one can predict subcomponents of the model

[9]:

# prediction using all regressors
pred = glm.predict(ts_long, betas, dm, channel_wise_regressors)

# prediction of all nuisance regressors, i.e. all regressors that don't start with 'HRF '
pred_wo_hrf = glm.predict(
    ts_long,
    betas.sel(regressor=~betas.regressor.str.startswith("HRF ")),
    dm,
    channel_wise_regressors,
)

# prediction of all HRF regressors, i.e. all regressors that start with 'HRF '
pred_hrf = glm.predict(
    ts_long,
    betas.sel(regressor=betas.regressor.str.startswith("HRF ")),
    dm,
    channel_wise_regressors,
)

Plot model predictions

[10]:

# plot the data and model prediction
#ch = "S6D7"
ch = "S1D3"
f, ax = p.subplots(1,1, figsize=(12, 4))
p.plot(ts_long.time, ts_long.sel(chromo="HbO", channel=ch), "r-", label="data HbO", alpha=.5)
p.plot(pred.time, pred.sel(chromo="HbO", channel=ch), "r-", label="model", lw=2 )
p.plot(pred.time, pred_wo_hrf.sel(chromo="HbO", channel=ch), "k:", label="model w/o HRF", alpha=.5)
plots.plot_stim_markers(ax, rec.stim, y=1)
p.xlim(60,300)
p.ylim(-.4,.4)
p.xlabel("time / s")
p.ylabel(r"$\Delta$  c / uM")
p.legend(ncol=4)


# subtract from data nuisance regressors and plot against predicted HRF components
f, ax = p.subplots(1,1, figsize=(12, 4))
p.plot(pred_hrf.time, pred_hrf.sel(chromo="HbO", channel=ch), "r-", label="HRF HbO")
p.plot(pred_hrf.time, pred_hrf.sel(chromo="HbR", channel=ch), "b-", label="HRF HbR")
p.plot(
    pred_hrf.time,
    ts_long.sel(chromo="HbO", channel=ch).pint.dequantify() - pred_wo_hrf.sel(chromo="HbO", channel=ch),
    "r-", label="data HbO - nuisance reg.", alpha=.5
)
p.plot(
    pred_hrf.time,
    ts_long.sel(chromo="HbR", channel=ch).pint.dequantify() - pred_wo_hrf.sel(chromo="HbR", channel=ch),
    "b-", label="data HbR - nuisance reg.", alpha=.5
)
plots.plot_stim_markers(ax, rec.stim, y=1)
p.legend(ncol=4, loc="lower right")

p.xlim(60,500)
p.xlabel("time / s")
p.ylabel(r"$\Delta$  c / uM");

../../_images/examples_modeling_32_glm_fingertapping_example_19_0.png

../../_images/examples_modeling_32_glm_fingertapping_example_19_1.png

Scalp plots

[11]:

f, ax = p.subplots(2, 3, figsize=(12, 8))
vlims = {"HbO" : [0.,0.3], "HbR" : [-0.1, 0.05]}
for i_chr, chromo in enumerate(betas.chromo.values):
    vmin, vmax = vlims[chromo]
    for i_reg, reg in enumerate(["HRF Tapping/Left", "HRF Tapping/Right", "HRF control"]):
        cedalion.plots.scalp_plot(
            rec["amp"],
            rec.geo3d,
            betas.sel(chromo=chromo, regressor=reg),
            ax[i_chr, i_reg],
            min_dist=1.5 * cedalion.units.cm,
            title=f"{chromo} {reg}",
            vmin=vmin,
            vmax=vmax,
            optode_labels=True,
            cmap="RdBu_r",
            cb_label=r"$\beta$"
        )
p.tight_layout()

../../_images/examples_modeling_32_glm_fingertapping_example_21_0.png