Modeling Cortical Responses to Naturalistic Speech

Author

Dr. Ole Bialas

Download Materials

import numpy as np
from matplotlib import pyplot as plt
import mne
import mtrf

If we want to study how the brain processes speech, we are faced with a fundamental limitation of traditional ERP-based approaches of analyzing EEG data: we cannot obtain average responses to repeated stimuli when presenting speech in a natural way (e.g. an audiobook). This is where linear models are useful - they allow us to estimate the relationship between the speech signal and neural response without having to epoch and average them.

The figure below illustrates the framework for modeling neural responses to naturalistic speech:

In the experiment, the participant listens to speech while their EEG is being recorded
From the speech signal, different features (envelope, spectrogram, phonemes) are extracted
These features are used as regressors to fit a linear model to the continuous neural recording
The estimated temporal response function (TRF) indicates the expected change in the neural response following a change in the predictor
By predicting the data, we can assess the model’s accuracy

From Crosse et al. (2021)

The model can be applied in two directions: as a forward (or encoding) model, predicting brain responses from speech features, or as a backward (or decoding) model, reconstructing stimulus features from the neural response. Forward models can be used to identify which speech features are encoded in the neural response as shown for example in Di Liberto et al. (2015). Backward models can be used to measure attention in a selective listening task by testing which auditory stream can be decoded more accurately from the neural response as shown in O’Sullivan et al. (2014).

In this notebook, you are going to use the mTRFpy package to estimate, visualize and evaluate TRF models in the context of naturalistic speech.

Setup

Data Preparation

We are going to use the sample data provided by mTRFpy. This data contains a 128-channel EEG recording of a participant listening to an audiobook version of Hemingway’s “The Old Man and the Sea” and the spectrogram of the corresponding audiobook segment at a sampling rate of 128 Hz. The data is split into 10 segments to allow for cross-validation when fitting the models.

spectrogram, response, fs = mtrf.load_sample_data(n_segments=10)
# crop all segments to the same length
min_len = min(r.shape[0] for r in response)
response = [r[:min_len] for r in response]
spectrogram = [s[:min_len] for s in spectrogram]

fs

np.uint8(128)

Let’s plot a segment of the EEG recording. It’s already preprocessed and ready for analysis.

t = np.linspace(0, 1000/fs, 1000)
fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(t, response[0][:1000], linewidth=0.4);
ax.set(xlabel="Time [s]", ylabel="Amplitude [μV]", xlim=(t.min(), t.max()));

The spectrogram contains the power of the speech signal divided into different bands across the frequency range that is relevant for human speech (80-8000 Hz).

fig, ax = plt.subplots(figsize=(9, 3))
im = ax.imshow(spectrogram[0][:1000].T, aspect="auto", origin="lower")
cbar = fig.colorbar(im, ax=ax)
cbar.set_label("Power [a.u.]")
xticks = np.linspace(0, 1000, 10).astype(int)
ax.set(xticks=xticks, xticklabels=np.round(xticks / fs, 1), xlabel="Time [s]", ylabel="Band");

We can also average the spectrogram across bands to obtain the acoustic envelope which describes the total energy in the speech signal across time.

envelope = [s.mean(axis=1) for s in spectrogram]
fig, ax = plt.subplots(figsize=(8, 3))
ax.plot(t, envelope[0][:1000]);
ax.set(xlabel="Time [s]", ylabel="Amplitude [a.u.]", xlim=(t.min(), t.max()));

Section 1: Encoding and Decoding Models with mTRFpy

Background

The mTRFpy package is a Python implementation of the Matlab mTRF toolbox and its purpose is to make it easy to fit, evaluate and visualize regularized linear models. The central element of the toolbox is the TRF() class which can be instantiated as a forward (default) or backward model (by setting direction=-1). This TRF can then be trained on the target stimulus and response across a range of time lags (in seconds). We can visualize the weights of the trained TRF to see the estimated temporal relationship between stimulus and response and use it to predict the response (forward modeling) or reconstruct the stimulus (backward modeling).

Exercises

In this section you are going to apply forward and backward TRF models using the mTRFpy sample data. You are going to train the models, visualize their weights and predict the original responses or stimuli while playing with the different TRF parameters to get a better understanding of the interface. Here are the essential code snippets:

