Getting Started: Fit, Select, and Test a Forward TRF¶

This notebook follows the minimum defensible workflow for a neuroscience analysis: simulate several trials, use training trials to select regularization, and evaluate once on an untouched trial. The same array conventions apply to speech features and EEG/MEG.

Imports¶

A standard installation already includes NumPy, SciPy, and Matplotlib. TRF(direction=1) means stimulus-to-response encoding.

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

Create data with a known temporal response¶

The simulation lets us check both prediction and kernel recovery. In an empirical analysis the true kernel is unknown, so only held-out prediction and reproducibility across observations remain available as validation.

In [2]:

def simulate_trials(*, seed=7, n_trials=6, n_samples=2048, noise_scale=0.12):
    rng = np.random.default_rng(seed)
    fs = 128.0
    tmin, tmax = 0.0, 0.25
    times = np.arange(round(tmin * fs), round(tmax * fs)) / fs
    true_kernel = (
        0.9 * np.exp(-0.5 * ((times - 0.055) / 0.014) ** 2)
        - 0.5 * np.exp(-0.5 * ((times - 0.125) / 0.024) ** 2)
    )

    stimulus, response = [], []
    for _ in range(n_trials):
        x = rng.standard_normal((n_samples, 1))
        y = fftconvolve(x[:, 0], true_kernel, mode="full")[:n_samples]
        y += noise_scale * rng.standard_normal(n_samples)
        stimulus.append(x)
        response.append(y[:, np.newaxis])
    return fs, tmin, tmax, times, true_kernel, stimulus, response

def simulate_trials(*, seed=7, n_trials=6, n_samples=2048, noise_scale=0.12): rng = np.random.default_rng(seed) fs = 128.0 tmin, tmax = 0.0, 0.25 times = np.arange(round(tmin * fs), round(tmax * fs)) / fs true_kernel = ( 0.9 * np.exp(-0.5 * ((times - 0.055) / 0.014) ** 2) - 0.5 * np.exp(-0.5 * ((times - 0.125) / 0.024) ** 2) ) stimulus, response = [], [] for _ in range(n_trials): x = rng.standard_normal((n_samples, 1)) y = fftconvolve(x[:, 0], true_kernel, mode="full")[:n_samples] y += noise_scale * rng.standard_normal(n_samples) stimulus.append(x) response.append(y[:, np.newaxis]) return fs, tmin, tmax, times, true_kernel, stimulus, response

In [3]:

fs, tmin, tmax, times, true_kernel, stimulus, response = simulate_trials()

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

print(f"{len(train_stimulus)} training trials")
print("one stimulus trial:", train_stimulus[0].shape)
print("one response trial:", train_response[0].shape)

fs, tmin, tmax, times, true_kernel, stimulus, response = simulate_trials() train_stimulus, test_stimulus = stimulus[:-1], stimulus[-1] train_response, test_response = response[:-1], response[-1] print(f"{len(train_stimulus)} training trials") print("one stimulus trial:", train_stimulus[0].shape) print("one response trial:", train_response[0].shape)

5 training trials
one stimulus trial: (2048, 1)
one response trial: (2048, 1)

Select ridge regularization using only training trials¶

Leave-one-trial-out CV treats the trial as the validation unit. The final held-out trial is not involved in choosing lambda. Whole-trial spectra are used here as the closest baseline to a conventional finite-lag TRF fit.

In [4]:

ridge_grid = np.logspace(-4, 1, 6)
model = TRF(direction=1, metric="pearsonr")
cv_scores = model.train(
    stimulus=train_stimulus,
    response=train_response,
    fs=fs,
    tmin=tmin,
    tmax=tmax,
    regularization=ridge_grid,
    k="loo",
    seed=7,
    segment_length=None,
    window=None,
)

print("selected lambda:", model.regularization)
print("CV score shape:", np.asarray(cv_scores).shape)

ridge_grid = np.logspace(-4, 1, 6) model = TRF(direction=1, metric="pearsonr") cv_scores = model.train( stimulus=train_stimulus, response=train_response, fs=fs, tmin=tmin, tmax=tmax, regularization=ridge_grid, k="loo", seed=7, segment_length=None, window=None, ) print("selected lambda:", model.regularization) print("CV score shape:", np.asarray(cv_scores).shape)

selected lambda: 1.0
CV score shape: (6,)

Evaluate once on the untouched trial¶

Pearson correlation measures tracking shape and is insensitive to scale. R² additionally penalizes scale and offset errors, so reporting both can reveal a poorly calibrated prediction.

In [5]:

prediction, heldout_r = model.predict(
    stimulus=test_stimulus,
    response=test_response,
)
heldout_r2 = float(r2_score(test_response, prediction).mean())
kernel_r = float(np.corrcoef(true_kernel, model.weights[0, :, 0])[0, 1])

print(f"held-out Pearson r: {float(heldout_r):.3f}")
print(f"held-out R²: {heldout_r2:.3f}")
print(f"kernel correlation to ground truth: {kernel_r:.3f}")

prediction, heldout_r = model.predict( stimulus=test_stimulus, response=test_response, ) heldout_r2 = float(r2_score(test_response, prediction).mean()) kernel_r = float(np.corrcoef(true_kernel, model.weights[0, :, 0])[0, 1]) print(f"held-out Pearson r: {float(heldout_r):.3f}") print(f"held-out R²: {heldout_r2:.3f}") print(f"kernel correlation to ground truth: {kernel_r:.3f}")

held-out Pearson r: 0.998
held-out R²: 0.996
kernel correlation to ground truth: 1.000

Inspect the recovered kernel and held-out prediction¶

The kernel plot answers when the predictor contributes to the response. The trace plot answers whether that fitted mapping generalizes to unseen data.

In [6]:

fig, axes = plt.subplots(2, 1, figsize=(9, 6), constrained_layout=True)
axes[0].plot(times * 1e3, true_kernel, "--", color="black", label="True")
axes[0].plot(model.times * 1e3, model.weights[0, :, 0], label="Recovered")
axes[0].set(xlabel="Lag (ms)", ylabel="Weight", title="Lag-domain kernel")
axes[0].legend()

display = slice(0, 400)
sample_time = np.arange(test_response.shape[0]) / fs
axes[1].plot(sample_time[display], test_response[display, 0], label="Observed")
axes[1].plot(sample_time[display], prediction[display, 0], label="Predicted")
axes[1].set(xlabel="Time (s)", ylabel="Response", title="Untouched test trial")
axes[1].legend()
plt.show()

fig, axes = plt.subplots(2, 1, figsize=(9, 6), constrained_layout=True) axes[0].plot(times * 1e3, true_kernel, "--", color="black", label="True") axes[0].plot(model.times * 1e3, model.weights[0, :, 0], label="Recovered") axes[0].set(xlabel="Lag (ms)", ylabel="Weight", title="Lag-domain kernel") axes[0].legend() display = slice(0, 400) sample_time = np.arange(test_response.shape[0]) / fs axes[1].plot(sample_time[display], test_response[display, 0], label="Observed") axes[1].plot(sample_time[display], prediction[display, 0], label="Predicted") axes[1].set(xlabel="Time (s)", ylabel="Response", title="Untouched test trial") axes[1].legend() plt.show()

Choose the next notebook from your question¶

• Multiple predictors or EEG outputs: Multiple Features and Channels
• Response-to-stimulus reconstruction: Backward Decoding
• Alternative ridge structures: Regularization
• Spectral interpretation: Frequency-Resolved Analysis or Spectral Diagnostics