Resting EEG Data
Contents
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'))
# 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'))
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
==================================================================================================
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'))
# 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'))
# 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)
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'))
# 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'))
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'))
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'))
# 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'))
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'))
# 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'))
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'))
# 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'))
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.