Description	Code
Create a forward TRF model	`trf = mtrf.TRF()`
Create a backward TRF model	`trf = mtrf.TRF(direction=-1)`
Train the TRF	`trf.train(stimulus, response, fs, tmin=0, tmax=0.4, regularization=100)`
Predict the response and compute the score	`pred, score = trf.predict(stimulus, response)`
Plot the TRF	`trf.plot()`
Plot the TRF for a specific channel	`trf.plot(channel=105)`
Plot the global field power (GFP) as an image	`trf.plot(channel="gfp", kind="image")`
Convert a backward to a forward model	`trf_fwd = trf.to_forward(response)`

Example: Train a TRF on the response using the acoustic envelope with time lags between tmin=0 and tmax=0.3 with $\lambda=1$ (regularization).

trf = mtrf.TRF()
trf.train(envelope, response, fs, tmin=0, tmax=0.3, regularization=1)
trf.plot();

Exercise: Fit the same model again with $\lambda=100$ and plot the results.

Solution

trf = mtrf.TRF()
trf.train(envelope, response, fs, tmin=0, tmax=0.3, regularization=100)
trf.plot();

Exercise: Fit the same model again but increase the number of time lags to span the range from tmin=-0.1 to tmax=0.5 and plot the result.

Solution

trf = mtrf.TRF()
trf.train(envelope, response, fs, tmin=-0.1, tmax=0.5, regularization=100)
trf.plot();

Example: Use the trained TRF to predict the response from the envelope and plot a segment of the prediction and original response for channel 1.

Note: The second output of trf.predict is the prediction score which we assign to _ because we do not need it here.

response_pred, _ = trf.predict(envelope, response)
ch_idx = 1
plt.plot(response[0][:500, ch_idx], label="Original")
plt.plot(response_pred[0][:500, ch_idx], label="Prediction")
plt.legend()

Exercise: The cell below fits the same model but with $\lambda=0.1$ . Predict using the envelope and plot a segment of the predicted and original response. How did reducing $\lambda$ affect the model’s prediction?

trf = mtrf.TRF()
trf.train(envelope, response, fs, tmin=-0.1, tmax=0.5, regularization=0.1)

Solution

response_pred, _ = trf.predict(envelope, response)
ch_idx = 1
plt.plot(response[0][:500, ch_idx], label="Original")
plt.plot(response_pred[0][:500, ch_idx], label="Prediction")
plt.legend();

Exercise: The code below fits a backward model with direction=-1 that reconstructs the envelope from the neural response. Use trf.predict to reconstruct the envelope and plot a segment of the original envelope together with the reconstruction.

Note: Fitting this model takes longer because there are 128 predictors (i.e. EEG channels) now.

trf = mtrf.TRF(direction=-1)
trf.train(envelope, response, fs, tmin=0, tmax=0.4, regularization=10)

Solution

envelope_pred, _ = trf.predict(envelope, response)
plt.plot(envelope[0][:500], label="Original")
plt.plot(envelope_pred[0][:500], label="Prediction")
plt.legend();

Exercise: Plot the TRF from the previous exercise. What does the warning say? Why are the time lags negative?

Solution

trf.plot();

WARNING: decoder weights are hard to interpret, consider using the `to_forward()` method

Exercise: The cell below transforms the backward to a forward model with interpretable weights (for details see Haufe et al. 2014). Plot the weights of trf_fwd.

trf_fwd = trf.to_forward(response)

Solution

trf_fwd.plot();

Exercise: Fit a forward TRF to predict the response from the spectrogram with tmin=0, tmax=0.4 and regularization=1000. Plot it with channel=105 to visualize the TRF at this (left auditory) channel for every spectral band.

Solution

trf = mtrf.TRF()
trf.train(spectrogram, response, fs, tmin=0, tmax=0.4, regularization=1000)
trf.plot(channel=105);

Exercise: Plot the TRF again with channel="gfp" to plot the global field power (i.e. the standard deviation across channels).

Solution

trf.plot(channel="gfp");

Exercise: Plot the TRF again with channel="gfp" and kind="image" to plot the global field power as an image which allows you to visually separate the frequency bands.

Solution

trf.plot(channel="gfp", kind="image");

Section 2: Optimizing Models and Controlling for Overfitting

Background

Because TRFs, especially with many features, tend to overfit on the training data, they give us an unrealistic picture of the model’s accuracy if we test it on the same data. To get an unbiased estimate of model accuracy, we can use cross-validation:

train the model on all but one segment of the data
test it on the withheld segment
rotate the test and training data and average accuracy across all segments

The same method can also be used to find the best value of the regularization parameter $\lambda$ from a list of candidate values. When optimizing the value of $\lambda$ it also matters what we use as a metric for evaluating the model — using Pearson’s correlation and using the explained variance $R^2$ may result in a different best $\lambda$ .

