Analyzing results¶
Below, we will walk through an example analysis of TRF models that predict EEG responses from acoustic features. The examples use a subset of data from a larger EEG study in which participants listened to an audiobook presented in quiet either monaurally and binaurally. We'll examine two common TRF outputs: prediction accuracy, quantified as Pearson’s correlation r between the predicted and observed EEG responses, and the fitted TRFs. Again, three models were estimated: a full model containing both the 8-band gammatone envelope and 8-band onset predictors, and two reduced models containing only one of these predictor types.
As you'll see, we will make frequent use of functions and utilities from the eelbrain package (e.g., plotting and mass-univariate permutation testing). If you're unsure about the usage of a specific function, consult the Eelbrain documentation. In addition, this tutorial is not meant to provide a comprehensive overview of TRF analyses or to define what is a right or standard analysis workflow. This question depends on things like the predictors you include, your experiment design, and your hypotheses. You're encouraged to check the relavant literature to decide what makes sense to do.
Loading in the results¶
Let’s begin by loading the TRF results, which are stored as an eelbrain.Dataset object. Check the comments in the code chunks below for what each line is doing:
# Load packages
import os
import eelbrain
import numpy as np
import pandas as pd
from pathlib import Path
import matplotlib.pyplot as plt
# Paths to the folder containing the results and to the metadata spreadsheet
results_path = 'examples/example_results'
metadata_path = 'examples/example_metadata.csv'
# The column in the metadata listing the EEG filenames
filename_col = 'filename'
# Columns from the meatdata to include in the results dataset
cols_to_include = ['filename', 'sub', 'session', 'intensity', 'modality']
# Load the metadata. We'll add the information here to the results dataset later.
metadata = pd.read_csv(metadata_path)
metadata.head()
| filename | sub | session | intensity | modality | stim | |
|---|---|---|---|---|---|---|
| 0 | sub-01_ses-2_task-60bi_acq-2.fif | 1 | 2 | 60 | bi | track16.wav+track17.wav+track18.wav+track19.wa... |
| 1 | sub-01_ses-2_task-60mo_acq-1.fif | 1 | 2 | 60 | mo | track1.wav+track2.wav+track3.wav+track4.wav+tr... |
| 2 | sub-02_ses-2_task-60bi_acq-3.fif | 2 | 2 | 60 | bi | track31.wav+track32.wav+track33.wav+track34.wa... |
| 3 | sub-02_ses-2_task-60mo_acq-4.fif | 2 | 2 | 60 | mo | track46.wav+track47.wav+track48.wav+track49.wa... |
| 4 | sub-03_ses-2_task-60bi_acq-1.fif | 3 | 2 | 60 | bi | track1.wav+track2.wav+track3.wav+track4.wav+tr... |
# Get the list of all result folders. Each folder should contain all the BoostingResults objects
results_path = Path(results_path)
result_folders = [str(x) for x in results_path.iterdir() if x.is_dir()]
print(f'Found {len(result_folders)} result folders.')
Found 40 result folders.
# Initialize an empty list to store the results datasets
ds_list = []
# A helper function to add a Case dimension with length of 1 to NDVar
def add_case_dim(ndvar):
return eelbrain.NDVar(ndvar.get_data()[np.newaxis, ...], (eelbrain.Case(1), *ndvar.dims))
for result_folder in result_folders:
# Get paths to:
# 1. The full model containng 8-band envelope and onsets
# 2. The model containing only the 8-band envelope
# 3. The model containing only the 8-band onsets
# The filename should have the format: 'BoostingResults-MODEL_NAME.pickle'
full_result_path = os.path.join(result_folder, 'BoostingResults-full.pickle')
envelope_only_result_path = os.path.join(result_folder, 'BoostingResults-envelope_only.pickle')
onsets_only_result_path = os.path.join(result_folder, 'BoostingResults-onsets_only.pickle')
# First, extract the prediction accuracy (correlation r) and fitted TRFs for the full model
full_results = eelbrain.load.unpickle(full_result_path)
full_r = full_results.r
# full_h contains the TRF weights for all predictors in the full model. Note that it's a tuple with two elements because our full model
# is specified as [gammatone-8, gammatone-on-8]. So to get the TRF weights for the 8-band envelope, we need to access full_h[0], and
# for the 8-band onsets, we need to access full_h[1]
full_h = full_results.h
full_h_envelope = full_h[0]
full_h_onsets = full_h[1]
# Also add TRF weights scaled to the predictors. Will explain the differences between h and h_scaled below
full_h_scaled_envelope = full_results.h_scaled[0]
full_h_scaled_onsets = full_results.h_scaled[1]
# Next, load the envelope-only model results
# For this, we just need the prediction accuracy
envelope_only_results = eelbrain.load.unpickle(envelope_only_result_path)
envelope_only_r = envelope_only_results.r
# Then load the onsets-only model results
# Again, we just need the prediction accuracy
onsets_only_results = eelbrain.load.unpickle(onsets_only_result_path)
onsets_only_r = onsets_only_results.r
# Finally, calculate the unique contribution to prediction accuracy for each predictor (Δr)
# Note that unique contribution of A is the r of full model minus the r of the reduced model "without A"
onsets_delta_r = full_r - envelope_only_r
envelope_delta_r = full_r - onsets_only_r
# Now build an Eelbrain dataset to store all results
# Initiate an empty dataset and add in the included metadata
ds = eelbrain.Dataset()
filename = os.path.basename(result_folder) + '.fif'
metadata_subset = metadata[metadata[filename_col] == filename]
for c in cols_to_include:
ds[c] = eelbrain.Factor(metadata_subset[c].values)
# Then the prediction accuracy r and TRFs.
# Because ds has only 1 row, call the add_case_dim() to add a Case dimension so that the datasests can be combined later
ds['full_r'] = add_case_dim(full_r)
ds['full_h_envelope'] = add_case_dim(full_h_envelope)
ds['full_h_onsets'] = add_case_dim(full_h_onsets)
ds['full_h_scaled_envelope'] = add_case_dim(full_h_scaled_envelope)
ds['full_h_scaled_onsets'] = add_case_dim(full_h_scaled_onsets)
ds['envelope_only_r'] = add_case_dim(envelope_only_r)
ds['onsets_delta_r'] = add_case_dim(onsets_delta_r)
ds['envelope_delta_r'] = add_case_dim(envelope_delta_r)
ds_list.append(ds)
# Combine all datasets
results = eelbrain.combine(ds_list)
print(f'Results are saved in an Eelbrain dataset with {results.n_cases} rows and {len(results)} columns:' )
print(results.head())
Results are saved in an Eelbrain dataset with 40 rows and 13 columns: # filename sub session intensity modality --------------------------------------------------------------------------- 0 sub-01_ses-2_task-60bi_acq-2.fif 1 2 60 bi 1 sub-01_ses-2_task-60mo_acq-1.fif 1 2 60 mo 2 sub-02_ses-2_task-60bi_acq-3.fif 2 2 60 bi 3 sub-02_ses-2_task-60mo_acq-4.fif 2 2 60 mo 4 sub-03_ses-2_task-60bi_acq-1.fif 3 2 60 bi 5 sub-03_ses-2_task-60mo_acq-2.fif 3 2 60 mo 6 sub-04_ses-2_task-60bi_acq-4.fif 4 2 60 bi 7 sub-04_ses-2_task-60mo_acq-3.fif 4 2 60 mo 8 sub-05_ses-2_task-60bi_acq-2.fif 5 2 60 bi 9 sub-05_ses-2_task-60mo_acq-1.fif 5 2 60 mo --------------------------------------------------------------------------- NDVars: full_r, full_h_envelope, full_h_onsets, full_h_scaled_envelope, full_h_scaled_onsets, envelope_only_r, onsets_delta_r, envelope_delta_r
Note that for an eelbrain.Dataset like results, the number of rows is given by results.n_cases. Do not use len(results) for this purpose because this returns the number of columns, including the columns for NDVar objects. This behavior is different from a Pandas DataFrame.
Check prediction accuracy and TRFs¶
Before comparing different conditions or models, first check whether participants show evidence of neural tracking of the acoustic features. Very low prediction accuracy or messy TRFs even for speech in quiet may indicate data quality or processing issues, such as noisy EEG or EEG-stimulus misalignment. To access the eelbrain.NDVar containing the prediction accuracy for the full acoustic model, use:
print(results['full_r'])
print(results['full_r'].x) # 2D Numpy array with shape (40, 32)
<NDVar 'full_r': 40 case, 32 sensor> [[0.05774414 0.0621944 0.04473119 ... 0.04176542 0.05329045 0.04373066] [0.05052079 0.06017615 0.05923666 ... 0.01592444 0.06833154 0.06632544] [0.12342843 0.09928786 0.14178462 ... 0.10233013 0.07506555 0.0928209 ] ... [0.02333837 0.03116023 0.03360338 ... 0.01928693 0.02275538 0.02582008] [0.04640964 0.04737036 0.04306908 ... 0.05139373 0.05308849 0.04245131] [0.02893636 0.03421895 0.03584145 ... 0.02786277 0.0413256 0.04483416]]
The NDVar gives you the r value for each of the 32 electrodes (sensors) and each row in the dataset. You can do topographies visualizing the distribution of prediction accuracy using plotting functions from eelbrain:
print('Grand average prediction accuracy (r):', results['full_r'].mean('case').mean('sensor'))
# Topomap of grand average prediction accuracy
p = eelbrain.plot.Topomap('full_r', data = results, clip = 'circle', cmap = 'Reds')
cbar = p.plot_colorbar(h = 1.2)
# Same plot, but divided by modality condition
p = eelbrain.plot.Topomap('full_r', 'modality', data = results, clip = 'circle', cmap = 'Reds',
axtitle = ['Binaural', 'Monaural'])
cbar = p.plot_colorbar(h = 1.2)
Grand average prediction accuracy (r): 0.04710825791123796
The overall prediction accuracy is 0.047, which looks good. The r values can vary a lot depending on several factors, but they usually lie somewhere between 0.02 and 0.1 in quiet conditions. Also, the strongest neural tracking is observed over the frontocentral electrodes, which is typical for neural tracking of envelope-based features.
Next, let’s inspect the fitted TRFs. In the code below, we first sum the fitted TRF weights across the eight frequency bands and then average them across electrodes for each acoustic predictor type. We then combine the results across the two predictor types. This provides a compact summary of the overall impulse response to the acoustic predictors, averaged across all participants.
results['avg_envlope_h'] = results['full_h_envelope'].sum('frequency').mean('sensor')
results['avg_onsets_h'] = results['full_h_onsets'].sum('frequency').mean('sensor')
results['acoustic_trf'] = results['avg_envlope_h'] + results['avg_onsets_h']
# Plot the acoustic TRF
p = eelbrain.plot.UTSStat(
y = 'acoustic_trf', data = results,
ylabel = 'TRF weights [a.u.]', axtitle = 'TRF to acoustics')
p.add_vline(x = 0.0, linewidth = 1.0, linestyle = '--', color = 'gray')
# Same plot, but divided by modality condition
p = eelbrain.plot.UTSStat(
y = 'acoustic_trf', x = 'modality', data = results, colors = ['#B51C1E', '#03749C'],
ylabel = 'TRF weights [a.u.]', axtitle = 'TRF to acoustics')
p.add_vline(x = 0.0, linewidth = 1.0, linestyle = '--', color = 'gray')
First, note that the TRF time range extends from −100 to 500 ms, corresponding to $\tau_{\min}$ and $\tau_{\max}$, respectively. Acoustic TRFs typically show a positive peak between 0 and 100 ms, followed by a negative deflection and a smaller later positive peak. The response should be relatively flat before 0 ms and after approximately 300 ms. It may be OK if your TRF morpholgy slightly deviate or doesn't have all components. But if it is very flat and appears very noisy, check whether this pattern is driven by specific participants or data files with poor quality or alignment issues.
The code above uses h from the boosting result objects rather than h_scaled. See the Eelbrain documentation here for a detailed explanation of the difference. In brief, when the model is estimated with scale_data = True (which is the default behavior and the setting used in the example results), the predictors are scaled before TRF estimation. Specifically, the mean is subtracted from each predictor (e.g., envelope at each frequency band), and the result is divided by either the standard deviation for L2 normalization or the mean absolute value for L1 normalization. The TRFs stored in h are then in the transformed predictor space and have arbitrary units. This is useful when you want to combine TRF weights across predictors to obtain a summary of the overall brain response to the acoustic features.
In contrast, h_scaled contains TRFs that have been scaled back to the original units of each predictor. This can be useful if you want to interpret TRF weights in relation to the original predictor values. However, if you plan to average or sum TRF weights across different predictor types, make sure these predictor types have the same units.
Mass-univariate tests¶
One way to compare prediction accuracy across groups or conditions is to run statistical tests (e.g., t-test, ANOVA) separately at each electrode. This can also help identify where on the scalp the effect occurs but testing each electrode separately introduces the multiple comparisons problem and can inflate the Type I error rate.
A solution is to use permutation-based mass-univariate tests. See the Eelbrain paper and documentation for details. In brief, mass-univariate tests evaluate effects across many electrodes and/or time points while controlling for multiple comparisons. Cluster-based tests do this by identifying clusters of contiguous electrodes and/or time points that show a consistent effect, and then evaluating the statistical significance of those clusters using permutation testing.
For our example results, we can run a mass-univarate paired t-test to comapre the prediction accuracy of the full model between the binaural ('bi') and monaural ('mo') conditions. We use paired/releated-measures tests here because each pariticpant provided observations in both modality conditions and hence modality is a within-subject factor (use testnd.TTestIndependent for between-subjects factors):
# Setting connectivity
results['full_r'].sensor.set_adjacency(connect_dist = 1.6)
test_res = eelbrain.testnd.TTestRelated(
'full_r', 'modality', c1 = 'bi', c0 = 'mo', data = results, pmin = 0.05,
match = 'sub', tail = 0 # Two-tailed test
)
sig_regions = test_res.find_clusters(0.05)
sig_regions
Permutation test: 100%|██████████| 10000/10000 [00:04<00:00, 2218.03 permutations/s]
| # | id | n_sensors | v | p | sig |
| 0 | 1 | 10 | 29.514 | 0.0264 | * |
The test reults tell you that there is a significant cluster with p = 0.0264 containing 10 electrodes. You can access this cluster (first row with ID = 1) and plot a map of the t values:
# Get the cluster NDVar containing the t-values
cluster_id = sig_regions[0, 'id']
cluster = test_res.cluster(cluster_id)
# Significant electrodes should have non-zero t-values
mask = cluster != 0
sig_electrodes = cluster.sensor.names[mask]
p = eelbrain.plot.Topomap(cluster, clip = 'circle', cmap = 'Reds', axtitle = 't-value map of sig. cluster')
p.mark_sensors(sig_electrodes, color = 'yellow', zorder = 2)
cbar = p.plot_colorbar(h = 1.2)
print('Electrodes in the sig. cluster:', sig_electrodes)
print('tmax =', np.max(cluster.x[mask]))
Electrodes in the sig. cluster: ['Fp1', 'F3', 'FC1', 'T8', 'FC6', 'F4', 'F8', 'AF4', 'Fp2', 'Fz'] tmax = 4.3402305568560315
Significant difference is mostly in the frontal region and the right hemisphere. You can report the largest t-value ($t_{max}$) as the effect size measure (Brodbeck et al., 2018).
Mass-univariate tests can similarly be run on TRFs, which have an additional time dimension. Here, we focus on TRFs to acoustic envelope and first average them across the frequency band dimension (note that h_scaled is used here):
results['envelope_trf_avg'] = results['full_h_scaled_envelope'].mean('frequency')
results['envelope_trf_avg'].sensor.set_adjacency(connect_dist = 1.6)
test_res = eelbrain.testnd.TTestRelated(
'envelope_trf_avg', 'modality', c1 = 'bi', c0 = 'mo', data = results, pmin = 0.05,
match = 'sub', tail = 0 # Two-tailed test
)
sig_regions = test_res.find_clusters(0.05)
sig_regions
Permutation test: 100%|██████████| 10000/10000 [00:05<00:00, 1718.98 permutations/s]
| # | id | n_sensors | tstart | tstop | duration | v | p | sig |
| 0 | 1 | 28 | 0.0015625 | 0.22031 | 0.21875 | 1032.8 | 0 | *** |
The results show that there spatiotempotal cluster containing 28 electrodes and extending from 0.002 to 0.220. You can visualize the t-values and significant cluster using the built-in Eelbrain function eelbrain.plot.TopoArray():
p = eelbrain.plot.TopoArray(
test_res, t = [0.055, 0.155], # Show the topomap at these two time points
head_radius = 0.35, cmap = 'RdBu_r', clip = 'circle',
title = test_res, axw = 3.0, axh = 4.2)
p_cb = p.plot_colorbar(label = 't-value')
Region-of-interest (ROI) analysis¶
Sometimes you may want to focus on a specific region of interest (ROI). For neural tracking of acoustics, one option is to define the ROI based on electrodes that show significantly above-zero prediction accuracy. This can be done using a mass-univariate one-sample t-test:
# Aggregate the results by subject across the two modality conditions
# This ensures that there will only one prediction accuracy observation for each subject at each electrode, avoiding artificially inflating the number of observations due to repeated measures
results_agg = results.aggregate('sub', drop_bad = True)
test_res = eelbrain.testnd.TTestOneSample(
'full_r', popmean = 0, match = 'sub', data = results_agg, tail = 1, # One-tailed test
pmin = 0.05)
test_res.find_clusters(0.05)
Permutation test: 100%|██████████| 10000/10000 [00:04<00:00, 2094.70 permutations/s]
| # | id | n_sensors | v | p | sig |
| 0 | 1 | 32 | 295.28 | 0 | *** |
# Get the cluster NDVar containing the t-values
cluster_id = sig_regions[0, 'id']
cluster = test_res.cluster(cluster_id)
# Significant electrodes should have non-zero t-values
mask = cluster != 0
sig_electrodes = cluster.sensor.names[mask]
p = eelbrain.plot.Topomap(cluster, clip = 'circle', cmap = 'Reds', axtitle = 't-value map of sig. cluster')
p.mark_sensors(sig_electrodes, color = 'yellow', zorder = 2)
cbar = p.plot_colorbar(h = 1.2)
print('Electrodes in the sig. cluster:', sig_electrodes)
print('tmax =', np.max(cluster.x[mask]))
Electrodes in the sig. cluster: ['Fp1', 'AF3', 'F7', 'F3', 'FC1', 'FC5', 'T7', 'C3', 'CP1', 'CP5', 'P7', 'P3', 'Pz', 'PO3', 'O1', 'Oz', 'O2', 'PO4', 'P4', 'P8', 'CP6', 'CP2', 'C4', 'T8', 'FC6', 'FC2', 'F4', 'F8', 'AF4', 'Fp2', 'Fz', 'Cz'] tmax = 13.448790763351411
You can see that all electrodes are included in the significant cluster, so the ROI includes the full set of electrodes.
Alternatively, if you have a predefined ROI, you can extract the average prediction accuracy within that ROI to obtain a single value for each row in the dataset. These values can then be analyzed using standard parametric methods such as mixed-effects regression in R:
# Frontocentral ROI
fronto_central_chs = ['Fz', 'F3', 'F4', 'FC1', 'FC2', 'Cz', 'C3', 'C4']
results['full_r_roi'] = results['full_r'].sub(sensor = fronto_central_chs).mean('sensor')
# Convert to Pandas dataframe and do a quick plot
results_df = results.as_dataframe()
# Color scheme
colors = {'bi': '#B51C1E', 'mo': '#03749C'}
modalities = ['bi', 'mo']
modality_labels = {'bi': 'Binaural', 'mo': 'Monaural'}
positions = [1, 2]
fig, ax = plt.subplots(figsize=(4, 3))
for pos, modality in zip(positions, modalities):
vals = results_df.loc[results_df['modality'] == modality, 'full_r_roi'].values
color = colors[modality]
# Boxplot without fill
ax.boxplot(
vals,
positions=[pos],
widths=0.5,
patch_artist=False,
showfliers=False,
boxprops=dict(color=color, linewidth=1.8),
whiskerprops=dict(color=color, linewidth=1.8),
capprops=dict(color=color, linewidth=1.8),
medianprops=dict(color=color, linewidth=1.8),
)
# Individual points with jitter
rng = np.random.default_rng(2026 + pos)
jitter = rng.uniform(-0.08, 0.08, size=len(vals))
ax.scatter(
np.full(len(vals), pos) + jitter,
vals, color=color, s=30, alpha=0.5, zorder=3
)
# Mean marker
ax.scatter(
pos, vals.mean(), color=color, s=90, marker="D", zorder=4
)
ax.set_xticks(positions)
ax.set_xticklabels([modality_labels[m] for m in modalities])
ax.set_xlabel('Modality')
ax.set_ylabel('Prediction accuracy (r)')
# Cleaner panel aesthetics
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.tight_layout()
plt.show()
Unique contribution¶
We also estimated the unique contribution of each acoustic predictor. As with overall prediction accuracy, several types of analyses can be performed on the delta r measures. For example, you can run mass-univariate one-sample t-tests to determine whether the unique contribution of a predictor differs significantly from zero. Note that the unique contribution can sometimes be negative, because adding non-informative predictors may reduce prediction accuracy.
Below is an example comparing the unique contribution of the 8-band envelope predictor across the two modality conditions. The results show two significant clusters. You can follow the steps above to identify where these clusters occur and how the conditions differ.
# Setting connectivity
results['envelope_delta_r'].sensor.set_adjacency(connect_dist = 1.6)
test_res = eelbrain.testnd.TTestRelated(
'envelope_delta_r', 'modality', c1 = 'bi', c0 = 'mo', data = results, pmin = 0.05,
match = 'sub', tail = 0 # Two-tailed test
)
sig_regions = test_res.find_clusters(0.05)
sig_regions
Permutation test: 100%|██████████| 10000/10000 [00:05<00:00, 1992.26 permutations/s]
| # | id | n_sensors | v | p | sig |
| 0 | 1 | 6 | 17.727 | 0.0336 | * |
| 1 | 7 | 7 | 24.21 | 0.0167 | * |
Conclusion¶
This tutorial covered only basic analyses of TRF results from neural tracking of acoustic features. There are many other analytic approaches and the same general workflow can also be extended to investigate neural tracking of higher-order linguistic features. Feel free to adapt and expand the exmaple code for your own EEG project, and consult the Eelbrain documentation and relevant literature when deciding which analysis approaches are most appropriate for your research questions.