Resting EEG Data

This notebook analyzes a dataset of EEG data collected from healthy young adults.

Dataset Details

This dataset is a dataset of extra-cranial eyes-closed resting-state EEG data recorded in the Voytek Lab.

# Setup notebook state
from nbutils import setup_notebook; setup_notebook()
import numpy as np
import matplotlib.pyplot as plt
from fooof import FOOOFGroup
from fooof.utils import trim_spectrum
from fooof.objs import combine_fooofs
from neurodsp.utils import create_times
from neurodsp.spectral import compute_spectrum
from neurodsp.plts import plot_time_series, plot_power_spectra
from neurodsp.plts.utils import make_axes
# Import custom project code
from apm.io import APMDB, get_files, load_pickle
from apm.io.data import load_eeg_demo_group_data, load_eeg_demo_info
from apm.analysis import (compute_avgs, unpack_corrs, compute_all_corrs,
                          compute_corrs_to_feature, compute_diffs_to_feature)
from apm.plts import plot_dots, plot_corr_matrix
from apm.plts.results import plot_topo
from apm.plts.multi import plot_results_all, plot_topo_row
from apm.plts.utils import figsaver, plot_colorbar
from apm.plts.settings import LABELS
from apm.utils import format_corr
# Import dataset settings from script
import sys; from pathlib import Path;
sys.path.append(str(Path('..').resolve() / 'scripts'))
from scripts.analyze_eeg1 import DATA_FOLDER, FS, N_SECONDS

Set Up Paths

# Define load path
db = APMDB()
LOADPATH = db.data_path / 'eeg1'
# Check the set of results files available for this dataset
get_files(LOADPATH)
['eeg1_results.p',
 'eeg1_results_peaks.p',
 'eeg1_spatial_alpha_corrs_diffs.p',
 'eeg1_spatial_corrs.p',
 'eeg1_spatial_corrs_alpha.p',
 'eeg1_spatial_corrs_exp.p',
 'specparam']

Settings

# Define times vector for data
times = create_times(N_SECONDS + 1 / FS, FS)
# Add plot kwargs
dot_kwargs = {
    'alpha' : 0.75,
}
# Settings for saving figures
SAVE_FIG = True
FIGPATH = db.figs_path / '51_eeg_data'

# Create helper function to manage figsaver settings
fsaver = figsaver(SAVE_FIG, FIGPATH)
# Define topography vlims for each measure
vlims = {
    'autocorr_decay_time' : (0.1, 1.2),
    'dfa' : (0.6, 1.5),
    'higuchi_fd' : (1.2, 1.9),
    'hjorth_complexity' : (3.5, 18),
    'lempelziv' : (180, 515),
    'sample_entropy' : (0.25, 1.15),
    'perm_entropy' : (2.4, 2.6),
    'irasa' : (0.70, 1.85),
    'specparam' : (0.70, 1.85),
}

Load Data

# Load group data
group_data = load_eeg_demo_group_data(DATA_FOLDER)
# Check data size
n_subjs, n_chs, n_times = group_data.shape
print('Number of subjects: {}'.format(n_subjs))
Number of subjects: 29
# Load MNE info object for the current dataset
info = load_eeg_demo_info(DATA_FOLDER)

Data Checks

# Set example subject index
subj_ind = 13
# Set example channel index
chi = 29
# Check example channel label
info.ch_names[chi]
'Oz'
# Plot a segment of time series data
plot_time_series(times, group_data[subj_ind, chi, :], lw=1., xlim=[5, 10], figsize=(12, 1))
plt.gca().axis('off');
if SAVE_FIG: plt.savefig(FIGPATH / ('eeg_timeseries.pdf'))
../_images/51-RestingEEGData_23_0.png
# Compute a power spectrum of an example
freqs, powers = compute_spectrum(group_data[subj_ind, chi, :], FS, nperseg=2*FS, noverlap=FS)
# Plot the power spectrum of the example data segment
plot_power_spectra(*trim_spectrum(freqs, powers, [3, 40]), minorticks=False, 
                   xticks=[5, 10, 25], xticklabels=[5, 10, 25], yticks=[],
                   figsize=(4.25, 4), **fsaver('rest_eeg_psd'))
