Kernel Uncertainty and Held-Out Significance¶

Bootstrap intervals and permutation tests answer different questions. Trial bootstrap intervals describe variability of the fitted kernel across sampled trials. A held-out permutation test asks whether prediction exceeds a prespecified surrogate null.

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

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(
    seed=66,
    n_trials=8,
    noise_scale=0.25,
)
train_x, test_x = stimulus[:-1], stimulus[-1]
train_y, test_y = response[:-1], response[-1]

fs, tmin, tmax, times, true_kernel, stimulus, response = simulate_trials( seed=66, n_trials=8, noise_scale=0.25, ) train_x, test_x = stimulus[:-1], stimulus[-1] train_y, test_y = response[:-1], response[-1]

Estimate a pointwise trial-bootstrap interval¶

The resampling unit is the trial. Forty resamples keep this rendered example quick; use substantially more for a reported analysis.

In [4]:

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",
    bootstrap_samples=40,
    bootstrap_level=0.95,
    bootstrap_seed=7,
)

interval, interval_times = model.bootstrap_interval_at()
print("interval:", interval.shape)

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", bootstrap_samples=40, bootstrap_level=0.95, bootstrap_seed=7, ) interval, interval_times = model.bootstrap_interval_at() print("interval:", interval.shape)

interval: (2, 1, 32, 1)

In [5]:

fig, ax = model.plot(
    show_bootstrap_interval=True,
    title="Pointwise trial-bootstrap interval",
)
ax.plot(times * 1e3, true_kernel, "--", color="black", label="True")
ax.legend()
plt.show()

fig, ax = model.plot( show_bootstrap_interval=True, title="Pointwise trial-bootstrap interval", ) ax.plot(times * 1e3, true_kernel, "--", color="black", label="True") ax.legend() plt.show()

The shaded region is pointwise across lags, not a simultaneous confidence band. It also does not provide population inference if the resampled trials all belong to one participant.

Test tracking on the untouched trial¶

Circular shifts preserve much of each signal's autocorrelation while breaking its temporal alignment. The model is held fixed; the test does not repeat model selection.

In [6]:

result = model.permutation_test(
    stimulus=test_x,
    response=test_y,
    n_permutations=99,
    surrogate="circular_shift",
    seed=7,
)
print(f"observed r: {float(result.observed_score):.3f}")
print(f"p-value: {float(result.p_value):.3f}")
print("minimum attainable p:", 1 / (result.n_permutations + 1))

result = model.permutation_test( stimulus=test_x, response=test_y, n_permutations=99, surrogate="circular_shift", seed=7, ) print(f"observed r: {float(result.observed_score):.3f}") print(f"p-value: {float(result.p_value):.3f}") print("minimum attainable p:", 1 / (result.n_permutations + 1))

observed r: 0.992
p-value: 0.010
minimum attainable p: 0.01

In [7]:

fig, ax = plt.subplots(figsize=(7, 4))
ax.hist(result.null_scores, bins=18, color="#9ecae1", edgecolor="white")
ax.axvline(result.observed_score, color="#c44e52", linewidth=2,
           label="Observed")
ax.set(xlabel="Held-out Pearson r", ylabel="Surrogates",
       title="Circular-shift null distribution")
ax.legend()
plt.show()

fig, ax = plt.subplots(figsize=(7, 4)) ax.hist(result.null_scores, bins=18, color="#9ecae1", edgecolor="white") ax.axvline(result.observed_score, color="#c44e52", linewidth=2, label="Observed") ax.set(xlabel="Held-out Pearson r", ylabel="Surrogates", title="Circular-shift null distribution") ax.legend() plt.show()

Prespecify the surrogate, shift restrictions, tail, and evaluation unit. See Significance Testing and Trial Weighting and Bootstrap.