Backward Decoding¶

A backward model reconstructs a stimulus feature from multichannel responses. This is useful for questions about decodability or stimulus tracking, but decoder weights are multivariate filters and should not be interpreted like forward neural response kernels.

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 one envelope and four response channels¶

In [2]:

rng = np.random.default_rng(21)
fs = 128.0
tmin, tmax = 0.0, 0.25
times = np.arange(round(tmax * fs)) / fs
channel_kernels = np.stack([
    amplitude * np.exp(-0.5 * ((times - delay) / 0.025) ** 2)
    for amplitude, delay in [(0.9, 0.05), (0.6, 0.08), (-0.5, 0.11), (0.4, 0.15)]
])

stimulus, response = [], []
for _ in range(6):
    envelope = rng.standard_normal((2048, 1))
    eeg = np.column_stack([
        fftconvolve(envelope[:, 0], kernel, mode="full")[:2048]
        for kernel in channel_kernels
    ])
    eeg += 0.35 * rng.standard_normal(eeg.shape)
    stimulus.append(envelope)
    response.append(eeg)

train_stimulus, test_stimulus = stimulus[:-1], stimulus[-1]
train_response, test_response = response[:-1], response[-1]

rng = np.random.default_rng(21) fs = 128.0 tmin, tmax = 0.0, 0.25 times = np.arange(round(tmax * fs)) / fs channel_kernels = np.stack([ amplitude * np.exp(-0.5 * ((times - delay) / 0.025) ** 2) for amplitude, delay in [(0.9, 0.05), (0.6, 0.08), (-0.5, 0.11), (0.4, 0.15)] ]) stimulus, response = [], [] for _ in range(6): envelope = rng.standard_normal((2048, 1)) eeg = np.column_stack([ fftconvolve(envelope[:, 0], kernel, mode="full")[:2048] for kernel in channel_kernels ]) eeg += 0.35 * rng.standard_normal(eeg.shape) stimulus.append(envelope) response.append(eeg) train_stimulus, test_stimulus = stimulus[:-1], stimulus[-1] train_response, test_response = response[:-1], response[-1]

Fit in the response-to-stimulus direction¶

The user-facing interval is specified as [0, tmax). To follow mTRF decoder semantics, model.times stores the reversed physical lags ending at zero.

In [3]:

decoder = TRF(direction=-1)
decoder.train(
    stimulus=train_stimulus,
    response=train_response,
    fs=fs,
    tmin=tmin,
    tmax=tmax,
    regularization=np.logspace(-3, 1, 5),
    k="loo",
    segment_length=512,
    overlap=0.5,
    window="hann",
)
decoded, heldout_r = decoder.predict(
    response=test_response,
    stimulus=test_stimulus,
)

print(f"held-out decoding r: {float(heldout_r):.3f}")
print("stored physical lag range:", decoder.times[[0, -1]])
print("decoder weights:", decoder.weights.shape)

decoder = TRF(direction=-1) decoder.train( stimulus=train_stimulus, response=train_response, fs=fs, tmin=tmin, tmax=tmax, regularization=np.logspace(-3, 1, 5), k="loo", segment_length=512, overlap=0.5, window="hann", ) decoded, heldout_r = decoder.predict( response=test_response, stimulus=test_stimulus, ) print(f"held-out decoding r: {float(heldout_r):.3f}") print("stored physical lag range:", decoder.times[[0, -1]]) print("decoder weights:", decoder.weights.shape)

held-out decoding r: 0.471
stored physical lag range: [-0.2421875  0.       ]
decoder weights: (4, 32, 1)

Plot prediction first, weights second¶

Held-out reconstruction is the primary decoder result. Correlated EEG predictors can yield different weight patterns with similar predictions.

In [4]:

fig, axes = plt.subplots(2, 1, figsize=(9, 6), constrained_layout=True)
display = slice(0, 400)
sample_time = np.arange(test_stimulus.shape[0]) / fs
axes[0].plot(sample_time[display], test_stimulus[display, 0], label="Observed")
axes[0].plot(sample_time[display], decoded[display, 0], label="Decoded")
axes[0].set(xlabel="Time (s)", ylabel="Envelope", title="Untouched test trial")
axes[0].legend()

for channel in range(decoder.weights.shape[0]):
    axes[1].plot(decoder.times * 1e3, decoder.weights[channel, :, 0],
                 label=f"Channel {channel + 1}")
axes[1].set(xlabel="Physical lag (ms)", ylabel="Weight",
            title="Multichannel decoder weights")
axes[1].legend(ncols=2)
plt.show()

fig, axes = plt.subplots(2, 1, figsize=(9, 6), constrained_layout=True) display = slice(0, 400) sample_time = np.arange(test_stimulus.shape[0]) / fs axes[0].plot(sample_time[display], test_stimulus[display, 0], label="Observed") axes[0].plot(sample_time[display], decoded[display, 0], label="Decoded") axes[0].set(xlabel="Time (s)", ylabel="Envelope", title="Untouched test trial") axes[0].legend() for channel in range(decoder.weights.shape[0]): axes[1].plot(decoder.times * 1e3, decoder.weights[channel, :, 0], label=f"Channel {channel + 1}") axes[1].set(xlabel="Physical lag (ms)", ylabel="Weight", title="Multichannel decoder weights") axes[1].legend(ncols=2) plt.show()