../_images/51-RestingEEGData_25_0.png

Check Spectral Fits

# Load specparam measures
all_fgs = []
temp = FOOOFGroup()
for fres in get_files(LOADPATH / 'specparam'):
    temp.load(fres, LOADPATH / 'specparam')
    all_fgs.append(temp.copy())
fg = combine_fooofs(all_fgs)
# Check overall group results
fg.print_results()
fg.plot()
==================================================================================================
                                                                                                  
                                       FOOOF - GROUP RESULTS                                      
                                                                                                  
                            Number of power spectra in the Group: 1856                            
                                                                                                  
                        The model was run on the frequency range 3 - 40 Hz                        
                                 Frequency Resolution is 0.50 Hz                                  
                                                                                                  
                              Power spectra were fit without a knee.                              
                                                                                                  
                                      Aperiodic Fit Values:                                       
                        Exponents - Min: -1.096, Max:  2.726, Mean: 1.379                         
                                                                                                  
                        In total 9302 peaks were extracted from the group                         
                                                                                                  
                                     Goodness of fit metrics:                                     
                            R2s -  Min:  0.565, Max:  0.995, Mean: 0.972                          
                         Errors -  Min:  0.034, Max:  0.139, Mean: 0.067                          
                                                                                                  
==================================================================================================
../_images/51-RestingEEGData_28_1.png

Load Results

# Load precomputed aperiodic measure results
group_results = load_pickle('eeg1_results', LOADPATH)
# Check size of computed results [n_subjs, n_chs]
group_results['dfa'].shape
(29, 64)
# Load precomputed peak results
group_results_peaks = load_pickle('eeg1_results_peaks', LOADPATH)

Check Measures & Labels

# Check list of computed measures
print(list(group_results.keys()))
['autocorr_decay_time', 'dfa', 'higuchi_fd', 'hjorth_complexity', 'lempelziv', 'sample_entropy', 'perm_entropy', 'irasa', 'specparam']
# Collect list of exponent & timeseries measure labels
exp_measures = ['specparam', 'irasa']
ts_measures = list(group_results.keys())
[ts_measures.remove(meas) for meas in exp_measures];
# Collect labels for time series measures
ts_labels = [LABELS[meas] for meas in ts_measures]

Extract Measures for Example Channel

# Sub-select results to channel of interest
results = {key : val[:, chi] for key, val in group_results.items()}
results_peaks = {key : val[:, chi] for key, val in group_results_peaks.items()}
# Update missing value for AC value at Cz
results['autocorr_decay_time'][9] = np.mean([group_results['autocorr_decay_time'][9, 28],
                                             group_results['autocorr_decay_time'][9, 30]])
# Compute correlations for selected channel
all_corrs = compute_all_corrs(results)

Compare Exponent Measures

# Plot the comparison of specparam and IRASA exponent estimations
plot_dots(results['specparam'], results['irasa'], **dot_kwargs, figsize=(5, 4),
          xlim=[0.5, 2.5], ylim=[0.5, 2.5], tposition='tl', expected=[0, 3],
          xlabel='Aperiodic Exponent (SP)', ylabel='Aperiodic Exponent (IR)',
          **fsaver('eeg_exp_exp_comp'))
../_images/51-RestingEEGData_42_0.png
# Check correlation between specparam and irasa exponent estimates
print('  SP-EXP & IR-EXP:  ', format_corr(*all_corrs['specparam']['irasa']))
  SP-EXP & IR-EXP:   r=+0.946  CI[+0.853, +0.979],  p=0.000

Compare Exponent to Time Series Measures

# Plot comparisons between exponent and time series measures
axes = make_axes(1, len(ts_measures), figsize=(30, 4), wspace=0.05)
for ind, meas in enumerate(ts_measures):
    plot_dots(results['specparam'], results[meas], **dot_kwargs,
              xlabel='Aperiodic Exponent', ylabel=LABELS[meas], ax=axes[ind])
axes[0].set_ylim([0, 1])
if SAVE_FIG: plt.savefig(FIGPATH / ('eeg1_exp_ts_scatters.pdf'))
../_images/51-RestingEEGData_45_0.png
# Check the correlations between time series and exponent measures
for meas in ts_measures:
    print(meas)
    print('    SP-EXP:  ', format_corr(*all_corrs['specparam'][meas]))
    print('    IR-EXP:  ', format_corr(*all_corrs['irasa'][meas]))