Exercises

In this section you are going to use cross-validation to get an unbiased estimate of model accuracy and optimize $\lambda$ for different models. To evaluate a given model, you’ll use mtrf.stats.crossval. To find the best value for $\lambda$ you can simply pass an array of candidate values to the regularization argument of trf.train(). Here are the code examples:

Description	Code
Predict and get the in-sample score	`pred, score = trf.predict(stimulus, response)`
Get the cross-validated score	`score = mtrf.stats.crossval(trf, stimulus, response, fs, tmin=0, tmax=0.4, regularization=10)`
Create a log-spaced array of candidate $\lambda$ values	`lambdas = np.logspace(-1, 3, 20)`
Train with cross-validation to find the best $\lambda$	`scores = trf.train(stimulus, response, fs, tmin=0, tmax=0.4, regularization=lambdas)`
Retrieve the best $\lambda$ after cross-validation training	`trf.regularization`
Use a custom evaluation metric	`trf = mtrf.TRF(metric=my_metric_fn)`

Example: Fit a backward TRF with $\lambda=10$ to reconstruct the envelope from the response and estimate the score (Pearson’s correlation).

To speed up computation, we are only selecting every 4th EEG channel here.

response_reduced = [r[:,::4] for r in response]

trf = mtrf.TRF(direction=-1)
trf.train(envelope, response_reduced, fs, tmin=0, tmax=0.3, regularization=10)
response_pred, score = trf.predict(envelope, response_reduced)
score

np.float64(0.3370924590778743)

Example: Use cross-validation to get an unbiased estimate of the model’s accuracy. What does the difference between this and the previous score tell you about the model?

mtrf.stats.crossval(trf, envelope, response_reduced, fs, tmin=0, tmax=0.3, regularization=10)

