GLM Fingertapping Example

[1]:

import matplotlib.pyplot as p
import numpy as np
import pandas as pd
import xarray as xr

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
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
fmax = 0 * 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
)

dms = (
    glm.design_matrix.hrf_regressors(
        ts_long, rec.stim, glm.Gamma(tau=0 * units.s, sigma=3 * units.s, T=3 * units.s)
    )
    & glm.design_matrix.drift_regressors(ts, drift_order=1)
    & glm.design_matrix.closest_short_channel_regressor(ts_long, ts_short, rec.geo3d)
)

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(dms)
display(dms.common)

DesignMatrix(universal=['HRF control','HRF Tapping/Left','HRF Tapping/Right','Drift 0','Drift 1'], channel_wise=['short'])

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(dms.channel_wise[0])

<xarray.DataArray 'concentration' (regressor: 1, chromo: 2, channel: 20,
                                   time: 23239)> Size: 7MB
[µM] -0.1671 -0.3522 -0.03758 0.1025 ... 0.06639 -0.001217 0.005855 -0.001648
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'

[7]:

dms.channel_wise[0] = dms.channel_wise[0].pint.dequantify()
dms.channel_wise[0] /= dms.channel_wise[0].max("time")

Fitting the model

[9]:

results = glm.fit(ts_long, dms, noise_model="ar_irls", max_jobs=1)

display(results)

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

[10]:

betas = results.sm.params
betas

[11]:

betas.rename("betas").to_dataframe()

[11]:

			source	detector	betas
channel	chromo	regressor
S1D1	HbO	index	S1	D1	1.986168e-07
		HRF control	S1	D1	1.244468e-02
		HRF Tapping/Left	S1	D1	2.439477e-01
		HRF Tapping/Right	S1	D1	3.886186e-01
		Drift 0	S1	D1	-2.321850e-02
...	...	...	...	...	...
S8D8	HbR	HRF Tapping/Left	S8	D8	-4.158567e-02
		HRF Tapping/Right	S8	D8	-2.980406e-02
		Drift 0	S8	D8	2.812522e-02
		Drift 1	S8	D8	1.210228e-06
		short	S8	D8	9.317619e-02

280 rows × 3 columns

[12]:

# best fit parameters + confidence intervals
pd.concat([results[0,0].item().conf_int(), results[0,0].item().params.rename("beta")], axis=1)

[12]:

	0	1	beta
index	-0.000004	3.936040e-06	1.986168e-07
HRF control	-0.030763	5.565234e-02	1.244468e-02
HRF Tapping/Left	0.200736	2.871596e-01	2.439477e-01
HRF Tapping/Right	0.345392	4.318451e-01	3.886186e-01
Drift 0	-0.025848	-2.058892e-02	-2.321850e-02
Drift 1	-0.000001	-8.860806e-07	-9.991095e-07
short	0.226738	2.462663e-01	2.365021e-01

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

[13]:

# prediction using all regressors
betas = results.sm.params
pred = glm.predict(ts_long, betas, dms)#, 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 ")),
    dms,
)

# 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 ")),
    dms,
)