{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# S4: Model-driven (GLM) Analysis\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 1, "id": "1", "metadata": { "execution": { "iopub.execute_input": "2026-01-13T08:14:12.618683Z", "iopub.status.busy": "2026-01-13T08:14:12.618437Z", "iopub.status.idle": "2026-01-13T08:14:12.625730Z", "shell.execute_reply": "2026-01-13T08:14:12.625188Z" } }, "outputs": [], "source": [ "# This cells setups the environment when executed in Google Colab.\n", "try:\n", " import google.colab\n", " !curl -s https://raw.githubusercontent.com/ibs-lab/cedalion/dev/scripts/colab_setup.py -o colab_setup.py\n", " # Select branch with --branch \"branch name\" (default is \"dev\")\n", " %run colab_setup.py\n", "except ImportError:\n", " pass" ] }, { "cell_type": "code", "execution_count": 2, "id": "2", "metadata": { "execution": { "iopub.execute_input": "2026-01-13T08:14:12.627205Z", "iopub.status.busy": "2026-01-13T08:14:12.627028Z", "iopub.status.idle": "2026-01-13T08:14:14.725941Z", "shell.execute_reply": "2026-01-13T08:14:14.725217Z" } }, "outputs": [], "source": [ "import cedalion\n", "import cedalion.sigproc.quality as quality\n", "import cedalion.sigproc.motion as motion\n", "import cedalion.sigproc.physio as physio\n", "\n", "\n", "from cedalion import units\n", "import cedalion.xrutils as xrutils\n", "import cedalion.models.glm as glm\n", "import cedalion.data\n", "\n", "import cedalion.vis.blocks as vbx\n", "from cedalion.vis.anatomy import scalp_plot, plot_montage3D\n", "from cedalion.vis.colors import p_values_cmap\n", "\n", "import numpy as np\n", "import xarray as xr\n", "\n", "import matplotlib.pyplot as p\n", "from statsmodels.stats.multitest import multipletests\n", "\n", "xr.set_options(display_expand_data=False)\n", "xrutils.unit_stripping_is_error()" ] }, { "cell_type": "markdown", "id": "3", "metadata": {}, "source": [ "## Load Data\n", "\n", "The dataset is loaded and events are renamed. Afterwards, the stimulus data frame is summarized." ] }, { "cell_type": "code", "execution_count": 3, "id": "4", "metadata": { "execution": { "iopub.execute_input": "2026-01-13T08:14:14.728875Z", "iopub.status.busy": "2026-01-13T08:14:14.728356Z", "iopub.status.idle": "2026-01-13T08:14:16.022187Z", "shell.execute_reply": "2026-01-13T08:14:16.021703Z" } }, "outputs": [], "source": [ "rec = cedalion.data.get_fingertappingDOT()\n", "\n", "# assign better trial_type labels\n", "rec.stim.cd.rename_events(\n", " {\n", " \"1\": \"Control\",\n", " \"2\": \"FTapping/Left\",\n", " \"3\": \"FTapping/Right\",\n", " \"4\": \"BallSqueezing/Left\",\n", " \"5\": \"BallSqueezing/Right\",\n", " }\n", ")" ] }, { "cell_type": "code", "execution_count": 4, "id": "5", "metadata": { "execution": { "iopub.execute_input": "2026-01-13T08:14:16.023795Z", "iopub.status.busy": "2026-01-13T08:14:16.023644Z", "iopub.status.idle": "2026-01-13T08:14:16.032582Z", "shell.execute_reply": "2026-01-13T08:14:16.032232Z" } }, "outputs": [ { "data": { "text/html": [ "
| \n", " | # trials | \n", "
|---|---|
| trial_type | \n", "\n", " |
| BallSqueezing/Left | \n", "17 | \n", "
| BallSqueezing/Right | \n", "16 | \n", "
| Control | \n", "65 | \n", "
| FTapping/Left | \n", "16 | \n", "
| FTapping/Right | \n", "16 | \n", "
| \n", " | # trials | \n", "
|---|---|
| trial_type | \n", "\n", " |
| Control | \n", "10 | \n", "
| FTapping/Left | \n", "3 | \n", "
| FTapping/Right | \n", "3 | \n", "
<xarray.DataArray 'concentration' (chromo: 2, channel: 98, time: 1352)> Size: 2MB\n",
"[µM] 0.8276 0.7548 0.7043 0.6868 0.6978 ... -0.08036 -0.07753 -0.07426 -0.07086\n",
"Coordinates:\n",
" * chromo (chromo) <U3 24B 'HbO' 'HbR'\n",
" * time (time) float64 11kB 5.046 5.276 5.505 5.734 ... 314.5 314.7 314.9\n",
" samples (time) int64 11kB 22 23 24 25 26 27 ... 1369 1370 1371 1372 1373\n",
" * channel (channel) object 784B 'S1D1' 'S1D2' 'S1D4' ... 'S14D31' 'S14D32'\n",
" source (channel) object 784B 'S1' 'S1' 'S1' 'S1' ... 'S14' 'S14' 'S14'\n",
" detector (channel) object 784B 'D1' 'D2' 'D4' 'D5' ... 'D29' 'D31' 'D32'<xarray.DataArray 'long channels' (chromo: 2, channel: 44, time: 1352)> Size: 952kB\n",
"[µM] 0.3715 0.3284 0.2975 0.2849 0.2892 ... 0.006362 0.008495 0.01122 0.01438\n",
"Coordinates:\n",
" * chromo (chromo) <U3 24B 'HbO' 'HbR'\n",
" * time (time) float64 11kB 5.046 5.276 5.505 5.734 ... 314.5 314.7 314.9\n",
" samples (time) int64 11kB 22 23 24 25 26 27 ... 1369 1370 1371 1372 1373\n",
" * channel (channel) object 352B 'S1D6' 'S1D8' 'S2D5' ... 'S14D25' 'S14D27'\n",
" source (channel) object 352B 'S1' 'S1' 'S2' 'S2' ... 'S13' 'S14' 'S14'\n",
" detector (channel) object 352B 'D6' 'D8' 'D5' 'D9' ... 'D28' 'D25' 'D27'<xarray.DataArray 'short channels' (chromo: 2, channel: 54, time: 1352)> Size: 1MB\n",
"[µM] 0.8276 0.7548 0.7043 0.6868 0.6978 ... -0.08036 -0.07753 -0.07426 -0.07086\n",
"Coordinates:\n",
" * chromo (chromo) <U3 24B 'HbO' 'HbR'\n",
" * time (time) float64 11kB 5.046 5.276 5.505 5.734 ... 314.5 314.7 314.9\n",
" samples (time) int64 11kB 22 23 24 25 26 27 ... 1369 1370 1371 1372 1373\n",
" * channel (channel) object 432B 'S1D1' 'S1D2' 'S1D4' ... 'S14D31' 'S14D32'\n",
" source (channel) object 432B 'S1' 'S1' 'S1' 'S1' ... 'S14' 'S14' 'S14'\n",
" detector (channel) object 432B 'D1' 'D2' 'D4' 'D5' ... 'D29' 'D31' 'D32'"
],
"text/plain": [
"<xarray.DataArray (time: 1352, regressor: 46, chromo: 2)> Size: 995kB\n",
"0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.9999 0.9999 -0.9999 -0.9999 0.4137 0.03177\n",
"Coordinates:\n",
" * time (time) float64 11kB 5.046 5.276 5.505 5.734 ... 314.5 314.7 314.9\n",
" * regressor (regressor) <U21 4kB 'HRF Control 00' ... 'short'\n",
" * chromo (chromo) <U3 24B 'HbO' 'HbR'\n",
" samples (time) int64 11kB 22 23 24 25 26 27 ... 1369 1370 1371 1372 1373<xarray.DataArray (time: 1352, regressor: 46, chromo: 2)> Size: 995kB\n",
"0.0 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.9999 0.9999 -0.9999 -0.9999 0.4137 0.03177\n",
"Coordinates:\n",
" * time (time) float64 11kB 5.046 5.276 5.505 5.734 ... 314.5 314.7 314.9\n",
" * regressor (regressor) <U21 4kB 'HRF Control 00' ... 'short'\n",
" * chromo (chromo) <U3 24B 'HbO' 'HbR'\n",
" samples (time) int64 11kB 22 23 24 25 26 27 ... 1369 1370 1371 1372 1373"
],
"text/plain": [
"<xarray.DataArray (channel: 44, chromo: 2)> Size: 704B\n",
"<statsmodels.robust.robust_linear_model.RLMResultsWrapper object at 0x76c0524...\n",
"Coordinates:\n",
" * chromo (chromo) <U3 24B 'HbO' 'HbR'\n",
" * channel (channel) object 352B 'S1D6' 'S1D8' 'S2D5' ... 'S14D25' 'S14D27'\n",
" source (channel) object 352B 'S1' 'S1' 'S2' 'S2' ... 'S13' 'S14' 'S14'\n",
" detector (channel) object 352B 'D6' 'D8' 'D5' 'D9' ... 'D28' 'D25' 'D27'\n",
"Attributes:\n",
" description: AR_IRLS<xarray.DataArray (channel: 44, chromo: 2, regressor: 46)> Size: 32kB\n",
"0.2767 -1.323 3.382 -5.959 7.91 ... -0.008269 -0.0007833 0.01253 0.01076 0.2799\n",
"Coordinates:\n",
" * regressor (regressor) object 368B 'HRF Control 00' ... 'short'\n",
" * chromo (chromo) <U3 24B 'HbO' 'HbR'\n",
" * channel (channel) object 352B 'S1D6' 'S1D8' 'S2D5' ... 'S14D25' 'S14D27'\n",
" source (channel) object 352B 'S1' 'S1' 'S2' 'S2' ... 'S13' 'S14' 'S14'\n",
" detector (channel) object 352B 'D6' 'D8' 'D5' 'D9' ... 'D28' 'D25' 'D27'\n",
"Attributes:\n",
" description: AR_IRLS"
],
"text/plain": [
"<xarray.DataArray (channel: 44, chromo: 2)> Size: 704B\n",
" Test for Constraints ...\n",
"Coordinates:\n",
" * chromo (chromo) <U3 24B 'HbO' 'HbR'\n",
" * channel (channel) object 352B 'S1D6' 'S1D8' 'S2D5' ... 'S14D25' 'S14D27'\n",
" source (channel) object 352B 'S1' 'S1' 'S2' 'S2' ... 'S13' 'S14' 'S14'\n",
" detector (channel) object 352B 'D6' 'D8' 'D5' 'D9' ... 'D28' 'D25' 'D27'\n",
"Attributes:\n",
" description: AR_IRLS<xarray.DataArray (time: 77, regressor: 33, chromo: 2)> Size: 41kB\n",
"1.0 1.0 0.8694 0.8694 0.5698 0.5698 ... 0.1581 0.3804 0.3804 0.6912 0.6912\n",
"Coordinates:\n",
" * time (time) float64 616B -2.294 -2.064 -1.835 ... 14.68 14.91 15.14\n",
" * regressor (regressor) <U21 3kB 'HRF Control 00' ... 'HRF FTapping/Right 10'\n",
" * chromo (chromo) <U3 24B 'HbO' 'HbR'