Cross-validating[..................................................] 0/10
Cross-validating[#####.............................................] 1/10
Cross-validating[##########........................................] 2/10
Cross-validating[###############...................................] 3/10
Cross-validating[####################..............................] 4/10
Cross-validating[#########################.........................] 5/10
Cross-validating[##############################....................] 6/10
Cross-validating[###################################...............] 7/10
Cross-validating[########################################..........] 8/10
Cross-validating[#############################################.....] 9/10
Cross-validating[##################################################] 10/10

np.float64(0.16396059598506602)

Exercise: Fit a forward TRF with $\lambda=1$ to predict the response from the envelope and estimate the score.

Solution

trf = mtrf.TRF()
trf.train(envelope, response, fs, tmin=0, tmax=0.3, regularization=1)
response_pred, score = trf.predict(envelope, response)
score

np.float64(0.08848668110279839)

Exercise: Use cross-validation to get an unbiased estimate of the model’s accuracy. Is this model more or less prone to overfitting compared to the backward model?

Solution

mtrf.stats.crossval(trf, envelope, response, fs, tmin=0, tmax=0.3, regularization=1)

Cross-validating[..................................................] 0/10
Cross-validating[#####.............................................] 1/10
Cross-validating[##########........................................] 2/10
Cross-validating[###############...................................] 3/10
Cross-validating[####################..............................] 4/10
Cross-validating[#########################.........................] 5/10
Cross-validating[##############################....................] 6/10
Cross-validating[###################################...............] 7/10
Cross-validating[########################################..........] 8/10
Cross-validating[#############################################.....] 9/10
Cross-validating[##################################################] 10/10

np.float64(0.05447728655513813)

Exercise: Fit a forward TRF with $\lambda=1000$ to predict the response from the spectrogram and estimate the score with and without cross-validation.

Solution

trf = mtrf.TRF()
trf.train(spectrogram, response, fs, tmin=0, tmax=0.3, regularization=1000)
response_pred, score = trf.predict(spectrogram, response)
score_crossval = mtrf.stats.crossval(trf, spectrogram, response, fs, tmin=0, tmax=0.3, regularization=1000)
score, score_crossval

Cross-validating[..................................................] 0/10
Cross-validating[#####.............................................] 1/10
Cross-validating[##########........................................] 2/10
Cross-validating[###############...................................] 3/10
Cross-validating[####################..............................] 4/10
Cross-validating[#########################.........................] 5/10
Cross-validating[##############################....................] 6/10
Cross-validating[###################################...............] 7/10
Cross-validating[########################################..........] 8/10
Cross-validating[#############################################.....] 9/10
Cross-validating[##################################################] 10/10

(np.float64(0.09652326916733048), np.float64(0.03045756899181804))

Example: Create a log-spaced array of lambdas between $10^{-1}$ and $10^2$ to train a forward TRF to predict the response from the envelope with cross-validation. Plot the scores and highlight the best $\lambda$ (trf.regularization).

Note: You can set verbose=False to suppress the progressbar.

lambdas = np.logspace(-1, 2, 10)
trf = mtrf.TRF()
scores = trf.train(
    envelope, response, fs, tmin=0, tmax=0.3, regularization=lambdas, verbose=False)

plt.semilogx(lambdas, scores)
plt.axvline(x=trf.regularization, color="red")

Exercise: Train a forward TRF to predict the response from the spectrogram with cross-validation and adjust the range of lambda as needed. Plot the scores and highlight the best $\lambda$ .

Solution

lambdas = np.logspace(-1, 5, 10)
trf = mtrf.TRF()
scores = trf.train(
    spectrogram, response, fs, tmin=0, tmax=0.3, regularization=lambdas, verbose=False)

plt.semilogx(lambdas, scores)
plt.axvline(x=trf.regularization, color="red")

Exercise: The cell below defines a function for calculating $R^2$ . Create a forward TRF with metric=rsquared to predict the response from the envelope and fit it with cross-validation. How does the curve of scores differ from the previous example using Pearson’s correlation?

def rsquared(y, y_pred):
    ss_res = np.sum((y - y_pred) ** 2)
    ss_tot = np.sum((y - y.mean()) ** 2)
    return 1 - ss_res / ss_tot

Solution

lambdas = np.logspace(-1, 3, 10)
trf = mtrf.TRF(metric=rsquared)
scores = trf.train(
    envelope, response, fs, tmin=0, tmax=0.3, regularization=lambdas, verbose=False)

plt.semilogx(lambdas, scores)
plt.axvline(x=trf.regularization, color="red")

Section 3: Interfacing with MNE

Background

The mtrf.TRF() class provides a method to convert the weights of a trained model to a list with one mne.Evoked object for every feature in the model (e.g. 16 for the spectrogram). This allows us to use MNE’s visualization methods to plot the TRF weights as ERP-like waveforms and topographical maps. We can also obtain a per-channel prediction score and plot it as a topographical map to see which channels are best predicted by the model.

Exercises

In the following exercises you are going to convert trained TRFs to mne.Evoked objects and visualize them. Here are the relevant code snippets:

Description	Code
Load a standard channel montage	`montage = mne.channels.make_standard_montage("biosemi128")`
Convert TRF weights to a list of Evoked objects	`evokeds = trf.to_mne_evoked(montage)`
Plot an Evoked object with topomaps at specific times	`evokeds[0].plot_joint(times=[0.1, 0.2])`
Get per-channel prediction scores	`pred, score = trf.predict(stimulus, response, average=False)`
Plot scores as a topographic map	`mne.viz.plot_topomap(score, pos=evokeds[0].info, size=3)`

Run the cell below to fit a forward TRF that predicts the response from the spectrogram.

trf = mtrf.TRF()
trf.train(spectrogram, response, fs, tmin=0, tmax=0.3, regularization=1000)

Example: Use trf.to_mne_evoked to return a list with one mne.Evoked object for every feature in the TRF and plot the first element in evokeds.

Note: The montage provides information on the channel locations.

montage = mne.channels.make_standard_montage("biosemi128")
evokeds = trf.to_mne_evoked(montage)
evokeds[0].plot_joint(times=[0.175, 0.26]);

No projector specified for this dataset. Please consider the method self.add_proj.

Exercise: Plot the 4th element in evokeds.

Solution

evokeds[3].plot_joint(times=[0.1, 0.14, 0.17]);

No projector specified for this dataset. Please consider the method self.add_proj.

Exercise: The cell below fits a forward TRF that predicts the response from the envelope. Convert it .to_mne_evoked() and plot it.

trf = mtrf.TRF()
trf.train(envelope, response, fs, tmin=0, tmax=0.3, regularization=1)

Solution

evokeds = trf.to_mne_evoked(montage)
evokeds[0].plot_joint(times=[0.115, 0.17, 0.22, 0.265]);

No projector specified for this dataset. Please consider the method self.add_proj.

Example: Estimate the accuracy of the trained TRF with average=False to get one score for every channel and plot it as a topographical map using the channel positions from evokeds[0].info. Where are you seeing the most accurate predictions?

response_pred, score = trf.predict(envelope, response, average=False)
mne.viz.plot_topomap(score, pos=evokeds[0].info, size=3);

Exercise: The cell below fits a forward TRF that predicts the response from the spectrogram and estimates the score for every channel. Plot the scores as a topographical map.

trf = mtrf.TRF()
trf.train(spectrogram, response, fs, tmin=0, tmax=0.3, regularization=1000)
response_pred, score = trf.predict(spectrogram, response, average=False)

Solution

mne.viz.plot_topomap(score, pos=evokeds[0].info, size=3)

Section 4: Statistical Inference

Background

You saw that the absolute correlation between the neural response and the models’ predictions is very low in absolute terms. This is to be expected since EEG data is noisy and only a small portion of the data can be predicted from the stimulus. However, this raises a question: is the model reliably predicting brain activity or are the results just due to chance? To answer this question, we need to perform a statistical test.

To estimate the significance of our results we can permute the data. We can randomly scramble the segments and apply our model to get the prediction accuracy for randomized data. By repeating this permutation many times we get a distribution of scores that is expected when there is no real relationship between stimulus and response. Comparing the observed score to the permutation distribution gives us the p-value for our result.

Exercises

In this section, you are going to use mtrf.stats.permutation_distribution() to test the significance of the model scores. Here are the code examples:

Description	Code
Set a random seed for reproducibility	`np.random.seed(42)`
Generate a permutation distribution	`scores_permute = mtrf.stats.permutation_distribution(trf, stimulus, response, fs, tmin=0, tmax=0.4, regularization=10, n_permute=100)`
Compute the p-value	`p = (scores_permute > score).sum() / len(scores_permute)`

Example: The code below fits a forward TRF that predicts the response from the envelope with cross-validation to find the best $\lambda$ and estimates the model’s score. Create a permutation_distribution() and compare the score to the distribution to get the p-value, then plot the result.

lambdas = np.logspace(-1, 2, 10)

trf = mtrf.TRF()
scores = trf.train(
    envelope, response, fs, tmin=0, tmax=0.3, regularization=lambdas, verbose=False)
score = scores.max()

score

np.float64(0.05447728655513813)

Solution

np.random.seed(42)
scores_permute = mtrf.stats.permutation_distribution(
    trf, envelope, response, fs, tmin=0, tmax=0.3, regularization=trf.regularization, n_permute=100, verbose=False
    )

plt.hist(scores_permute, label="permutation scores");
plt.axvline(score, color="red", label="observed")
plt.xlabel("scores")
plt.ylabel("N observations")
plt.title(f"p = {(scores_permute>score).sum()/len(scores_permute)}");

Exercise: The code below fits a forward TRF that predicts the response from the spectrogram with cross-validation to find the best $\lambda$ and estimates the model’s score. Create a permutation_distribution() and compare the score to the distribution to get the p-value, then plot the result. Why does the result differ from the previous exercise?

lambdas = np.logspace(-1, 5, 10)

trf = mtrf.TRF()
scores = trf.train(
    spectrogram, response, fs, tmin=0, tmax=0.3, regularization=lambdas, verbose=False)
score = scores.max()

score

np.float64(0.03269672164172045)

Solution

np.random.seed(42)
scores_permute = mtrf.stats.permutation_distribution(
    trf, spectrogram, response, fs, tmin=0, tmax=0.3, regularization=trf.regularization, n_permute=100, verbose=False
    )

plt.hist(scores_permute, label="permutation scores");
plt.axvline(score, color="red", label="observed")
plt.xlabel("scores")
plt.ylabel("N observations")
plt.title(f"p = {(scores_permute>score).sum()/len(scores_permute)}");

Exercise: The code below fits a backward TRF that reconstructs the envelope from the response with cross-validation to find the best $\lambda$ and estimates the model’s score. Create a permutation_distribution() and compare the score to the distribution to get the p-value, then plot the result.

response_reduced = [r[:,::4] for r in response]  # reduce size to speed up computation
lambdas = np.logspace(-1, 5, 10)

trf = mtrf.TRF(direction=-1)
scores = trf.train(
    envelope, response_reduced, fs, tmin=0, tmax=0.3, regularization=lambdas, verbose=False)
score = scores.max()
score

np.float64(0.1639605959850661)

Solution

np.random.seed(42)
scores_permute = mtrf.stats.permutation_distribution(
    trf, envelope, response_reduced, fs, tmin=0, tmax=0.3, regularization=trf.regularization, n_permute=100, verbose=False
    )

plt.hist(scores_permute, label="permutation scores");
plt.axvline(score, color="red", label="observed")
plt.xlabel("scores")
plt.ylabel("N observations")
plt.title(f"p = {(scores_permute>score).sum()/len(scores_permute)}");