autocorr_decay_time
    SP-EXP:   r=+0.261  CI[-0.131, +0.600],  p=0.172
    IR-EXP:   r=+0.130  CI[-0.265, +0.480],  p=0.501
dfa
    SP-EXP:   r=+0.017  CI[-0.384, +0.435],  p=0.931
    IR-EXP:   r=-0.098  CI[-0.480, +0.315],  p=0.613
higuchi_fd
    SP-EXP:   r=-0.517  CI[-0.770, -0.159],  p=0.004
    IR-EXP:   r=-0.638  CI[-0.801, -0.372],  p=0.000
hjorth_complexity
    SP-EXP:   r=+0.723  CI[+0.463, +0.878],  p=0.000
    IR-EXP:   r=+0.693  CI[+0.424, +0.847],  p=0.000
lempelziv
    SP-EXP:   r=-0.605  CI[-0.807, -0.298],  p=0.001
    IR-EXP:   r=-0.572  CI[-0.767, -0.241],  p=0.001
sample_entropy
    SP-EXP:   r=-0.628  CI[-0.827, -0.315],  p=0.000
    IR-EXP:   r=-0.636  CI[-0.800, -0.357],  p=0.000
perm_entropy
    SP-EXP:   r=-0.190  CI[-0.468, +0.135],  p=0.323
    IR-EXP:   r=-0.309  CI[-0.554, -0.002],  p=0.102

Compare Time Series Measures to Each Other

# Plot multi-panel plot comparing all time series measures to each other
plot_results_all(results, ts_measures)
../_images/51-RestingEEGData_48_0.png

Correlations

# Subselect time domain measures
all_corrs_ts = {ke : va for ke, va in all_corrs.items() if ke not in exp_measures}
# Plot the correlations matrix across all time series measures
plot_corr_matrix(unpack_corrs(all_corrs_ts), cbar=False, figsize=(5, 5),
                 xticklabels=ts_labels, yticklabels=ts_labels, **fsaver('subj_corrs_ts'))
../_images/51-RestingEEGData_51_0.png
# Extract the correlations between specparam and time domain measures
exp_corrs_subjs = np.atleast_2d([all_corrs['specparam'][label][0] for label in ts_measures]).T
# Plot correlations between exponent and time domain measures
plot_corr_matrix(exp_corrs_subjs, cbar=False, figsize=(1.5, 5), **fsaver('subj_corrs_exp'))
../_images/51-RestingEEGData_53_0.png

Compare to Alpha Power

# Compute correlations between measures at selected electrode and alpha
alpha_corrs = compute_corrs_to_feature(results, results_peaks['alpha_power'])
# # Compute differences between correlations to alpha
# alpha_corr_diffs = compute_diffs_to_feature(results, results_peaks['alpha_power'])
# Check the correlations between alpha power and aperiodic measures
print('Correlations with alpha:')
for label in alpha_corrs.keys():
    print('     {:20s}:  '.format(label), format_corr(*alpha_corrs[label]))
Correlations with alpha:
     autocorr_decay_time :   r=-0.464  CI[-0.718, -0.111],  p=0.011
     dfa                 :   r=-0.636  CI[-0.808, -0.362],  p=0.000
     higuchi_fd          :   r=-0.635  CI[-0.832, -0.334],  p=0.000
     hjorth_complexity   :   r=+0.317  CI[-0.084, +0.630],  p=0.094
     lempelziv           :   r=-0.324  CI[-0.606, +0.064],  p=0.086
     sample_entropy      :   r=-0.388  CI[-0.671, +0.014],  p=0.038
     perm_entropy        :   r=-0.595  CI[-0.807, -0.276],  p=0.001
     irasa               :   r=+0.309  CI[-0.062, +0.628],  p=0.103
     specparam           :   r=+0.287  CI[-0.082, +0.592],  p=0.132
# Organize correlations between alpha and time domain measures
alpha_corrs_ts = np.atleast_2d([alpha_corrs[label][0] for label in ts_measures]).T
# Plot correlations between alpha and time domain measures
plot_corr_matrix(alpha_corrs_ts, cbar=False, figsize=(1.5, 5), **fsaver('subj_corrs_alpha'))
../_images/51-RestingEEGData_59_0.png

