Multiple Predictor Features and EEG Channels¶

Use a multivariate forward model when the scientific question asks how several stimulus features jointly predict several response channels. Each kernel then represents one predictor-output pair while controlling for the other predictors.

In [1]:

import matplotlib
matplotlib.use("module://matplotlib_inline.backend_inline")
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import fftconvolve

from fftrf import TRF, pearsonr, r2_score

import matplotlib matplotlib.use("module://matplotlib_inline.backend_inline") import matplotlib.pyplot as plt import numpy as np from scipy.signal import fftconvolve from fftrf import TRF, pearsonr, r2_score

Simulate two stimulus features and two response channels¶

Time remains axis 0. One trial has shape (samples, 2) on both the predictor and response sides.

In [2]:

rng = np.random.default_rng(12)
fs = 128.0
tmin, tmax = 0.0, 0.25
times = np.arange(round(tmax * fs)) / fs
true_weights = np.zeros((2, times.size, 2))
true_weights[0, :, 0] = 0.8 * np.exp(-0.5 * ((times - 0.05) / 0.015) ** 2)
true_weights[1, :, 0] = -0.5 * np.exp(-0.5 * ((times - 0.12) / 0.025) ** 2)
true_weights[0, :, 1] = 0.45 * np.exp(-0.5 * ((times - 0.09) / 0.02) ** 2)
true_weights[1, :, 1] = 0.65 * np.exp(-0.5 * ((times - 0.17) / 0.03) ** 2)

stimulus, response = [], []
for _ in range(6):
    x = rng.standard_normal((2048, 2))
    y = np.column_stack([
        sum(fftconvolve(x[:, i], true_weights[i, :, o], mode="full")[:2048]
            for i in range(2))
        for o in range(2)
    ])
    y += 0.15 * rng.standard_normal(y.shape)
    stimulus.append(x)
    response.append(y)

train_x, test_x = stimulus[:-1], stimulus[-1]
train_y, test_y = response[:-1], response[-1]

rng = np.random.default_rng(12) fs = 128.0 tmin, tmax = 0.0, 0.25 times = np.arange(round(tmax * fs)) / fs true_weights = np.zeros((2, times.size, 2)) true_weights[0, :, 0] = 0.8 * np.exp(-0.5 * ((times - 0.05) / 0.015) ** 2) true_weights[1, :, 0] = -0.5 * np.exp(-0.5 * ((times - 0.12) / 0.025) ** 2) true_weights[0, :, 1] = 0.45 * np.exp(-0.5 * ((times - 0.09) / 0.02) ** 2) true_weights[1, :, 1] = 0.65 * np.exp(-0.5 * ((times - 0.17) / 0.03) ** 2) stimulus, response = [], [] for _ in range(6): x = rng.standard_normal((2048, 2)) y = np.column_stack([ sum(fftconvolve(x[:, i], true_weights[i, :, o], mode="full")[:2048] for i in range(2)) for o in range(2) ]) y += 0.15 * rng.standard_normal(y.shape) stimulus.append(x) response.append(y) train_x, test_x = stimulus[:-1], stimulus[-1] train_y, test_y = response[:-1], response[-1]

Fit one joint model¶

A scalar ridge value applies to every predictor. Use feature-banded regularization only when separate penalties are justified and selected within training data.

In [3]:

model = TRF(direction=1)
model.train(
    train_x,
    train_y,
    fs=fs,
    tmin=tmin,
    tmax=tmax,
    regularization=1e-2,
    segment_length=512,
    overlap=0.5,
    window="hann",
)
prediction, channel_r = model.predict(test_x, test_y, average=False)

print("weights:", model.weights.shape, "(inputs, lags, outputs)")
print("held-out r by output:", np.asarray(channel_r))

model = TRF(direction=1) model.train( train_x, train_y, fs=fs, tmin=tmin, tmax=tmax, regularization=1e-2, segment_length=512, overlap=0.5, window="hann", ) prediction, channel_r = model.predict(test_x, test_y, average=False) print("weights:", model.weights.shape, "(inputs, lags, outputs)") print("held-out r by output:", np.asarray(channel_r))

weights: (2, 32, 2) (inputs, lags, outputs)
held-out r by output: [0.9965767  0.99653847]

Inspect the full kernel bank¶

Do not interpret a single predictor kernel as a marginal effect. It is the linear contribution conditional on the other included predictors.

In [4]:

fig, axes = model.plot_grid(
    input_labels=["Envelope", "Onset-like feature"],
    output_labels=["EEG channel 1", "EEG channel 2"],
    sharey=False,
    title="Jointly estimated predictor-output kernels",
)
for index, ax in enumerate(np.asarray(axes).ravel()):
    input_index, output_index = divmod(index, 2)
    ax.plot(
        times * 1e3,
        true_weights[input_index, :, output_index],
        "--",
        color="black",
        linewidth=1,
    )
plt.show()

fig, axes = model.plot_grid( input_labels=["Envelope", "Onset-like feature"], output_labels=["EEG channel 1", "EEG channel 2"], sharey=False, title="Jointly estimated predictor-output kernels", ) for index, ax in enumerate(np.asarray(axes).ravel()): input_index, output_index = divmod(index, 2) ax.plot( times * 1e3, true_weights[input_index, :, output_index], "--", color="black", linewidth=1, ) plt.show()

Continue with Regularization if differently scaled predictor groups need separately selected penalties.