Spatial Analyses

Compute Measures Across Channels

# Compute the average across the group
group_avg = compute_avgs(group_results, np.nanmean)

Exponent Topograghies

# Plot the group average topographies for the exponent measures
plot_topo_row(group_avg, exp_measures, info, vlims, **fsaver('exp_topos'))
../_images/51-RestingEEGData_64_0.png
# Create colorbar for the aperiodic exponent topographies
plot_colorbar('Aperiodic Exponent', *vlims['specparam'], 
              **fsaver('colorbar_eeg1_exp'), close=True)
# Compute correlation between exponent topographies
exp_spatial_corr = compute_all_corrs({meas: group_avg[meas] for meas in exp_measures})
# Check correlation between specparam and irasa exponent estimates
print('  SP-EXP & IR-EXP:  ', format_corr(*exp_spatial_corr['specparam']['irasa']))
  SP-EXP & IR-EXP:   r=+0.972  CI[+0.938, +0.985],  p=0.000

Time Series Measure Topographies

# Plot the group average topographies for the time domain measures
plot_topo_row(group_avg, ts_measures, info, vlims, **fsaver('eeg1_topo_row'))
../_images/51-RestingEEGData_69_0.png

Spatial Correlations

# Load precomputed group correlation results
group_corrs = load_pickle('eeg1_spatial_corrs', LOADPATH)
# Plot the correlation matrix of spatial topographies
plot_corr_matrix(unpack_corrs(group_corrs), cbar=False, figsize=(5, 5),
                 xticklabels=ts_labels, yticklabels=ts_labels, **fsaver('space_corrs_ts'))
../_images/51-RestingEEGData_72_0.png
# Load precomputed group exponent correlation results
group_exp_corrs = load_pickle('eeg1_spatial_corrs_exp', LOADPATH)
# Organize the correlations between the exponent and time domain measures
exp_corr_mat = np.atleast_2d([group_exp_corrs[label][0] for label in group_exp_corrs]).T
# Plot the correlations between exponent and time domain measures
plot_corr_matrix(exp_corr_mat, cbar=False, figsize=(1.5, 5), **fsaver('space_corrs_exp'))
../_images/51-RestingEEGData_75_0.png

Alpha Power Topography

# Compute the average alpha topography across the group
group_avg_peaks = compute_avgs(group_results_peaks, np.nanmean)
# Check range of alpha power values
amin, amax = np.min(group_avg_peaks['alpha_power']), np.max(group_avg_peaks['alpha_power'])
amin, amax
(-11.505512552707579, -10.321239762330212)
# Set alpha vlimits
alpha_vlim = (-11.5, -10.5)
# Plot the average alpha topography across the group
plot_topo(group_avg_peaks['alpha_power'], info, vlim=alpha_vlim, **fsaver('alpha_topo'))
../_images/51-RestingEEGData_80_0.png
# Create colorbar for the aperiodic exponent topographies
plot_colorbar('Alpha Power (uV**2)', *alpha_vlim, figsize=(2.65, 7),
              **fsaver('colorbar_eeg1_alpha'), close=True)
# Load precomputed correlations between aperiodic measures and alpha power
group_alpha_corrs = load_pickle('eeg1_spatial_corrs_alpha', LOADPATH)
# Organize the alpha correlations
alpha_corr_mat = np.atleast_2d([group_alpha_corrs[label][0] for label in group_alpha_corrs]).T
# Plot the correlations between alpha activity and time domain measures
plot_corr_matrix(alpha_corr_mat, cbar=False, figsize=(1.5, 5), **fsaver('space_corrs_alpha'))
../_images/51-RestingEEGData_84_0.png

Conclusions

Conclusions of this empirical data analysis thus far:

  • in this EEG data, the specparam & IRASA exponent estimates are highly comparable

  • in this empirical data, aperiodic exponent is moderately correlated with the time domain methods

Note that this dataset is analyzed as a small / pilot, such that these interim conclusions are not considered decisive, and are best considered as initial results to be further explored and replicated in subsequent larger datasets.