Long-term reliability and stability of parameterized resting state EEG: evidence from a five-year follow-up.

Code ↔ Paper

4 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 4 matches

[1] § Method › Physiological recording and data reduction ↔ notebooks/IfAdo_preprocessing.ipynb, lines 86–149 · score 0.95 · Power Spectral Density, peak width limits, 3–40 Hz, peak height, peak threshold, aperiodic mode
[2] § Method › Physiological recording and data reduction ↔ src/preprocessing.py, lines 195–312 · score 0.94 · Power Spectral Density, peak width limits, peak height, peak threshold, aperiodic mode, 1–8 Hz
[3] § Method › Physiological recording and data reduction ↔ src/preprocessing.py, lines 195–312 · score 0.63 · MNE, segmented, rejected, overlapping, epochs, filtered
[4] § Method › Physiological recording and data reduction ↔ notebooks/IfAdo_preprocessing.ipynb, lines 86–149 · score 0.57 · Bad channels, rejected, artifact, overlapping, epochs, filtered

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Jupyter notebook · 1,986 lines · 168 KB · MIT · 2 matches

# %% [markdown]
# Long-term reliability and stability of parameterised resting state EEG
# ======================================================================
#
# Analysis pipeline for 5-year test-retest reliability study.
#
# Data: Dortmund Vital Study (Getzmann et al., 2024)
# https://doi.org/10.18112/openneuro.ds005385.v1.0.2
#
# Measures extracted:
# - Aperiodic exponent
# - Aperiodic offset
# - Individual alpha peak frequency (IAPF)
# - Alpha power (via FOOOF)
# - Absolute Alpha power
#
# Repository: https://github.com/MindSpaceLab/Aperiodic_Test_Retest_5year
# %%
# Core Imports
import os
import numpy as np
import pandas as pd
import mne as mne
import scipy as scipy
import matplotlib.pyplot as plt
import seaborn as sns
import time
from mne_icalabel import label_components
from mne_bids import (
BIDSPath,
find_matching_paths,
get_entity_vals,
make_report,
print_dir_tree,
read_raw_bids,
)
from pyprep.find_noisy_channels import NoisyChannels
import mne_bids as mne_bids
import neurodsp as ndsp
import fooof as fooof
from pathlib import Path
from sklearn.decomposition import PCA
from sklearn.metrics import silhouette_score
from sklearn.cluster import KMeans
from scipy.spatial import ConvexHull
from scipy.stats import pearsonr
import stargazer
#from statsmodels.stats.multitest import multipletests
import warnings
warnings.filterwarnings('ignore')
# EEG processing
mne.set_log_level('WARNING')
# Add custom modules path
import sys
sys.path.insert(0, str(Path.cwd() / 'src'))
print("✓ Core libraries loaded")
%matplotlib qt
# %%
# Custom modules
import sys
import importlib
sys.path.append("../src")
module_name = ["preprocessing", 'utils','viz']
for mod in module_name:
if mod in sys.modules:
del sys.modules[mod]
try:
import preprocessing as preproc
import utils as utils
import viz as viz
print("✓ Custom modules loaded")
except ImportError as e:
print(f"⚠ Custom modules not found: {e}")
print(" Create modules in src/ directory as needed")
# %%
# Configuration
DATA_PATH = Path(os.path.realpath(
os.path.join(os.path.abspath(''), '..', 'Data')))
FIGURES_PATH = Path(os.path.realpath(
os.path.join(os.path.abspath(''), '..', 'Figures')))
RESULTS_PATH = Path(os.path.realpath(
os.path.join(os.path.abspath(''), '..', 'Results')))
# Create directories if they don't exist
FIGURES_PATH.mkdir(parents=True, exist_ok=True)
RESULTS_PATH.mkdir(parents=True, exist_ok=True)
# Preprocessing parameters
FILTER_LOW = 0.1
FILTER_HIGH = 40
EPOCH_DUR = 2.0 # seconds
SAMPLE_RATE = 250
REFERENCE = ['TP9', 'TP10']
# Artifact rejection parameters
THRESHOLD = 150e-6 # 100 µV
# Parameters for PSD analysis
BANDWIDTH = 2 # Multitaper bandwidth
# FOOOF parameter settings
SPEC_PARAM_SETTINGS = {
'freq_range': [3,40], # Frequency range to fit
'peak_width_limits': [1, 8], # Frequency resolution * 2 recommended
'min_peak_height': 0.1,
#'max_n_peaks': 6,
'peak_threshold': 2,
'aperiodic_mode': 'fixed', # 'fixed' or 'knee'
}
N_JOBS = -1 # Number of jobs to run in parallel
PSD_SETTINGS = {
'filter_low': 0.1,
'filter_high': 40,
'epoch_duration': 2.0,
'epoch_overlap': 0,
'reject_threshold': 200e-6, # 200 µV peak-to-peak, equivalent to ±100 µV
'fmin': 1,
'fmax': 40.0,
'n_fft': 1000,
'n_per_seg': 500,
'n_overlap': 250,
'window': 'hamming'
}
#Pre-ICC and HLM parameters
FIT_THRESHOLD = 0.9 # Minimum fit quality for FOOOF model to be included
BAD_CHAN_THRESHOLD = 0.5 # Maximum proportion of bad channels allowed per subject
print(f"✓ Config set - Data: {DATA_PATH}, Filter: {FILTER_LOW}-{FILTER_HIGH} Hz, Reference: {REFERENCE}, Epoch: {EPOCH_DUR} s, "
f"Sample rate: {SAMPLE_RATE} Hz, Artifact threshold: {THRESHOLD}, "
f"PSD bandwidth: {BANDWIDTH}, "
f"FOOOF settings: {SPEC_PARAM_SETTINGS},"
f"Number of jobs: {N_JOBS},"
f"PSD settings: {PSD_SETTINGS},"
f"Pre-ICC and HLM settings: {{'fit_threshold': {FIT_THRESHOLD}, 'bad_chan_threshold': {BAD_CHAN_THRESHOLD}}}"
)
# %%
#Although the data are stored in BIDs format, there are some issues with importing using MNE_BIDS to load data coming from EEGLAB - the formatting of the events causes some issues.
#But, we'll still use it for managing the files.
sessions = get_entity_vals(DATA_PATH, "session", ignore_sessions="off")
datatype = "eeg" #Define the datatype as EEG
extensions = [".edf"] # ignore .json files
bids_paths = find_matching_paths(
DATA_PATH, datatypes=datatype, sessions=sessions, extensions=extensions
)
subject_ids = get_entity_vals(DATA_PATH, "subject", ignore_subjects="on")
#here we manually remove subject 485, as one of the ransac bad channels is a reference channel (Tp9). If you want to verify this, you can run the preprocessing on this subject alone.
subject_ids.remove('485')
tasks = ['EyesClosed', 'EyesOpen']
sessions = ['1', '2']
acquisitions = ['pre'] #,'post'] We only want pre-task data.
task_combinations = list(enumerate([(task, ses, acq) for task in tasks for ses in sessions for acq in acquisitions]))
print(f"✓ Found {len(subject_ids)} subjects with {len(bids_paths)} files in BIDS format and {len(task_combinations)} task combinations")
for idx, (task, ses, acq) in task_combinations:
print(f" {idx}: Task={task}, Session={ses}, Acquisition={acq}")
# %% [markdown]
# # 1. Perform initial preprocessing of the IfAdo dataset.
# This includes resampling, re-referencing, RANSAC, ICA decomposition and removal of artifactual components.
# %% [markdown]
# ## 1.1 Preprocessing
# %%
#Comment out for now to avoid re-running full preprocessing on all subjects
ica_removed_df, channels_removed_df, failed = preproc.process_data_bids(subject_ids, DATA_PATH, task_combinations, sample_rate=SAMPLE_RATE, reference=REFERENCE, overwrite=True, skip_existing=False)
# %% [markdown]
# ## 1.2 Epoching and FOOOOF
# %%
#Comment out so it doesn't re-run
# Epoching and FOOOF processing - uncomment to run
failed_files, epoch_nums = preproc.epoch_psd_fooof(subject_ids, DATA_PATH, task_combinations=task_combinations, psd_settings=PSD_SETTINGS, fooof_settings=SPEC_PARAM_SETTINGS, overwrite=True)
# Save failed preprocessing subjects
if failed_files:
failed_df = pd.DataFrame(failed_files, columns=['subject'])
failed_df.to_csv(DATA_PATH / 'failed_fooof.csv', index=False)
print(f"✓ Saved failed fooof subjects to {DATA_PATH / 'failed_fooof.csv'}")
# Save failed files as a csv
failed_files_df = pd.DataFrame({'Failed Files': failed_files})
failed_files_df.to_csv('failed_files.csv', index=False)
# Save the epoch numbers as a csv
epoch_nums.to_csv('epoch_numbers.csv', index=False)
# Function to compare failed files with scanned files, just in case the failures were due to reasons other than missing files
def compare_failed_files(failed_files_df, all_files,filepath):
missing_files = []
other_failures = []
for failed_file in failed_files_df['Failed Files']:
raw_file = failed_file.replace('_proc-preproc_eeg.fif', '.edf')
raw_file_path = os.path.join(filepath, raw_file)
if raw_file_path not in all_files:
missing_files.append(failed_file)
else:
other_failures.append(failed_file)
return missing_files, other_failures
# Scan the data directory
all_files = utils.scan_data_directory(DATA_PATH)
#Load the failed files
failed_files_df = pd.read_csv('failed_files.csv')
# Compare the failed files to the scanned files
missing_files, other_failures = compare_failed_files(failed_files_df, all_files,DATA_PATH)
# %% [markdown]
# # 2. Extract FOOOF parameters from each file.
# %%
import json
from fooof.bands import Bands
bands = Bands({'alpha': [8, 13]})
from fooof.analysis import get_band_peak_fg
from fooof.utils import trim_spectrum
def extract_fooof_params(subject_ids, filepath):
fooof_params_dict = {}
for task in task_combinations:
task_name = task[1][0]
session = task[1][1]
acquisition = task[1][2]
# Initialize empty arrays for the outputs at the task level
exponents = {}
offsets = {}
fits = {}
errors = {}
alpha_freq = {}
alpha_amp = {}
alpha_absolute = {}
for subject in subject_ids:
try:
bids_path = BIDSPath(subject=subject, session=session, task=task_name, acquisition=acquisition, datatype=datatype, root=filepath, suffix='eeg', extension='.fif', processing='preproc')
fooof_path = str(bids_path.fpath)
fooof_path = fooof_path.replace('proc-preproc_eeg.fif', 'FOOOF.json')
fg = fooof.FOOOFGroup()
fg.load(fooof_path)
subject = 'sub-' + subject
exponents[subject] = fg.get_params('aperiodic_params', 'exponent')
offsets[subject] = fg.get_params('aperiodic_params', 'offset')
fits[subject] = fg.get_params('r_squared')
errors[subject] = fg.get_params('error')
alpha_fooof = get_band_peak_fg(fg, bands.alpha)
alpha_freq[subject] = np.array(alpha_fooof)[:,0]
alpha_amp[subject] = np.array(alpha_fooof)[:,1]
af, asp = trim_spectrum(fg.freqs, fg.power_spectra, bands.alpha)
alpha_absolute[subject] = np.log10(np.trapz(10**asp, af, axis=1))
except Exception as e:
print(f'Error processing {bids_path.basename}: {e}')
fooof_params_dict[f'ses-{session}_task-{task_name}_acq-{acquisition}'] = {
'exponents': exponents,
'offsets': offsets,
'fits': fits,
'errors': errors,
'alpha_freq': alpha_freq,
'alpha_amp': alpha_amp,
'alpha_absolute': alpha_absolute
}
return fooof_params_dict
# %%
# Extracts FOOOF parameters for all subjects and tasks, saving to the data path.
# Note that it will return an error if any subjects are missing data for any task, which most will be, since only ~200 have session 2 data
fooof_params_dict_master = extract_fooof_params(subject_ids, DATA_PATH)
#Splits the dictionary into separate dictionaries for each task
for task in task_combinations:
task_name = task[1][0]
session = task[1][1]
acquisition = task[1][2]
fooof_params_dict = fooof_params_dict_master[f'ses-{session}_task-{task_name}_acq-{acquisition}']
np.save(RESULTS_PATH / f'fooof_params_dict_{task_name}_{session}_{acquisition}.npy', fooof_params_dict)
# %% [markdown]
# ## 3. Load FOOOF parameters, demographics, and start cleaning data.
# %%
#Load the demographic data - participants.tsv
participants = pd.read_csv(os.path.join(DATA_PATH,'participants.tsv'),delimiter='\t')
#Because I had to do the main processing in stages, there are duplicate csvs recording the removed channels and ICA components. We need to merge these into a single dataframe.
channels_files = utils.scan_data_directory(os.path.join(os.path.abspath('..'),'notebooks'))
channels_files = [file for file in channels_files if 'channels_removed' in file]
channels_removed = pd.concat([pd.read_csv(file) for file in channels_files])
channels_removed['ID'] = channels_removed['ID'].str.split('_')
channels_removed['Task'] = channels_removed['ID'].str[2]
channels_removed['Session'] = channels_removed['ID'].str[1]
channels_removed['Acquisition'] = channels_removed['ID'].str[3]
channels_removed['ID'] = channels_removed['ID'].str[0]
channels_removed = channels_removed.pivot_table(index='ID',columns=['Task','Session','Acquisition'],values='Channels Removed',aggfunc='mean')
channels_removed.columns = ['_'.join(col).strip() for col in channels_removed.columns.values]
channels_removed.columns = ['Channels_Removed_' + col for col in channels_removed.columns]
ica_files = utils.scan_data_directory(os.path.join(os.path.abspath('..'),'notebooks'))
ica_files = [file for file in ica_files if 'ica_removed' in file]
ica_removed = pd.concat([pd.read_csv(file) for file in ica_files])
ica_removed['ID'] = ica_removed['ID'].str.split('_')
ica_removed['Task'] = ica_removed['ID'].str[2]
ica_removed['Session'] = ica_removed['ID'].str[1]
ica_removed['Acquisition'] = ica_removed['ID'].str[3]
ica_removed['ID'] = ica_removed['ID'].str[0]
ica_removed = ica_removed.pivot_table(index='ID',columns=['Task','Session','Acquisition'],values='ICAs Removed',aggfunc='mean')
ica_removed.columns = ['_'.join(col).strip() for col in ica_removed.columns.values]
ica_removed.columns = ['ICA_Removed_' + col for col in ica_removed.columns]
demographic_data = participants.merge(channels_removed, left_on='participant_id', right_on='ID')
demographic_data = demographic_data.merge(ica_removed, left_on='participant_id', right_on='ID')
# Load epoch counts
epoch_nums = pd.read_csv('epoch_numbers.csv')
fname_col = epoch_nums.columns[0]
epoch_nums_parsed = epoch_nums.copy()
epoch_nums_parsed[['subject', 'session', 'task_acq', 'rest']] = (
epoch_nums_parsed[fname_col].astype(str).str.split('_', n=3, expand=True)
)
epoch_nums_parsed['task'] = epoch_nums_parsed['task_acq'].str.replace('task-', '', regex=False)
epoch_nums_parsed['session'] = epoch_nums_parsed['session'].str.replace('ses-', '', regex=False)
epoch_nums_parsed['acquisition'] = epoch_nums_parsed['rest'].str.extract(r'acq-([A-Za-z0-9]+)')
epoch_nums_parsed['task_session_acq'] = (
epoch_nums_parsed['task'] + '_' + epoch_nums_parsed['session'] + '_' + epoch_nums_parsed['acquisition']
)
needed_cols = ['subject', 'task_session_acq', 'Original_N', 'Retained_N']
m = epoch_nums_parsed.melt(
id_vars=['subject', 'task_session_acq'],
value_vars=['Original_N', 'Retained_N'],
var_name='metric', value_name='N'
)
m['metric_clean'] = m['metric'].str.replace('_N', '', regex=False) # 'Original' or 'Retained'
m['wide_col'] = 'Epochs' + m['metric_clean'] + '_' + m['task_session_acq']
epoch_nums_wide = (
m.pivot_table(index='subject', columns='wide_col', values='N', aggfunc='first')
.reset_index()
.rename(columns={'subject': 'participant_id'})
)
demographic_data = demographic_data.merge(epoch_nums_wide, on='participant_id', how='left')
for task in task_combinations:
task_name = task[1][0]
session = task[1][1]
acquisition = task[1][2]
retained_col = f'EpochsRetained_{task_name}_{session}_{acquisition}'
original_col = f'EpochsOriginal_{task_name}_{session}_{acquisition}'
demographic_data[f'EpochsProportion_{task_name}_{session}_{acquisition}'] = demographic_data[retained_col] / demographic_data[original_col]
#Here, we create some simple flags for later filtering.
# Eyes Closed: >50% retention in BOTH session 1 and session 2 pre tasks
demographic_data['Good_Retention_EyesClosed'] = (
(demographic_data['session2'] == 'yes') &
(demographic_data['EpochsProportion_EyesClosed_1_pre'] > 0.5) &
(demographic_data['EpochsProportion_EyesClosed_2_pre'] > 0.5)
)
# Eyes Open: >50% retention in BOTH session 1 and session 2 pre tasks
demographic_data['Good_Retention_EyesOpen'] = (
(demographic_data['session2'] == 'yes') &
(demographic_data['EpochsProportion_EyesOpen_1_pre'] > 0.5) &
(demographic_data['EpochsProportion_EyesOpen_2_pre'] > 0.5)
)
# Print summary
print("Summary of good retention (>50% epochs retained in both sessions):")
print(f"Eyes Closed - Good retention: {demographic_data['Good_Retention_EyesClosed'].sum()} participants")
print(f"Eyes Open - Good retention: {demographic_data['Good_Retention_EyesOpen'].sum()} participants")
print(f"\nParticipants with session2='yes': {(demographic_data['session2'] == 'yes').sum()}")
print(f"Eyes Closed good retention (session2=yes): {demographic_data['Good_Retention_EyesClosed'].sum()} participants")
print(f"Eyes Open good retention (session2=yes): {demographic_data['Good_Retention_EyesOpen'].sum()} participants")
print("Participants with >50% channels removed in any task/timepoint:")
print(demographic_data[((demographic_data['Channels_Removed_task-EyesClosed_ses-1_acq-pre'] > 31) & (demographic_data['Channels_Removed_task-EyesClosed_ses-1_acq-pre'] < 64)) | ((demographic_data['Channels_Removed_task-EyesClosed_ses-2_acq-pre'] > 31) & (demographic_data['Channels_Removed_task-EyesClosed_ses-2_acq-pre'] < 64)) | ((demographic_data['Channels_Removed_task-EyesOpen_ses-1_acq-pre'] > 31) & (demographic_data['Channels_Removed_task-EyesOpen_ses-1_acq-pre'] < 64)) | ((demographic_data['Channels_Removed_task-EyesOpen_ses-2_acq-pre'] > 31) & (demographic_data['Channels_Removed_task-EyesOpen_ses-2_acq-pre'] < 64))]['participant_id'].tolist())
demographic_data['Good_Ref_EyesClosed'] = (
(demographic_data['Channels_Removed_task-EyesClosed_ses-1_acq-pre'] < 99) &
(demographic_data['Channels_Removed_task-EyesClosed_ses-2_acq-pre'] < 99))
demographic_data['Good_Ref_EyesOpen'] = (
(demographic_data['Channels_Removed_task-EyesOpen_ses-1_acq-pre'] < 99) &
(demographic_data['Channels_Removed_task-EyesOpen_ses-2_acq-pre'] < 99)
)
print("Summary of good reference (no bad reference channels at both time points):")
print(f"Eyes Closed - Good reference: {demographic_data['Good_Ref_EyesClosed'].sum()} participants")
print(f"Eyes Open - Good reference: {demographic_data['Good_Ref_EyesOpen'].sum()} participants")
# %%
#For each task, load the FOOOF results into its own dictionary
fooof_params_dict_EyesClosed_1_pre = np.load(RESULTS_PATH / 'fooof_params_dict_EyesClosed_1_pre.npy', allow_pickle=True)[()]
fooof_params_dict_EyesClosed_1_post = np.load(RESULTS_PATH / 'fooof_params_dict_EyesClosed_1_post.npy', allow_pickle=True)[()]
fooof_params_dict_EyesClosed_2_pre = np.load(RESULTS_PATH / 'fooof_params_dict_EyesClosed_2_pre.npy', allow_pickle=True)[()]
fooof_params_dict_EyesClosed_2_post = np.load(RESULTS_PATH / 'fooof_params_dict_EyesClosed_2_post.npy', allow_pickle=True)[()]
fooof_params_dict_EyesOpen_1_pre = np.load(RESULTS_PATH / 'fooof_params_dict_EyesOpen_1_pre.npy', allow_pickle=True)[()]
fooof_params_dict_EyesOpen_1_post = np.load(RESULTS_PATH / 'fooof_params_dict_EyesOpen_1_post.npy', allow_pickle=True)[()]
fooof_params_dict_EyesOpen_2_pre = np.load(RESULTS_PATH / 'fooof_params_dict_EyesOpen_2_pre.npy', allow_pickle=True)[()]
fooof_params_dict_EyesOpen_2_post = np.load(RESULTS_PATH / 'fooof_params_dict_EyesOpen_2_post.npy', allow_pickle=True)[()]
# %%
#Generate and save the good_mask for each task dictonary
from utils import analyze_bad_fits, analyze_missing_alpha
fooof_params_dict_EyesClosed_1_pre['good_mask'] = analyze_bad_fits(fooof_params_dict_EyesClosed_1_pre, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesClosed_1_post['good_mask'] = analyze_bad_fits(fooof_params_dict_EyesClosed_1_post, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesClosed_2_pre['good_mask'] = analyze_bad_fits(fooof_params_dict_EyesClosed_2_pre, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesClosed_2_post['good_mask'] = analyze_bad_fits(fooof_params_dict_EyesClosed_2_post, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesOpen_1_pre['good_mask'] = analyze_bad_fits(fooof_params_dict_EyesOpen_1_pre, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesOpen_1_post['good_mask'] = analyze_bad_fits(fooof_params_dict_EyesOpen_1_post, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesOpen_2_pre['good_mask'] = analyze_bad_fits(fooof_params_dict_EyesOpen_2_pre, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesOpen_2_post['good_mask'] = analyze_bad_fits(fooof_params_dict_EyesOpen_2_post, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
#And generate and save the nans for each task dictonary
fooof_params_dict_EyesClosed_1_pre['nans'] = analyze_bad_fits(fooof_params_dict_EyesClosed_1_pre, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
fooof_params_dict_EyesClosed_1_post['nans'] = analyze_bad_fits(fooof_params_dict_EyesClosed_1_post, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
fooof_params_dict_EyesClosed_2_pre['nans'] = analyze_bad_fits(fooof_params_dict_EyesClosed_2_pre, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
fooof_params_dict_EyesClosed_2_post['nans'] = analyze_bad_fits(fooof_params_dict_EyesClosed_2_post, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
fooof_params_dict_EyesOpen_1_pre['nans'] = analyze_bad_fits(fooof_params_dict_EyesOpen_1_pre, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
fooof_params_dict_EyesOpen_1_post['nans'] = analyze_bad_fits(fooof_params_dict_EyesOpen_1_post, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
fooof_params_dict_EyesOpen_2_pre['nans'] = analyze_bad_fits(fooof_params_dict_EyesOpen_2_pre, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
fooof_params_dict_EyesOpen_2_post['nans'] = analyze_bad_fits(fooof_params_dict_EyesOpen_2_post, fit_threshold=FIT_THRESHOLD, bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
#Generate and save the good_mask for each task dictonary based on missing alpha peaks
fooof_params_dict_EyesClosed_1_pre['good_mask_alpha'] = analyze_missing_alpha(fooof_params_dict_EyesClosed_1_pre, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesClosed_1_post['good_mask_alpha'] = analyze_missing_alpha(fooof_params_dict_EyesClosed_1_post, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesClosed_2_pre['good_mask_alpha'] = analyze_missing_alpha(fooof_params_dict_EyesClosed_2_pre, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesClosed_2_post['good_mask_alpha'] = analyze_missing_alpha(fooof_params_dict_EyesClosed_2_post, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesOpen_1_pre['good_mask_alpha'] = analyze_missing_alpha(fooof_params_dict_EyesOpen_1_pre, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesOpen_1_post['good_mask_alpha'] = analyze_missing_alpha(fooof_params_dict_EyesOpen_1_post, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesOpen_2_pre['good_mask_alpha'] = analyze_missing_alpha(fooof_params_dict_EyesOpen_2_pre, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesOpen_2_post['good_mask_alpha'] = analyze_missing_alpha(fooof_params_dict_EyesOpen_2_post, bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
# %%
#Generated and save good_mask based on exponent values being >0. If they're <0, we replace them with NaN.
from utils import analyze_bad_exponents
fooof_params_dict_EyesClosed_1_pre['good_mask_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesClosed_1_pre,bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesClosed_1_post['good_mask_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesClosed_1_post,bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesClosed_2_pre['good_mask_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesClosed_2_pre,bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesClosed_2_post['good_mask_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesClosed_2_post,bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesOpen_1_pre['good_mask_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesOpen_1_pre,bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesOpen_1_post['good_mask_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesOpen_1_post,bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesOpen_2_pre['good_mask_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesOpen_2_pre,bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
fooof_params_dict_EyesOpen_2_post['good_mask_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesOpen_2_post,bad_channel_prop=BAD_CHAN_THRESHOLD)['good_mask']
#and generate nans for bad exponents - but not overwriting existing nans, but adding to them.
fooof_params_dict_EyesClosed_1_pre['nans_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesClosed_1_pre,bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
fooof_params_dict_EyesClosed_1_post['nans_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesClosed_1_post,bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
fooof_params_dict_EyesClosed_2_pre['nans_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesClosed_2_pre,bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
fooof_params_dict_EyesClosed_2_post['nans_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesClosed_2_post,bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
fooof_params_dict_EyesOpen_1_pre['nans_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesOpen_1_pre,bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
fooof_params_dict_EyesOpen_1_post['nans_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesOpen_1_post,bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
fooof_params_dict_EyesOpen_2_pre['nans_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesOpen_2_pre,bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
fooof_params_dict_EyesOpen_2_post['nans_exponent'] = analyze_bad_exponents(fooof_params_dict_EyesOpen_2_post,bad_channel_prop=BAD_CHAN_THRESHOLD)['nans']
# %%
#count how many subjects were removed for each task
removed_counts = []
for task_dict, task_name in [
(fooof_params_dict_EyesClosed_1_pre, "EyesClosed_1_pre"),
(fooof_params_dict_EyesClosed_1_post, "EyesClosed_1_post"),
(fooof_params_dict_EyesClosed_2_pre, "EyesClosed_2_pre"),
(fooof_params_dict_EyesClosed_2_post, "EyesClosed_2_post"),
(fooof_params_dict_EyesOpen_1_pre, "EyesOpen_1_pre"),
(fooof_params_dict_EyesOpen_1_post, "EyesOpen_1_post"),
(fooof_params_dict_EyesOpen_2_pre, "EyesOpen_2_pre"),
(fooof_params_dict_EyesOpen_2_post, "EyesOpen_2_post")
]:
removed = sum(not v for v in task_dict['good_mask'].values())
removed_counts.append({'Task': task_name, 'Type': 'FOOOF_Fits', 'Subjects_Removed': removed})
print(f"{task_name}: {removed} subjects removed")
print()
#And of those, howmany removed due to exponents <0
for task_dict, task_name in [
(fooof_params_dict_EyesClosed_1_pre, "EyesClosed_1_pre"),
(fooof_params_dict_EyesClosed_1_post, "EyesClosed_1_post"),
(fooof_params_dict_EyesClosed_2_pre, "EyesClosed_2_pre"),
(fooof_params_dict_EyesClosed_2_post, "EyesClosed_2_post"),
(fooof_params_dict_EyesOpen_1_pre, "EyesOpen_1_pre"),
(fooof_params_dict_EyesOpen_1_post, "EyesOpen_1_post"),
(fooof_params_dict_EyesOpen_2_pre, "EyesOpen_2_pre"),
(fooof_params_dict_EyesOpen_2_post, "EyesOpen_2_post")
]:
removed = sum(not v for v in task_dict['good_mask_exponent'].values())
removed_counts.append({'Task': task_name, 'Type': 'Bad_Exponents', 'Subjects_Removed': removed})
print(f"{task_name} (exponent): {removed} subjects removed")
print()
#Count how many subjects were removed for each task based on missing alpha peaks
for task_dict, task_name in [
(fooof_params_dict_EyesClosed_1_pre, "EyesClosed_1_pre"),
(fooof_params_dict_EyesClosed_1_post, "EyesClosed_1_post"),
(fooof_params_dict_EyesClosed_2_pre, "EyesClosed_2_pre"),
(fooof_params_dict_EyesClosed_2_post, "EyesClosed_2_post"),
(fooof_params_dict_EyesOpen_1_pre, "EyesOpen_1_pre"),
(fooof_params_dict_EyesOpen_1_post, "EyesOpen_1_post"),
(fooof_params_dict_EyesOpen_2_pre, "EyesOpen_2_pre"),
(fooof_params_dict_EyesOpen_2_post, "EyesOpen_2_post")
]:
removed = sum(not v for v in task_dict['good_mask_alpha'].values())
removed_counts.append({'Task': task_name, 'Type': 'Alpha_Peaks', 'Subjects_Removed': removed})
print(f"{task_name} (alpha): {removed} subjects removed")
# Convert to DataFrame and save
removed_counts_df = pd.DataFrame(removed_counts)
removed_counts_df.to_csv(RESULTS_PATH / 'subjects_removed_summary.csv', index=False)
print(f"\nSummary saved to {RESULTS_PATH / 'subjects_removed_summary.csv'}")
# %%
#append the good_mask information to the demographic data frame.
for key, value in fooof_params_dict_EyesClosed_1_pre['good_mask'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Fits_EyesClosed_1_pre'] = value
for key, value in fooof_params_dict_EyesClosed_1_post['good_mask'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Fits_EyesClosed_1_post'] = value
for key, value in fooof_params_dict_EyesClosed_2_pre['good_mask'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Fits_EyesClosed_2_pre'] = value
for key, value in fooof_params_dict_EyesClosed_2_post['good_mask'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Fits_EyesClosed_2_post'] = value
for key, value in fooof_params_dict_EyesOpen_1_pre['good_mask'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Fits_EyesOpen_1_pre'] = value
for key, value in fooof_params_dict_EyesOpen_1_post['good_mask'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Fits_EyesOpen_1_post'] = value
for key, value in fooof_params_dict_EyesOpen_2_pre['good_mask'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Fits_EyesOpen_2_pre'] = value
for key, value in fooof_params_dict_EyesOpen_2_post['good_mask'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Fits_EyesOpen_2_post'] = value
for key, value in fooof_params_dict_EyesClosed_1_pre['good_mask_alpha'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Alpha_EyesClosed_1_pre'] = value
for key, value in fooof_params_dict_EyesClosed_1_post['good_mask_alpha'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Alpha_EyesClosed_1_post'] = value
for key, value in fooof_params_dict_EyesClosed_2_pre['good_mask_alpha'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Alpha_EyesClosed_2_pre'] = value
for key, value in fooof_params_dict_EyesClosed_2_post['good_mask_alpha'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Alpha_EyesClosed_2_post'] = value
for key, value in fooof_params_dict_EyesOpen_1_pre['good_mask_alpha'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Alpha_EyesOpen_1_pre'] = value
for key, value in fooof_params_dict_EyesOpen_1_post['good_mask_alpha'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Alpha_EyesOpen_1_post'] = value
for key, value in fooof_params_dict_EyesOpen_2_pre['good_mask_alpha'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Alpha_EyesOpen_2_pre'] = value
for key, value in fooof_params_dict_EyesOpen_2_post['good_mask_alpha'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Alpha_EyesOpen_2_post'] = value
for key, value in fooof_params_dict_EyesClosed_1_pre['good_mask_exponent'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Exponent_EyesClosed_1_pre'] = value
for key, value in fooof_params_dict_EyesClosed_1_post['good_mask_exponent'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Exponent_EyesClosed_1_post'] = value
for key, value in fooof_params_dict_EyesClosed_2_pre['good_mask_exponent'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Exponent_EyesClosed_2_pre'] = value
for key, value in fooof_params_dict_EyesClosed_2_post['good_mask_exponent'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Exponent_EyesClosed_2_post'] = value
for key, value in fooof_params_dict_EyesOpen_1_pre['good_mask_exponent'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Exponent_EyesOpen_1_pre'] = value
for key, value in fooof_params_dict_EyesOpen_1_post['good_mask_exponent'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Exponent_EyesOpen_1_post'] = value
for key, value in fooof_params_dict_EyesOpen_2_pre['good_mask_exponent'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Exponent_EyesOpen_2_pre'] = value
for key, value in fooof_params_dict_EyesOpen_2_post['good_mask_exponent'].items():
demographic_data.loc[demographic_data['participant_id'] == key, 'Good_Exponent_EyesOpen_2_post'] = value
#calculate sum of good fits across all tasks for each subject and each task type
demographic_data['Total_Good_Fits'] = demographic_data[['Good_Fits_EyesClosed_1_pre', 'Good_Fits_EyesClosed_1_post', 'Good_Fits_EyesClosed_2_pre', 'Good_Fits_EyesClosed_2_post', 'Good_Fits_EyesOpen_1_pre', 'Good_Fits_EyesOpen_1_post', 'Good_Fits_EyesOpen_2_pre', 'Good_Fits_EyesOpen_2_post']].sum(axis=1)
demographic_data['Total_Good_Fits_Pre'] = demographic_data[['Good_Fits_EyesClosed_1_pre', 'Good_Fits_EyesClosed_2_pre', 'Good_Fits_EyesOpen_1_pre', 'Good_Fits_EyesOpen_2_pre']].sum(axis=1)
demographic_data['Total_Good_Fits_Post'] = demographic_data[['Good_Fits_EyesClosed_1_post', 'Good_Fits_EyesClosed_2_post', 'Good_Fits_EyesOpen_1_post', 'Good_Fits_EyesOpen_2_post']].sum(axis=1)
#and for alpha metrics
demographic_data['Total_Good_Alpha'] = demographic_data[['Good_Alpha_EyesClosed_1_pre', 'Good_Alpha_EyesClosed_1_post', 'Good_Alpha_EyesClosed_2_pre', 'Good_Alpha_EyesClosed_2_post', 'Good_Alpha_EyesOpen_1_pre', 'Good_Alpha_EyesOpen_1_post',
'Good_Alpha_EyesOpen_2_pre', 'Good_Alpha_EyesOpen_2_post']].sum(axis=1)
demographic_data['Total_Good_Alpha_Pre'] = demographic_data[['Good_Alpha_EyesClosed_1_pre', 'Good_Alpha_EyesClosed_2_pre', 'Good_Alpha_EyesOpen_1_pre', 'Good_Alpha_EyesOpen_2_pre']].sum(axis=1)
demographic_data['Total_Good_Alpha_Post'] = demographic_data[['Good_Alpha_EyesClosed_1_post', 'Good_Alpha_EyesClosed_2_post', 'Good_Alpha_EyesOpen_1_post', 'Good_Alpha_EyesOpen_2_post']].sum(axis=1)
#and for exponent metrics
demographic_data['Total_Good_Exponent'] = demographic_data[['Good_Exponent_EyesClosed_1_pre', 'Good_Exponent_EyesClosed_1_post', 'Good_Exponent_EyesClosed_2_pre', 'Good_Exponent_EyesClosed_2_post', 'Good_Exponent_EyesOpen_1_pre', 'Good_Exponent_EyesOpen_1_post',
'Good_Exponent_EyesOpen_2_pre', 'Good_Exponent_EyesOpen_2_post']].sum(axis=1)
demographic_data['Total_Good_Exponent_Pre'] = demographic_data[['Good_Exponent_EyesClosed_1_pre', 'Good_Exponent_EyesClosed_2_pre', 'Good_Exponent_EyesOpen_1_pre', 'Good_Exponent_EyesOpen_2_pre']].sum(axis=1)
demographic_data['Total_Good_Exponent_Post'] = demographic_data[['Good_Exponent_EyesClosed_1_post', 'Good_Exponent_EyesClosed_2_post', 'Good_Exponent_EyesOpen_1_post', 'Good_Exponent_EyesOpen_2_post']].sum(axis=1)
# %%
#Here, we create interpolated mean values for each FOOOF parameter in each task dictionary
from utils import replace_nans_with_mean_alpha_new, replace_nans_with_mean_new, replace_nans_with_mean_exponent
fooof_params_dict_EyesClosed_1_pre['exponent_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_1_pre, 'exponents')
fooof_params_dict_EyesClosed_1_post['exponent_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_1_post, 'exponents')
fooof_params_dict_EyesClosed_2_pre['exponent_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_2_pre, 'exponents')
fooof_params_dict_EyesClosed_2_post['exponent_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_2_post, 'exponents')
fooof_params_dict_EyesOpen_1_pre['exponent_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_1_pre, 'exponents')
fooof_params_dict_EyesOpen_1_post['exponent_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_1_post, 'exponents')
fooof_params_dict_EyesOpen_2_pre['exponent_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_2_pre, 'exponents')
fooof_params_dict_EyesOpen_2_post['exponent_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_2_post, 'exponents')
fooof_params_dict_EyesClosed_1_pre['offset_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_1_pre, 'offsets')
fooof_params_dict_EyesClosed_1_post['offset_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_1_post, 'offsets')
fooof_params_dict_EyesClosed_2_pre['offset_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_2_pre, 'offsets')
fooof_params_dict_EyesClosed_2_post['offset_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_2_post, 'offsets')
fooof_params_dict_EyesOpen_1_pre['offset_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_1_pre, 'offsets')
fooof_params_dict_EyesOpen_1_post['offset_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_1_post, 'offsets')
fooof_params_dict_EyesOpen_2_pre['offset_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_2_pre, 'offsets')
fooof_params_dict_EyesOpen_2_post['offset_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_2_post, 'offsets')
fooof_params_dict_EyesClosed_1_pre['fits_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_1_pre, 'fits')
fooof_params_dict_EyesClosed_1_post['fits_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_1_post, 'fits')
fooof_params_dict_EyesClosed_2_pre['fits_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_2_pre, 'fits')
fooof_params_dict_EyesClosed_2_post['fits_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_2_post, 'fits')
fooof_params_dict_EyesOpen_1_pre['fits_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_1_pre, 'fits')
fooof_params_dict_EyesOpen_1_post['fits_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_1_post, 'fits')
fooof_params_dict_EyesOpen_2_pre['fits_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_2_pre, 'fits')
fooof_params_dict_EyesOpen_2_post['fits_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_2_post, 'fits')
fooof_params_dict_EyesClosed_1_pre['error_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_1_pre, 'errors')
fooof_params_dict_EyesClosed_1_post['error_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_1_post, 'errors')
fooof_params_dict_EyesClosed_2_pre['error_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_2_pre, 'errors')
fooof_params_dict_EyesClosed_2_post['error_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesClosed_2_post, 'errors')
fooof_params_dict_EyesOpen_1_pre['error_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_1_pre, 'errors')
fooof_params_dict_EyesOpen_1_post['error_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_1_post, 'errors')
fooof_params_dict_EyesOpen_2_pre['error_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_2_pre, 'errors')
fooof_params_dict_EyesOpen_2_post['error_mean_interp'] = replace_nans_with_mean_new(fooof_params_dict_EyesOpen_2_post, 'errors')
fooof_params_dict_EyesClosed_1_pre['alpha_freq_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesClosed_1_pre, 'alpha_freq')
fooof_params_dict_EyesClosed_1_post['alpha_freq_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesClosed_1_post, 'alpha_freq')
fooof_params_dict_EyesClosed_2_pre['alpha_freq_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesClosed_2_pre, 'alpha_freq')
fooof_params_dict_EyesClosed_2_post['alpha_freq_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesClosed_2_post, 'alpha_freq')
fooof_params_dict_EyesOpen_1_pre['alpha_freq_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesOpen_1_pre, 'alpha_freq')
fooof_params_dict_EyesOpen_1_post['alpha_freq_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesOpen_1_post, 'alpha_freq')
fooof_params_dict_EyesOpen_2_pre['alpha_freq_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesOpen_2_pre, 'alpha_freq')
fooof_params_dict_EyesOpen_2_post['alpha_freq_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesOpen_2_post, 'alpha_freq')
fooof_params_dict_EyesClosed_1_pre['alpha_amp_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesClosed_1_pre, 'alpha_amp')
fooof_params_dict_EyesClosed_1_post['alpha_amp_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesClosed_1_post, 'alpha_amp')
fooof_params_dict_EyesClosed_2_pre['alpha_amp_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesClosed_2_pre, 'alpha_amp')
fooof_params_dict_EyesClosed_2_post['alpha_amp_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesClosed_2_post, 'alpha_amp')
fooof_params_dict_EyesOpen_1_pre['alpha_amp_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesOpen_1_pre, 'alpha_amp')
fooof_params_dict_EyesOpen_1_post['alpha_amp_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesOpen_1_post, 'alpha_amp')
fooof_params_dict_EyesOpen_2_pre['alpha_amp_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesOpen_2_pre, 'alpha_amp')
fooof_params_dict_EyesOpen_2_post['alpha_amp_mean_interp'] = replace_nans_with_mean_alpha_new(fooof_params_dict_EyesOpen_2_post, 'alpha_amp')
#and for the exponent, replace chans with bad exponents with the mean value across that subject.
fooof_params_dict_EyesClosed_1_pre['exponent_mean_interp'] = replace_nans_with_mean_exponent(fooof_params_dict_EyesClosed_1_pre, 'exponents')
fooof_params_dict_EyesClosed_1_post['exponent_mean_interp'] = replace_nans_with_mean_exponent(fooof_params_dict_EyesClosed_1_post, 'exponents')
fooof_params_dict_EyesClosed_2_pre['exponent_mean_interp'] = replace_nans_with_mean_exponent(fooof_params_dict_EyesClosed_2_pre, 'exponents')
fooof_params_dict_EyesClosed_2_post['exponent_mean_interp'] = replace_nans_with_mean_exponent(fooof_params_dict_EyesClosed_2_post, 'exponents')
fooof_params_dict_EyesOpen_1_pre['exponent_mean_interp'] = replace_nans_with_mean_exponent(fooof_params_dict_EyesOpen_1_pre, 'exponents')
fooof_params_dict_EyesOpen_1_post['exponent_mean_interp'] = replace_nans_with_mean_exponent(fooof_params_dict_EyesOpen_1_post, 'exponents')
fooof_params_dict_EyesOpen_2_pre['exponent_mean_interp'] = replace_nans_with_mean_exponent(fooof_params_dict_EyesOpen_2_pre, 'exponents')
fooof_params_dict_EyesOpen_2_post['exponent_mean_interp'] = replace_nans_with_mean_exponent(fooof_params_dict_EyesOpen_2_post, 'exponents')
# %%
# Calculate mean fit and error for each subject for each task and session, adding to demographic_data. This time, we're using the interpolated values
for task in tasks:
for session in ['1', '2']:
for acq in ['pre']:
task_key = f'{task}_{session}_{acq}'
fooof_dict_name = f'fooof_params_dict_{task}_{session}_{acq}'
if fooof_dict_name in globals():
fooof_dict = globals()[fooof_dict_name]
col_name = f'Mean_Fits_interp_{task_key}'
demographic_data[col_name] = demographic_data['participant_id'].map(fooof_dict['fits_mean_interp']).apply(
lambda x: np.nanmean([v for v in x if not np.isnan(v)]) if isinstance(x, (list, np.ndarray)) else np.nan
)
col_name_err = f'Mean_Error_interp_{task_key}'
demographic_data[col_name_err] = demographic_data['participant_id'].map(fooof_dict['error_mean_interp']).apply(
lambda x: np.nanmean([v for v in x if not np.isnan(v)]) if isinstance(x, (list, np.ndarray)) else np.nan
)
#and again, but for non-interpolated values
for task in tasks:
for session in ['1', '2']:
for acq in ['pre']:
task_key = f'{task}_{session}_{acq}'
fooof_dict_name = f'fooof_params_dict_{task}_{session}_{acq}'
if fooof_dict_name in globals():
fooof_dict = globals()[fooof_dict_name]
col_name = f'Mean_Fits_{task_key}'
demographic_data[col_name] = demographic_data['participant_id'].map(fooof_dict['fits']).apply(
lambda x: np.nanmean([v for v in x if not np.isnan(v)]) if isinstance(x, (list, np.ndarray)) else np.nan
)
col_name_err = f'Mean_Error_{task_key}'
demographic_data[col_name_err] = demographic_data['participant_id'].map(fooof_dict['errors']).apply(
lambda x: np.nanmean([v for v in x if not np.isnan(v)]) if isinstance(x, (list, np.ndarray)) else np.nan
)
# %% [markdown]
# ## 4. PCA and KMeans clustering for each task and measure.
# For these, we only include subjects with good fits in all 4 tasks that we care about (EyesClosed_1_pre, EyesClosed_2_pre, EyesOpen_1_pre, EyesOpen_2_pre).
#
# Due to multiple time points, we will only do pca and clustering on time 1 data, but separetly for each measure and eyes open vs closed.
#
# To do, comment code, build functions for the reorganizing of the data for pca and kmeans, and build functions for those.
# %%
#Load subject 101 epoched data to get channel names
chan_info = mne.io.read_info(os.path.join(os.path.abspath('..'),'Data','sub-101','ses-1','eeg','sub-101_ses-1_task-EyesClosed_acq-pre_proc-epo_eeg.fif'))
# %%
#Create dataframes for each measure, which are a subset of subjects with good fits/peaksfor the tasks of interest. We won't use these for the actual ICCs/HLMs, but to identify clusters.
from turtle import left
exponent_df_EyesClosed_1_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesClosed_1_pre['exponent_mean_interp'], orient='index', columns=[f'Exponent_EyesClosed_1_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
exponent_df_EyesClosed_1_pre = exponent_df_EyesClosed_1_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesClosed','Good_Ref_EyesClosed']], left_index=True, right_on='participant_id')
exponent_df_EyesClosed_1_pre = exponent_df_EyesClosed_1_pre[(exponent_df_EyesClosed_1_pre['Total_Good_Fits_Pre'] == 4) & (exponent_df_EyesClosed_1_pre['Good_Retention_EyesClosed'] == 1) & (exponent_df_EyesClosed_1_pre['Good_Ref_EyesClosed'] == 1)]
exponent_df_EyesClosed_1_pre = exponent_df_EyesClosed_1_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesClosed','Good_Ref_EyesClosed'])
exponent_df_EyesClosed_2_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesClosed_2_pre['exponent_mean_interp'], orient='index', columns=[f'Exponent_EyesClosed_2_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
exponent_df_EyesClosed_2_pre = exponent_df_EyesClosed_2_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesClosed','Good_Ref_EyesClosed']], left_index=True, right_on='participant_id')
exponent_df_EyesClosed_2_pre = exponent_df_EyesClosed_2_pre[(exponent_df_EyesClosed_2_pre['Total_Good_Fits_Pre'] == 4) & (exponent_df_EyesClosed_2_pre['Good_Retention_EyesClosed'] == 1) & (exponent_df_EyesClosed_2_pre['Good_Ref_EyesClosed'] == 1)]
exponent_df_EyesClosed_2_pre = exponent_df_EyesClosed_2_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesClosed','Good_Ref_EyesClosed'])
exponent_df_EyesOpen_1_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesOpen_1_pre['exponent_mean_interp'], orient='index', columns=[f'Exponent_EyesOpen_1_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
exponent_df_EyesOpen_1_pre = exponent_df_EyesOpen_1_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesOpen','Good_Ref_EyesOpen']], left_index=True, right_on='participant_id')
exponent_df_EyesOpen_1_pre = exponent_df_EyesOpen_1_pre[(exponent_df_EyesOpen_1_pre['Total_Good_Fits_Pre'] == 4) & (exponent_df_EyesOpen_1_pre['Good_Retention_EyesOpen'] == 1) & (exponent_df_EyesOpen_1_pre['Good_Ref_EyesOpen'] == 1)]
exponent_df_EyesOpen_1_pre = exponent_df_EyesOpen_1_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesOpen','Good_Ref_EyesOpen'])
exponent_df_EyesOpen_2_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesOpen_2_pre['exponent_mean_interp'], orient='index', columns=[f'Exponent_EyesOpen_2_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
exponent_df_EyesOpen_2_pre = exponent_df_EyesOpen_2_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesOpen','Good_Ref_EyesOpen']], left_index=True, right_on='participant_id')
exponent_df_EyesOpen_2_pre = exponent_df_EyesOpen_2_pre[(exponent_df_EyesOpen_2_pre['Total_Good_Fits_Pre'] == 4) & (exponent_df_EyesOpen_2_pre['Good_Retention_EyesOpen'] == 1) & (exponent_df_EyesOpen_2_pre['Good_Ref_EyesOpen'] == 1)]
exponent_df_EyesOpen_2_pre = exponent_df_EyesOpen_2_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesOpen','Good_Ref_EyesOpen'])
offset_df_EyesClosed_1_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesClosed_1_pre['offset_mean_interp'], orient='index', columns=[f'Offset_EyesClosed_1_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
offset_df_EyesClosed_1_pre = offset_df_EyesClosed_1_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesClosed','Good_Ref_EyesClosed']], left_index=True, right_on='participant_id')
offset_df_EyesClosed_1_pre = offset_df_EyesClosed_1_pre[(offset_df_EyesClosed_1_pre['Total_Good_Fits_Pre'] == 4) & (offset_df_EyesClosed_1_pre['Good_Retention_EyesClosed'] == 1) & (offset_df_EyesClosed_1_pre['Good_Ref_EyesClosed'] == 1)]
offset_df_EyesClosed_1_pre = offset_df_EyesClosed_1_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesClosed','Good_Ref_EyesClosed'])
offset_df_EyesClosed_2_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesClosed_2_pre['offset_mean_interp'], orient='index', columns=[f'Offset_EyesClosed_2_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
offset_df_EyesClosed_2_pre = offset_df_EyesClosed_2_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesClosed','Good_Ref_EyesClosed']], left_index=True, right_on='participant_id')
offset_df_EyesClosed_2_pre = offset_df_EyesClosed_2_pre[(offset_df_EyesClosed_2_pre['Total_Good_Fits_Pre'] == 4) & (offset_df_EyesClosed_2_pre['Good_Retention_EyesClosed'] == 1) & (offset_df_EyesClosed_2_pre['Good_Ref_EyesClosed'] == 1)]
offset_df_EyesClosed_2_pre = offset_df_EyesClosed_2_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesClosed','Good_Ref_EyesClosed'])
offset_df_EyesOpen_1_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesOpen_1_pre['offset_mean_interp'], orient='index', columns=[f'Offset_EyesOpen_1_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
offset_df_EyesOpen_1_pre = offset_df_EyesOpen_1_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesOpen','Good_Ref_EyesOpen']], left_index=True, right_on='participant_id')
offset_df_EyesOpen_1_pre = offset_df_EyesOpen_1_pre[(offset_df_EyesOpen_1_pre['Total_Good_Fits_Pre'] == 4) & (offset_df_EyesOpen_1_pre['Good_Retention_EyesOpen'] == 1) & (offset_df_EyesOpen_1_pre['Good_Ref_EyesOpen'] == 1)]
offset_df_EyesOpen_1_pre = offset_df_EyesOpen_1_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesOpen','Good_Ref_EyesOpen'])
offset_df_EyesOpen_2_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesOpen_2_pre['offset_mean_interp'], orient='index', columns=[f'Offset_EyesOpen_2_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
offset_df_EyesOpen_2_pre = offset_df_EyesOpen_2_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesOpen','Good_Ref_EyesOpen']], left_index=True, right_on='participant_id')
offset_df_EyesOpen_2_pre = offset_df_EyesOpen_2_pre[(offset_df_EyesOpen_2_pre['Total_Good_Fits_Pre'] == 4) & (offset_df_EyesOpen_2_pre['Good_Retention_EyesOpen'] == 1) & (offset_df_EyesOpen_2_pre['Good_Ref_EyesOpen'] == 1)]
offset_df_EyesOpen_2_pre = offset_df_EyesOpen_2_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre','Good_Retention_EyesOpen','Good_Ref_EyesOpen'])
alpha_freq_df_EyesClosed_1_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesClosed_1_pre['alpha_freq_mean_interp'], orient='index', columns=[f'AlphaFreq_EyesClosed_1_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
alpha_freq_df_EyesClosed_1_pre = alpha_freq_df_EyesClosed_1_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre','Good_Retention_EyesClosed','Good_Ref_EyesClosed']], left_index=True, right_on='participant_id')
alpha_freq_df_EyesClosed_1_pre = alpha_freq_df_EyesClosed_1_pre[(alpha_freq_df_EyesClosed_1_pre['Total_Good_Fits_Pre'] == 4) & (alpha_freq_df_EyesClosed_1_pre['Total_Good_Alpha_Pre'] == 4) & (alpha_freq_df_EyesClosed_1_pre['Good_Retention_EyesClosed'] == 1) & (alpha_freq_df_EyesClosed_1_pre['Good_Ref_EyesClosed'] == 1)]
alpha_freq_df_EyesClosed_1_pre = alpha_freq_df_EyesClosed_1_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre', 'Good_Retention_EyesClosed','Good_Ref_EyesClosed'])
alpha_freq_df_EyesClosed_2_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesClosed_2_pre['alpha_freq_mean_interp'], orient='index', columns=[f'AlphaFreq_EyesClosed_2_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
alpha_freq_df_EyesClosed_2_pre = alpha_freq_df_EyesClosed_2_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre','Good_Retention_EyesClosed','Good_Ref_EyesClosed']], left_index=True, right_on='participant_id')
alpha_freq_df_EyesClosed_2_pre = alpha_freq_df_EyesClosed_2_pre[(alpha_freq_df_EyesClosed_2_pre['Total_Good_Fits_Pre'] == 4) & (alpha_freq_df_EyesClosed_2_pre['Total_Good_Alpha_Pre'] == 4) & (alpha_freq_df_EyesClosed_2_pre['Good_Retention_EyesClosed'] == 1) & (alpha_freq_df_EyesClosed_2_pre['Good_Ref_EyesClosed'] == 1)]
alpha_freq_df_EyesClosed_2_pre = alpha_freq_df_EyesClosed_2_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre', 'Good_Retention_EyesClosed','Good_Ref_EyesClosed'])
alpha_freq_df_EyesOpen_1_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesOpen_1_pre['alpha_freq_mean_interp'], orient='index', columns=[f'AlphaFreq_EyesOpen_1_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
alpha_freq_df_EyesOpen_1_pre = alpha_freq_df_EyesOpen_1_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre','Good_Retention_EyesOpen','Good_Ref_EyesOpen']], left_index=True, right_on='participant_id')
alpha_freq_df_EyesOpen_1_pre = alpha_freq_df_EyesOpen_1_pre[(alpha_freq_df_EyesOpen_1_pre['Total_Good_Fits_Pre'] == 4) & (alpha_freq_df_EyesOpen_1_pre['Total_Good_Alpha_Pre'] == 4) & (alpha_freq_df_EyesOpen_1_pre['Good_Retention_EyesOpen'] == 1) & (alpha_freq_df_EyesOpen_1_pre['Good_Ref_EyesOpen'] == 1)]
alpha_freq_df_EyesOpen_1_pre = alpha_freq_df_EyesOpen_1_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre', 'Good_Retention_EyesOpen','Good_Ref_EyesOpen'])
alpha_freq_df_EyesOpen_2_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesOpen_2_pre['alpha_freq_mean_interp'], orient='index', columns=[f'AlphaFreq_EyesOpen_2_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
alpha_freq_df_EyesOpen_2_pre = alpha_freq_df_EyesOpen_2_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre','Good_Retention_EyesOpen','Good_Ref_EyesOpen']], left_index=True, right_on='participant_id')
alpha_freq_df_EyesOpen_2_pre = alpha_freq_df_EyesOpen_2_pre[(alpha_freq_df_EyesOpen_2_pre['Total_Good_Fits_Pre'] == 4) & (alpha_freq_df_EyesOpen_2_pre['Total_Good_Alpha_Pre'] == 4) & (alpha_freq_df_EyesOpen_2_pre['Good_Retention_EyesOpen'] == 1) & (alpha_freq_df_EyesOpen_2_pre['Good_Ref_EyesOpen'] == 1)]
alpha_freq_df_EyesOpen_2_pre = alpha_freq_df_EyesOpen_2_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre', 'Good_Retention_EyesOpen','Good_Ref_EyesOpen'])
alpha_amp_df_EyesClosed_1_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesClosed_1_pre['alpha_amp_mean_interp'], orient='index', columns=[f'AlphaAmp_EyesClosed_1_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
alpha_amp_df_EyesClosed_1_pre = alpha_amp_df_EyesClosed_1_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre','Good_Retention_EyesClosed','Good_Ref_EyesClosed']], left_index=True, right_on='participant_id')
alpha_amp_df_EyesClosed_1_pre = alpha_amp_df_EyesClosed_1_pre[(alpha_amp_df_EyesClosed_1_pre['Total_Good_Fits_Pre'] == 4) & (alpha_amp_df_EyesClosed_1_pre['Total_Good_Alpha_Pre'] == 4) & (alpha_amp_df_EyesClosed_1_pre['Good_Retention_EyesClosed'] == 1) & (alpha_amp_df_EyesClosed_1_pre['Good_Ref_EyesClosed'] == 1)]
alpha_amp_df_EyesClosed_1_pre = alpha_amp_df_EyesClosed_1_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre', 'Good_Retention_EyesClosed','Good_Ref_EyesClosed'])
alpha_amp_df_EyesClosed_2_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesClosed_2_pre['alpha_amp_mean_interp'], orient='index', columns=[f'AlphaAmp_EyesClosed_2_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
alpha_amp_df_EyesClosed_2_pre = alpha_amp_df_EyesClosed_2_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre','Good_Retention_EyesClosed','Good_Ref_EyesClosed' ]], left_index=True, right_on='participant_id')
alpha_amp_df_EyesClosed_2_pre = alpha_amp_df_EyesClosed_2_pre[(alpha_amp_df_EyesClosed_2_pre['Total_Good_Fits_Pre'] == 4) & (alpha_amp_df_EyesClosed_2_pre['Total_Good_Alpha_Pre'] == 4) & (alpha_amp_df_EyesClosed_2_pre['Good_Retention_EyesClosed'] == 1) & (alpha_amp_df_EyesClosed_2_pre['Good_Ref_EyesClosed'] == 1)]
alpha_amp_df_EyesClosed_2_pre = alpha_amp_df_EyesClosed_2_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre', 'Good_Retention_EyesClosed','Good_Ref_EyesClosed'])
alpha_amp_df_EyesOpen_1_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesOpen_1_pre['alpha_amp_mean_interp'], orient='index', columns=[f'AlphaAmp_EyesOpen_1_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
alpha_amp_df_EyesOpen_1_pre = alpha_amp_df_EyesOpen_1_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre','Good_Retention_EyesOpen','Good_Ref_EyesOpen']], left_index=True, right_on='participant_id')
alpha_amp_df_EyesOpen_1_pre = alpha_amp_df_EyesOpen_1_pre[(alpha_amp_df_EyesOpen_1_pre['Total_Good_Fits_Pre'] == 4) & (alpha_amp_df_EyesOpen_1_pre['Total_Good_Alpha_Pre'] == 4) & (alpha_amp_df_EyesOpen_1_pre['Good_Retention_EyesOpen'] == 1) & (alpha_amp_df_EyesOpen_1_pre['Good_Ref_EyesOpen'] == 1)]
alpha_amp_df_EyesOpen_1_pre = alpha_amp_df_EyesOpen_1_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre', 'Good_Retention_EyesOpen','Good_Ref_EyesOpen'])
alpha_amp_df_EyesOpen_2_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesOpen_2_pre['alpha_amp_mean_interp'], orient='index', columns=[f'AlphaAmp_EyesOpen_2_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
alpha_amp_df_EyesOpen_2_pre = alpha_amp_df_EyesOpen_2_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre','Good_Retention_EyesOpen','Good_Ref_EyesOpen']], left_index=True, right_on='participant_id')
alpha_amp_df_EyesOpen_2_pre = alpha_amp_df_EyesOpen_2_pre[(alpha_amp_df_EyesOpen_2_pre['Total_Good_Fits_Pre'] == 4) & (alpha_amp_df_EyesOpen_2_pre['Total_Good_Alpha_Pre'] == 4) & (alpha_amp_df_EyesOpen_2_pre['Good_Retention_EyesOpen'] == 1) & (alpha_amp_df_EyesOpen_2_pre['Good_Ref_EyesOpen'] == 1)]
alpha_amp_df_EyesOpen_2_pre = alpha_amp_df_EyesOpen_2_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre', 'Good_Retention_EyesOpen','Good_Ref_EyesOpen'])
#Absolute alpha power *could* use the entire sample of participants, as FOOOF paramaters don't really factor in. However, doing so would make it less comparable to the other alpha metrics.
#Therefore, the same subset of participants with good fits/peaks for the alpha metrics will be used for absolute alpha power, to keep it consistent across measures.
abs_alpha_df_EyesClosed_1_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesClosed_1_pre['alpha_absolute'], orient='index', columns=[f'AbsAlpha_EyesClosed_1_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
abs_alpha_df_EyesClosed_1_pre = abs_alpha_df_EyesClosed_1_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre','Good_Retention_EyesClosed','Good_Ref_EyesClosed']], left_index=True, right_on='participant_id')
abs_alpha_df_EyesClosed_1_pre = abs_alpha_df_EyesClosed_1_pre[(abs_alpha_df_EyesClosed_1_pre['Total_Good_Fits_Pre'] == 4) & (abs_alpha_df_EyesClosed_1_pre['Total_Good_Alpha_Pre'] == 4) & (abs_alpha_df_EyesClosed_1_pre['Good_Retention_EyesClosed'] == 1) & (abs_alpha_df_EyesClosed_1_pre['Good_Ref_EyesClosed'] == 1)]
abs_alpha_df_EyesClosed_1_pre = abs_alpha_df_EyesClosed_1_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre', 'Good_Retention_EyesClosed','Good_Ref_EyesClosed'])
abs_alpha_df_EyesClosed_2_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesClosed_2_pre['alpha_absolute'], orient='index', columns=[f'AbsAlpha_EyesClosed_2_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
abs_alpha_df_EyesClosed_2_pre = abs_alpha_df_EyesClosed_2_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre','Good_Retention_EyesClosed','Good_Ref_EyesClosed']], left_index=True, right_on='participant_id')
abs_alpha_df_EyesClosed_2_pre = abs_alpha_df_EyesClosed_2_pre[(abs_alpha_df_EyesClosed_2_pre['Total_Good_Fits_Pre'] == 4) & (abs_alpha_df_EyesClosed_2_pre['Total_Good_Alpha_Pre'] == 4) & (abs_alpha_df_EyesClosed_2_pre['Good_Retention_EyesClosed'] == 1) & (abs_alpha_df_EyesClosed_2_pre['Good_Ref_EyesClosed'] == 1)]
abs_alpha_df_EyesClosed_2_pre = abs_alpha_df_EyesClosed_2_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre', 'Good_Retention_EyesClosed','Good_Ref_EyesClosed'])
abs_alpha_df_EyesOpen_1_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesOpen_1_pre['alpha_absolute'], orient='index', columns=[f'AbsAlpha_EyesOpen_1_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
abs_alpha_df_EyesOpen_1_pre = abs_alpha_df_EyesOpen_1_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre','Good_Retention_EyesOpen','Good_Ref_EyesOpen']], left_index=True, right_on='participant_id')
abs_alpha_df_EyesOpen_1_pre = abs_alpha_df_EyesOpen_1_pre[(abs_alpha_df_EyesOpen_1_pre['Total_Good_Fits_Pre'] == 4) & (abs_alpha_df_EyesOpen_1_pre['Total_Good_Alpha_Pre'] == 4) & (abs_alpha_df_EyesOpen_1_pre['Good_Retention_EyesOpen'] == 1) & (abs_alpha_df_EyesOpen_1_pre['Good_Ref_EyesOpen'] == 1)]
abs_alpha_df_EyesOpen_1_pre = abs_alpha_df_EyesOpen_1_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre', 'Good_Retention_EyesOpen','Good_Ref_EyesOpen'])
abs_alpha_df_EyesOpen_2_pre = pd.DataFrame.from_dict(fooof_params_dict_EyesOpen_2_pre['alpha_absolute'], orient='index', columns=[f'AbsAlpha_EyesOpen_2_pre_Ch{ch+1}' for ch in range(len(chan_info['ch_names']))])
abs_alpha_df_EyesOpen_2_pre = abs_alpha_df_EyesOpen_2_pre.merge(demographic_data[['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre','Good_Retention_EyesOpen','Good_Ref_EyesOpen']], left_index=True, right_on='participant_id')
abs_alpha_df_EyesOpen_2_pre = abs_alpha_df_EyesOpen_2_pre[(abs_alpha_df_EyesOpen_2_pre['Total_Good_Fits_Pre'] == 4) & (abs_alpha_df_EyesOpen_2_pre['Total_Good_Alpha_Pre'] == 4) & (abs_alpha_df_EyesOpen_2_pre['Good_Retention_EyesOpen'] == 1) & (abs_alpha_df_EyesOpen_2_pre['Good_Ref_EyesOpen'] == 1)]
abs_alpha_df_EyesOpen_2_pre = abs_alpha_df_EyesOpen_2_pre.drop(columns=['participant_id', 'Total_Good_Fits_Pre', 'Total_Good_Alpha_Pre', 'Good_Retention_EyesOpen','Good_Ref_EyesOpen'])
# %%
#Run PCA on each dataframe and general plots.
#Looking at the scree plots, it looks like 2 components is a good choice for all measures and both tasks.
from utils import run_pca_and_plot
plot_pca = False
pca_exp_EyesClosed_1_pre = run_pca_and_plot(exponent_df_EyesClosed_1_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='Exponent EC1 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_exp_EyesClosed_2_pre = run_pca_and_plot(exponent_df_EyesClosed_2_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='Exponent EC2 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_exp_EyesOpen_1_pre = run_pca_and_plot(exponent_df_EyesOpen_1_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='Exponent EO1 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_exp_EyesOpen_2_pre = run_pca_and_plot(exponent_df_EyesOpen_2_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='Exponent EO2 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_off_EyesClosed_1_pre = run_pca_and_plot(offset_df_EyesClosed_1_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='Offset EC1 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_off_EyesClosed_2_pre = run_pca_and_plot(offset_df_EyesClosed_2_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='Offset EC2 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_off_EyesOpen_1_pre = run_pca_and_plot(offset_df_EyesOpen_1_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='Offset EO1 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_off_EyesOpen_2_pre = run_pca_and_plot(offset_df_EyesOpen_2_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='Offset EO2 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_af_EyesClosed_1_pre = run_pca_and_plot(alpha_freq_df_EyesClosed_1_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='AlphaFreq EC1 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_af_EyesClosed_2_pre = run_pca_and_plot(alpha_freq_df_EyesClosed_2_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='AlphaFreq EC2 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_af_EyesOpen_1_pre = run_pca_and_plot(alpha_freq_df_EyesOpen_1_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='AlphaFreq EO1 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_af_EyesOpen_2_pre = run_pca_and_plot(alpha_freq_df_EyesOpen_2_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='AlphaFreq EO2 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_aa_EyesClosed_1_pre = run_pca_and_plot(alpha_amp_df_EyesClosed_1_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='AlphaAmp EC1 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_aa_EyesClosed_2_pre = run_pca_and_plot(alpha_amp_df_EyesClosed_2_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='AlphaAmp EC2 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_aa_EyesOpen_1_pre = run_pca_and_plot(alpha_amp_df_EyesOpen_1_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='AlphaAmp EO1 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_aa_EyesOpen_2_pre = run_pca_and_plot(alpha_amp_df_EyesOpen_2_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='AlphaAmp EO2 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_aaf_EyesClosed_1_pre = run_pca_and_plot(abs_alpha_df_EyesClosed_1_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='AbsAlpha EC1 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_aaf_EyesClosed_2_pre = run_pca_and_plot(abs_alpha_df_EyesClosed_2_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='AbsAlpha EC2 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_aaf_EyesOpen_1_pre = run_pca_and_plot(abs_alpha_df_EyesOpen_1_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='AbsAlpha EO1 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
pca_aaf_EyesOpen_2_pre = run_pca_and_plot(abs_alpha_df_EyesOpen_2_pre, chan_names=chan_info['ch_names'], n_components=10, title_prefix='AbsAlpha EO2 Pre',plot=plot_pca,save_path=str(FIGURES_PATH))
# %%
#Get kmeans diagnostics to help choose k for clustering channels
from utils import kmeans_diagnostics
k_range = range(1,11)
res_diag_exponent_EyesClosed_1_pre = kmeans_diagnostics(exponent_df_EyesClosed_1_pre, k_values=k_range, transpose=True, save_prefix='Exponent_EyesClosed_1_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_exponent_EyesClosed_2_pre = kmeans_diagnostics(exponent_df_EyesClosed_2_pre, k_values=k_range, transpose=True, save_prefix='Exponent_EyesClosed_2_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_exponent_EyesOpen_1_pre = kmeans_diagnostics(exponent_df_EyesOpen_1_pre, k_values=k_range, transpose=True, save_prefix='Exponent_EyesOpen_1_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_exponent_EyesOpen_2_pre = kmeans_diagnostics(exponent_df_EyesOpen_2_pre, k_values=k_range, transpose=True, save_prefix='Exponent_EyesOpen_2_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_offset_EyesClosed_1_pre = kmeans_diagnostics(offset_df_EyesClosed_1_pre, k_values=k_range, transpose=True, save_prefix='Offset_EyesClosed_1_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_offset_EyesClosed_2_pre = kmeans_diagnostics(offset_df_EyesClosed_2_pre, k_values=k_range, transpose=True, save_prefix='Offset_EyesClosed_2_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_offset_EyesOpen_1_pre = kmeans_diagnostics(offset_df_EyesOpen_1_pre, k_values=k_range, transpose=True, save_prefix='Offset_EyesOpen_1_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_offset_EyesOpen_2_pre = kmeans_diagnostics(offset_df_EyesOpen_2_pre, k_values=k_range, transpose=True, save_prefix='Offset_EyesOpen_2_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_alpha_freq_EyesClosed_1_pre = kmeans_diagnostics(alpha_freq_df_EyesClosed_1_pre, k_values=k_range, transpose=True, save_prefix='AlphaFreq_EyesClosed_1_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_alpha_freq_EyesClosed_2_pre = kmeans_diagnostics(alpha_freq_df_EyesClosed_2_pre, k_values=k_range, transpose=True, save_prefix='AlphaFreq_EyesClosed_2_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_alpha_freq_EyesOpen_1_pre = kmeans_diagnostics(alpha_freq_df_EyesOpen_1_pre, k_values=k_range, transpose=True, save_prefix='AlphaFreq_EyesOpen_1_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_alpha_freq_EyesOpen_2_pre = kmeans_diagnostics(alpha_freq_df_EyesOpen_2_pre, k_values=k_range, transpose=True, save_prefix='AlphaFreq_EyesOpen_2_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_alpha_amp_EyesClosed_1_pre = kmeans_diagnostics(alpha_amp_df_EyesClosed_1_pre, k_values=k_range, transpose=True, save_prefix='AlphaAmp_EyesClosed_1_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_alpha_amp_EyesClosed_2_pre = kmeans_diagnostics(alpha_amp_df_EyesClosed_2_pre, k_values=k_range, transpose=True, save_prefix='AlphaAmp_EyesClosed_2_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_alpha_amp_EyesOpen_1_pre = kmeans_diagnostics(alpha_amp_df_EyesOpen_1_pre, k_values=k_range, transpose=True, save_prefix='AlphaAmp_EyesOpen_1_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_alpha_amp_EyesOpen_2_pre = kmeans_diagnostics(alpha_amp_df_EyesOpen_2_pre, k_values=k_range, transpose=True, save_prefix='AlphaAmp_EyesOpen_2_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_abs_alpha_EyesClosed_1_pre = kmeans_diagnostics(abs_alpha_df_EyesClosed_1_pre, k_values=k_range, transpose=True, save_prefix='AbsAlpha_EyesClosed_1_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_abs_alpha_EyesClosed_2_pre = kmeans_diagnostics(abs_alpha_df_EyesClosed_2_pre, k_values=k_range, transpose=True, save_prefix='AbsAlpha_EyesClosed_2_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_abs_alpha_EyesOpen_1_pre = kmeans_diagnostics(abs_alpha_df_EyesOpen_1_pre, k_values=k_range, transpose=True, save_prefix='AbsAlpha_EyesOpen_1_pre', show=False,save_path=str(FIGURES_PATH))
res_diag_abs_alpha_EyesOpen_2_pre = kmeans_diagnostics(abs_alpha_df_EyesOpen_2_pre, k_values=k_range, transpose=True, save_prefix='AbsAlpha_EyesOpen_2_pre', show=False,save_path=str(FIGURES_PATH))
# %% [markdown]
# After clustering from 1 to 11, we can get some diagnostics - we're looking for k's with the highest silhouette scores - bascially, its 2 clusters that is best for everything. The exception is alpha amp eyes open at time 2, but since we're using time 1 for clustering, we can ignore that. As we will see in the HLMs, it there isn't some unusal pattern to the clustering in terms of age or time effects.
# %%
# Print silhouette scores for k=2 to 10 for each measure and task.
print('Silhouette scores for k=2 to 10 for each measure and task:')
print('Exponent_EyesClosed_1_pre:', np.round(res_diag_exponent_EyesClosed_1_pre['silhouette'][1:9], 3))
print('Exponent_EyesClosed_2_pre:', np.round(res_diag_exponent_EyesClosed_2_pre['silhouette'][1:9], 3))
print('Exponent_EyesOpen_1_pre:', np.round(res_diag_exponent_EyesOpen_1_pre['silhouette'][1:9], 3))
print('Exponent_EyesOpen_2_pre:', np.round(res_diag_exponent_EyesOpen_2_pre['silhouette'][1:9], 3))
print('Offset_EyesClosed_1_pre:', np.round(res_diag_offset_EyesClosed_1_pre['silhouette'][1:9], 3))
print('Offset_EyesClosed_2_pre:', np.round(res_diag_offset_EyesClosed_2_pre['silhouette'][1:9], 3))
print('Offset_EyesOpen_1_pre:', np.round(res_diag_offset_EyesOpen_1_pre['silhouette'][1:9], 3))
print('Offset_EyesOpen_2_pre:', np.round(res_diag_offset_EyesOpen_2_pre['silhouette'][1:9], 3))
print('AlphaFreq_EyesClosed_1_pre:', np.round(res_diag_alpha_freq_EyesClosed_1_pre['silhouette'][1:9], 3))
print('AlphaFreq_EyesClosed_2_pre:', np.round(res_diag_alpha_freq_EyesClosed_2_pre['silhouette'][1:9], 3))
print('AlphaFreq_EyesOpen_1_pre:', np.round(res_diag_alpha_freq_EyesOpen_1_pre['silhouette'][1:9], 3))
print('AlphaFreq_EyesOpen_2_pre:', np.round(res_diag_alpha_freq_EyesOpen_2_pre['silhouette'][1:9], 3))
print('AlphaAmp_EyesClosed_1_pre:', np.round(res_diag_alpha_amp_EyesClosed_1_pre['silhouette'][1:9], 3))
print('AlphaAmp_EyesClosed_2_pre:', np.round(res_diag_alpha_amp_EyesClosed_2_pre['silhouette'][1:9], 3))
print('AlphaAmp_EyesOpen_1_pre:', np.round(res_diag_alpha_amp_EyesOpen_1_pre['silhouette'][1:9], 3))
print('AlphaAmp_EyesOpen_2_pre:', np.round(res_diag_alpha_amp_EyesOpen_2_pre['silhouette'][1:9], 3))
print('AbsAlpha_EyesClosed_1_pre:', np.round(res_diag_abs_alpha_EyesClosed_1_pre['silhouette'][1:9], 3))
print('AbsAlpha_EyesClosed_2_pre:', np.round(res_diag_abs_alpha_EyesClosed_2_pre['silhouette'][1:9], 3))
print('AbsAlpha_EyesOpen_1_pre:', np.round(res_diag_abs_alpha_EyesOpen_1_pre['silhouette'][1:9], 3))
print('AbsAlpha_EyesOpen_2_pre:', np.round(res_diag_abs_alpha_EyesOpen_2_pre['silhouette'][1:9],3))
diag_dfs = []
for res_diag, name in zip([res_diag_exponent_EyesClosed_1_pre, res_diag_exponent_EyesClosed_2_pre, res_diag_exponent_EyesOpen_1_pre, res_diag_exponent_EyesOpen_2_pre,
res_diag_offset_EyesClosed_1_pre, res_diag_offset_EyesClosed_2_pre, res_diag_offset_EyesOpen_1_pre, res_diag_offset_EyesOpen_2_pre,
res_diag_alpha_freq_EyesClosed_1_pre, res_diag_alpha_freq_EyesClosed_2_pre, res_diag_alpha_freq_EyesOpen_1_pre, res_diag_alpha_freq_EyesOpen_2_pre,
res_diag_alpha_amp_EyesClosed_1_pre, res_diag_alpha_amp_EyesClosed_2_pre, res_diag_alpha_amp_EyesOpen_1_pre, res_diag_alpha_amp_EyesOpen_2_pre,
res_diag_abs_alpha_EyesClosed_1_pre, res_diag_abs_alpha_EyesClosed_2_pre, res_diag_abs_alpha_EyesOpen_1_pre, res_diag_abs_alpha_EyesOpen_2_pre],
['Exponent_EyesClosed_1_pre', 'Exponent_EyesClosed_2_pre', 'Exponent_EyesOpen_1_pre', 'Exponent_EyesOpen_2_pre',
'Offset_EyesClosed_1_pre', 'Offset_EyesClosed_2_pre', 'Offset_EyesOpen_1_pre', 'Offset_EyesOpen_2_pre',
'AlphaFreq_EyesClosed_1_pre', 'AlphaFreq_EyesClosed_2_pre', 'AlphaFreq_EyesOpen_1_pre', 'AlphaFreq_EyesOpen_2_pre',
'AlphaAmp_EyesClosed_1_pre', 'AlphaAmp_EyesClosed_2_pre', 'AlphaAmp_EyesOpen_1_pre', 'AlphaAmp_EyesOpen_2_pre',
'AbsAlpha_EyesClosed_1_pre', 'AbsAlpha_EyesClosed_2_pre', 'AbsAlpha_EyesOpen_1_pre', 'AbsAlpha_EyesOpen_2_pre']):
sil = np.array(res_diag['silhouette'], dtype=float)[1:]
row = pd.DataFrame([np.round(sil, 3)],
columns=[f'k={k}' for k in range(2, 2 + len(sil))])
row.insert(0, 'measure_task', name)
diag_dfs.append(row)
diag_results_df = pd.concat(diag_dfs, ignore_index=True)
diag_results_df.to_csv(RESULTS_PATH / 'kmeans_diagnostics.csv', index=False)
# %%
#Assign channels to k=2 clusters based on kmeans clustering
from utils import kmeans_cluster_and_save
res_km_exponent_EyesClosed_1_pre = kmeans_cluster_and_save(exponent_df_EyesClosed_1_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='Exponent_EyesClosed_1_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_exponent_EyesClosed_2_pre = kmeans_cluster_and_save(exponent_df_EyesClosed_2_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='Exponent_EyesClosed_2_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_exponent_EyesOpen_1_pre = kmeans_cluster_and_save(exponent_df_EyesOpen_1_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='Exponent_EyesOpen_1_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_exponent_EyesOpen_2_pre = kmeans_cluster_and_save(exponent_df_EyesOpen_2_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='Exponent_EyesOpen_2_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_offset_EyesClosed_1_pre = kmeans_cluster_and_save(offset_df_EyesClosed_1_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='Offset_EyesClosed_1_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_offset_EyesClosed_2_pre = kmeans_cluster_and_save(offset_df_EyesClosed_2_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='Offset_EyesClosed_2_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_offset_EyesOpen_1_pre = kmeans_cluster_and_save(offset_df_EyesOpen_1_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='Offset_EyesOpen_1_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_offset_EyesOpen_2_pre = kmeans_cluster_and_save(offset_df_EyesOpen_2_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='Offset_EyesOpen_2_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_alpha_freq_EyesClosed_1_pre = kmeans_cluster_and_save(alpha_freq_df_EyesClosed_1_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='AlphaFreq_EyesClosed_1_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_alpha_freq_EyesClosed_2_pre = kmeans_cluster_and_save(alpha_freq_df_EyesClosed_2_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='AlphaFreq_EyesClosed_2_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_alpha_freq_EyesOpen_1_pre = kmeans_cluster_and_save(alpha_freq_df_EyesOpen_1_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='AlphaFreq_EyesOpen_1_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_alpha_freq_EyesOpen_2_pre = kmeans_cluster_and_save(alpha_freq_df_EyesOpen_2_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='AlphaFreq_EyesOpen_2_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_alpha_amp_EyesClosed_1_pre = kmeans_cluster_and_save(alpha_amp_df_EyesClosed_1_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='AlphaAmp_EyesClosed_1_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_alpha_amp_EyesClosed_2_pre = kmeans_cluster_and_save(alpha_amp_df_EyesClosed_2_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='AlphaAmp_EyesClosed_2_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_alpha_amp_EyesOpen_1_pre = kmeans_cluster_and_save(alpha_amp_df_EyesOpen_1_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='AlphaAmp_EyesOpen_1_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_alpha_amp_EyesOpen_2_pre = kmeans_cluster_and_save(alpha_amp_df_EyesOpen_2_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='AlphaAmp_EyesOpen_2_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_abs_alpha_EyesClosed_1_pre = kmeans_cluster_and_save(abs_alpha_df_EyesClosed_1_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='AbsAlpha_EyesClosed_1_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_abs_alpha_EyesClosed_2_pre = kmeans_cluster_and_save(abs_alpha_df_EyesClosed_2_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='AbsAlpha_EyesClosed_2_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_abs_alpha_EyesOpen_1_pre = kmeans_cluster_and_save(abs_alpha_df_EyesOpen_1_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='AbsAlpha_EyesOpen_1_pre',
save_figs=True, save_path=str(FIGURES_PATH))
res_km_abs_alpha_EyesOpen_2_pre = kmeans_cluster_and_save(abs_alpha_df_EyesOpen_2_pre, k=2, transpose=True,
chan_names=chan_info['ch_names'], save_prefix='AbsAlpha_EyesOpen_2_pre',
save_figs=True, save_path=str(FIGURES_PATH))
# %% [markdown]
# Here we generate some topomaps of the clusters for each measure and task.
# %%
import numpy as np
import matplotlib.pyplot as plt
import mne
from matplotlib.colors import ListedColormap, BoundaryNorm
def plot_dual_cluster_topomap(cluster1_mask, cluster2_mask, title_prefix, chan_info):
assert np.all(cluster1_mask ^ cluster2_mask), "Masks must be disjoint and cover all channels"
data = np.where(cluster2_mask, 1, 0)
cmap = ListedColormap(["#d33ce732", "#2e86de5a"])
fig, ax = plt.subplots(figsize=(6, 5))
im, _ = mne.viz.plot_topomap(
data,
chan_info,
axes=ax,
contours=0,
image_interp='nearest',
cmap=cmap,
show=False,
names=chan_info['ch_names'],
sensors = False,
extrapolate='box',
size=15
)
ax.set_title(f'{title_prefix}')
for text in ax.texts:
text.set_fontsize(6)
plt.tight_layout()
fig.savefig(FIGURES_PATH / f'{title_prefix}_Both_Clusters_Topomap.png', dpi=600)
plt.close()
#Exponent
plot_dual_cluster_topomap(res_km_exponent_EyesClosed_1_pre['labels'] == 0, res_km_exponent_EyesClosed_1_pre['labels'] == 1, 'Exponent EC', chan_info)
plot_dual_cluster_topomap(res_km_exponent_EyesOpen_1_pre['labels'] == 0, res_km_exponent_EyesOpen_1_pre['labels'] == 1, 'Exponent EO', chan_info)
#Offset
plot_dual_cluster_topomap(res_km_offset_EyesClosed_1_pre['labels'] == 0, res_km_offset_EyesClosed_1_pre['labels'] == 1, 'Offset EC', chan_info)
plot_dual_cluster_topomap(res_km_offset_EyesOpen_1_pre['labels'] == 0, res_km_offset_EyesOpen_1_pre['labels'] == 1, 'Offset EO', chan_info)
#Alpha Frequency
plot_dual_cluster_topomap(res_km_alpha_freq_EyesClosed_1_pre['labels'] == 0, res_km_alpha_freq_EyesClosed_1_pre['labels'] == 1, 'IAPF EC', chan_info)
plot_dual_cluster_topomap(res_km_alpha_freq_EyesOpen_1_pre['labels'] == 0, res_km_alpha_freq_EyesOpen_1_pre['labels'] == 1, 'IAPF EO', chan_info)
#Alpha Amplitude
plot_dual_cluster_topomap(res_km_alpha_amp_EyesClosed_1_pre['labels'] == 0, res_km_alpha_amp_EyesClosed_1_pre['labels'] == 1, 'Alpha Power EC', chan_info)
plot_dual_cluster_topomap(res_km_alpha_amp_EyesOpen_1_pre['labels'] == 0, res_km_alpha_amp_EyesOpen_1_pre['labels'] == 1, 'Alpha Power EO', chan_info)
#Absolute Alpha Power
plot_dual_cluster_topomap(res_km_abs_alpha_EyesClosed_1_pre['labels'] == 0, res_km_abs_alpha_EyesClosed_1_pre['labels'] == 1, 'Absolute \u03B1 Power EC', chan_info)
plot_dual_cluster_topomap(res_km_abs_alpha_EyesOpen_1_pre['labels'] == 0, res_km_abs_alpha_EyesOpen_1_pre['labels'] == 1, 'Absolute \u03B1 Power EO', chan_info)
def plot_dual_cluster_topomap_subplot(cluster1_mask, cluster2_mask, title_prefix, chan_info, ax):
assert np.all(cluster1_mask ^ cluster2_mask), "Masks must be disjoint and cover all channels"
data = np.where(cluster2_mask, 1, 0)
cmap = ListedColormap(["#d33ce732", "#2e86de5a"])
im, _ = mne.viz.plot_topomap(
data,
chan_info,
axes=ax,
contours=0,
image_interp='nearest',
cmap=cmap,
show=False,
names=chan_info['ch_names'],
sensors=False,
extrapolate='box',
size=15
)
ax.set_title(title_prefix)
for text in ax.texts:
text.set_fontsize(6)
fig, axs = plt.subplots(2, 4, figsize=(16, 8))
plot_dual_cluster_topomap_subplot(res_km_exponent_EyesClosed_1_pre['labels'] == 0, res_km_exponent_EyesClosed_1_pre['labels'] == 1, 'Exponent EC', chan_info, axs[0,0])
plot_dual_cluster_topomap_subplot(res_km_exponent_EyesOpen_1_pre['labels'] == 0, res_km_exponent_EyesOpen_1_pre['labels'] == 1, 'Exponent EO', chan_info, axs[0,1])
plot_dual_cluster_topomap_subplot(res_km_offset_EyesClosed_1_pre['labels'] == 0, res_km_offset_EyesClosed_1_pre['labels'] == 1, 'Offset EC', chan_info, axs[0,2])
plot_dual_cluster_topomap_subplot(res_km_offset_EyesOpen_1_pre['labels'] == 0, res_km_offset_EyesOpen_1_pre['labels'] == 1, 'Offset EO', chan_info, axs[0,3])
plot_dual_cluster_topomap_subplot(res_km_alpha_freq_EyesClosed_1_pre['labels'] == 0, res_km_alpha_freq_EyesClosed_1_pre['labels'] == 1, 'IAPF EC', chan_info, axs[1,0])
plot_dual_cluster_topomap_subplot(res_km_alpha_freq_EyesOpen_1_pre['labels'] == 0, res_km_alpha_freq_EyesOpen_1_pre['labels'] == 1, 'IAPF EO', chan_info, axs[1,1])
plot_dual_cluster_topomap_subplot(res_km_alpha_amp_EyesClosed_1_pre['labels'] == 0, res_km_alpha_amp_EyesClosed_1_pre['labels'] == 1, ' \u03B1 PW EC', chan_info, axs[1,2])
plot_dual_cluster_topomap_subplot(res_km_alpha_amp_EyesOpen_1_pre['labels'] == 0, res_km_alpha_amp_EyesOpen_1_pre['labels'] == 1, '\u03B1 PW EO', chan_info, axs[1,3])
plt.tight_layout()
plt.savefig(FIGURES_PATH / 'All_Measures_Clusters_Topomap.png', dpi=600)
plt.close()
# %% [markdown]
# Lastly, we create dataframes that put each channel into an appropriate cluster, for each measure and task.
# This is applied to all subjects, because we 1) didn't use bad subjects to define clusters, and 2) we want to use the same clusters for all subjects,
# and 3) we will drop them from the ICC and HLM data anyway
# %%
from utils import cluster_params_to_dfs
#Do time 1 first
exponent_EyesClosed_1_pre_clusters = cluster_params_to_dfs(res_km_exponent_EyesClosed_1_pre, fooof_params_dict_EyesClosed_1_pre['exponent_mean_interp'])
exponent_EyesOpen_1_pre_clusters = cluster_params_to_dfs(res_km_exponent_EyesOpen_1_pre, fooof_params_dict_EyesOpen_1_pre['exponent_mean_interp'])
offset_EyesClosed_1_pre_clusters = cluster_params_to_dfs(res_km_offset_EyesClosed_1_pre, fooof_params_dict_EyesClosed_1_pre['offset_mean_interp'])
offset_EyesOpen_1_pre_clusters = cluster_params_to_dfs(res_km_offset_EyesOpen_1_pre, fooof_params_dict_EyesOpen_1_pre['offset_mean_interp'])
alpha_freq_EyesClosed_1_pre_clusters = cluster_params_to_dfs(res_km_alpha_freq_EyesClosed_1_pre, fooof_params_dict_EyesClosed_1_pre['alpha_freq_mean_interp'])
alpha_freq_EyesOpen_1_pre_clusters = cluster_params_to_dfs(res_km_alpha_freq_EyesOpen_1_pre, fooof_params_dict_EyesOpen_1_pre['alpha_freq_mean_interp'])
alpha_amp_EyesClosed_1_pre_clusters = cluster_params_to_dfs(res_km_alpha_amp_EyesClosed_1_pre, fooof_params_dict_EyesClosed_1_pre['alpha_amp_mean_interp'])
alpha_amp_EyesOpen_1_pre_clusters = cluster_params_to_dfs(res_km_alpha_amp_EyesOpen_1_pre, fooof_params_dict_EyesOpen_1_pre['alpha_amp_mean_interp'])
alpha_abs_amp_EyesClosed_1_pre_clusters = cluster_params_to_dfs(res_km_abs_alpha_EyesClosed_1_pre, fooof_params_dict_EyesClosed_1_pre['alpha_absolute'])
alpha_abs_amp_EyesOpen_1_pre_clusters = cluster_params_to_dfs(res_km_abs_alpha_EyesOpen_1_pre, fooof_params_dict_EyesOpen_1_pre['alpha_absolute'])
#Then time 2, using time 1 cluster info
exponent_EyesClosed_2_pre_clusters = cluster_params_to_dfs(res_km_exponent_EyesClosed_1_pre, fooof_params_dict_EyesClosed_2_pre['exponent_mean_interp'])
exponent_EyesOpen_2_pre_clusters = cluster_params_to_dfs(res_km_exponent_EyesOpen_1_pre, fooof_params_dict_EyesOpen_2_pre['exponent_mean_interp'])
offset_EyesClosed_2_pre_clusters = cluster_params_to_dfs(res_km_offset_EyesClosed_1_pre, fooof_params_dict_EyesClosed_2_pre['offset_mean_interp'])
offset_EyesOpen_2_pre_clusters = cluster_params_to_dfs(res_km_offset_EyesOpen_1_pre, fooof_params_dict_EyesOpen_2_pre['offset_mean_interp'])
alpha_freq_EyesClosed_2_pre_clusters = cluster_params_to_dfs(res_km_alpha_freq_EyesClosed_1_pre, fooof_params_dict_EyesClosed_2_pre['alpha_freq_mean_interp'])
alpha_freq_EyesOpen_2_pre_clusters = cluster_params_to_dfs(res_km_alpha_freq_EyesOpen_1_pre, fooof_params_dict_EyesOpen_2_pre['alpha_freq_mean_interp'])
alpha_amp_EyesClosed_2_pre_clusters = cluster_params_to_dfs(res_km_alpha_amp_EyesClosed_1_pre, fooof_params_dict_EyesClosed_2_pre['alpha_amp_mean_interp'])
alpha_amp_EyesOpen_2_pre_clusters = cluster_params_to_dfs(res_km_alpha_amp_EyesOpen_1_pre, fooof_params_dict_EyesOpen_2_pre['alpha_amp_mean_interp'])
alpha_abs_amp_EyesClosed_2_pre_clusters = cluster_params_to_dfs(res_km_abs_alpha_EyesClosed_2_pre, fooof_params_dict_EyesClosed_2_pre['alpha_absolute'])
alpha_abs_amp_EyesOpen_2_pre_clusters = cluster_params_to_dfs(res_km_abs_alpha_EyesOpen_2_pre, fooof_params_dict_EyesOpen_2_pre['alpha_absolute'])
#calculate row means for each cluster, for each measure, for each task, for each timepoint
exponent_EyesClosed_1_pre_cluster1_means = exponent_EyesClosed_1_pre_clusters[0].mean(axis=1)
exponent_EyesClosed_1_pre_cluster2_means = exponent_EyesClosed_1_pre_clusters[1].mean(axis=1)
exponent_EyesClosed_2_pre_cluster1_means = exponent_EyesClosed_2_pre_clusters[0].mean(axis=1)
exponent_EyesClosed_2_pre_cluster2_means = exponent_EyesClosed_2_pre_clusters[1].mean(axis=1)
exponent_EyesOpen_1_pre_cluster1_means = exponent_EyesOpen_1_pre_clusters[0].mean(axis=1)
exponent_EyesOpen_1_pre_cluster2_means = exponent_EyesOpen_1_pre_clusters[1].mean(axis=1)
exponent_EyesOpen_2_pre_cluster1_means = exponent_EyesOpen_2_pre_clusters[0].mean(axis=1)
exponent_EyesOpen_2_pre_cluster2_means = exponent_EyesOpen_2_pre_clusters[1].mean(axis=1)
offset_EyesClosed_1_pre_cluster1_means = offset_EyesClosed_1_pre_clusters[0].mean(axis=1)
offset_EyesClosed_1_pre_cluster2_means = offset_EyesClosed_1_pre_clusters[1].mean(axis=1)
offset_EyesClosed_2_pre_cluster1_means = offset_EyesClosed_2_pre_clusters[0].mean(axis=1)
offset_EyesClosed_2_pre_cluster2_means = offset_EyesClosed_2_pre_clusters[1].mean(axis=1)
offset_EyesOpen_1_pre_cluster1_means = offset_EyesOpen_1_pre_clusters[0].mean(axis=1)
offset_EyesOpen_1_pre_cluster2_means = offset_EyesOpen_1_pre_clusters[1].mean(axis=1)
offset_EyesOpen_2_pre_cluster1_means = offset_EyesOpen_2_pre_clusters[0].mean(axis=1)
offset_EyesOpen_2_pre_cluster2_means = offset_EyesOpen_2_pre_clusters[1].mean(axis=1)
alpha_freq_EyesClosed_1_pre_cluster1_means = alpha_freq_EyesClosed_1_pre_clusters[0].mean(axis=1)
alpha_freq_EyesClosed_1_pre_cluster2_means = alpha_freq_EyesClosed_1_pre_clusters[1].mean(axis=1)
alpha_freq_EyesClosed_2_pre_cluster1_means = alpha_freq_EyesClosed_2_pre_clusters[0].mean(axis=1)
alpha_freq_EyesClosed_2_pre_cluster2_means = alpha_freq_EyesClosed_2_pre_clusters[1].mean(axis=1)
alpha_freq_EyesOpen_1_pre_cluster1_means = alpha_freq_EyesOpen_1_pre_clusters[0].mean(axis=1)
alpha_freq_EyesOpen_1_pre_cluster2_means = alpha_freq_EyesOpen_1_pre_clusters[1].mean(axis=1)
alpha_freq_EyesOpen_2_pre_cluster1_means = alpha_freq_EyesOpen_2_pre_clusters[0].mean(axis=1)
alpha_freq_EyesOpen_2_pre_cluster2_means = alpha_freq_EyesOpen_2_pre_clusters[1].mean(axis=1)
alpha_amp_EyesClosed_1_pre_cluster1_means = alpha_amp_EyesClosed_1_pre_clusters[0].mean(axis=1)
alpha_amp_EyesClosed_1_pre_cluster2_means = alpha_amp_EyesClosed_1_pre_clusters[1].mean(axis=1)
alpha_amp_EyesClosed_2_pre_cluster1_means = alpha_amp_EyesClosed_2_pre_clusters[0].mean(axis=1)
alpha_amp_EyesClosed_2_pre_cluster2_means = alpha_amp_EyesClosed_2_pre_clusters[1].mean(axis=1)
alpha_amp_EyesOpen_1_pre_cluster1_means = alpha_amp_EyesOpen_1_pre_clusters[0].mean(axis=1)
alpha_amp_EyesOpen_1_pre_cluster2_means = alpha_amp_EyesOpen_1_pre_clusters[1].mean(axis=1)
alpha_amp_EyesOpen_2_pre_cluster1_means = alpha_amp_EyesOpen_2_pre_clusters[0].mean(axis=1)
alpha_amp_EyesOpen_2_pre_cluster2_means = alpha_amp_EyesOpen_2_pre_clusters[1].mean(axis=1)
alpha_abs_amp_EyesClosed_1_pre_cluster1_means = alpha_abs_amp_EyesClosed_1_pre_clusters[0].mean(axis=1)
alpha_abs_amp_EyesClosed_1_pre_cluster2_means = alpha_abs_amp_EyesClosed_1_pre_clusters[1].mean(axis=1)
alpha_abs_amp_EyesClosed_2_pre_cluster1_means = alpha_abs_amp_EyesClosed_2_pre_clusters[0].mean(axis=1)
alpha_abs_amp_EyesClosed_2_pre_cluster2_means = alpha_abs_amp_EyesClosed_2_pre_clusters[1].mean(axis=1)
alpha_abs_amp_EyesOpen_1_pre_cluster1_means = alpha_abs_amp_EyesOpen_1_pre_clusters[0].mean(axis=1)
alpha_abs_amp_EyesOpen_1_pre_cluster2_means = alpha_abs_amp_EyesOpen_1_pre_clusters[1].mean(axis=1)
alpha_abs_amp_EyesOpen_2_pre_cluster1_means = alpha_abs_amp_EyesOpen_2_pre_clusters[0].mean(axis=1)
alpha_abs_amp_EyesOpen_2_pre_cluster2_means = alpha_abs_amp_EyesOpen_2_pre_clusters[1].mean(axis=1)
# %%
#Merge everything.
demographic_data_hlm = demographic_data.copy()
demographic_data_hlm.set_index('participant_id', inplace=True)
demographic_data_hlm = demographic_data_hlm.merge(exponent_EyesClosed_1_pre_cluster1_means.rename('Exponent_EyesClosed_1_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(exponent_EyesClosed_1_pre_cluster2_means.rename('Exponent_EyesClosed_1_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(exponent_EyesClosed_2_pre_cluster1_means.rename('Exponent_EyesClosed_2_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(exponent_EyesClosed_2_pre_cluster2_means.rename('Exponent_EyesClosed_2_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(exponent_EyesOpen_1_pre_cluster1_means.rename('Exponent_EyesOpen_1_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(exponent_EyesOpen_1_pre_cluster2_means.rename('Exponent_EyesOpen_1_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(exponent_EyesOpen_2_pre_cluster1_means.rename('Exponent_EyesOpen_2_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(exponent_EyesOpen_2_pre_cluster2_means.rename('Exponent_EyesOpen_2_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(offset_EyesClosed_1_pre_cluster1_means.rename('Offset_EyesClosed_1_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(offset_EyesClosed_1_pre_cluster2_means.rename('Offset_EyesClosed_1_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(offset_EyesClosed_2_pre_cluster1_means.rename('Offset_EyesClosed_2_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(offset_EyesClosed_2_pre_cluster2_means.rename('Offset_EyesClosed_2_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(offset_EyesOpen_1_pre_cluster1_means.rename('Offset_EyesOpen_1_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(offset_EyesOpen_1_pre_cluster2_means.rename('Offset_EyesOpen_1_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(offset_EyesOpen_2_pre_cluster1_means.rename('Offset_EyesOpen_2_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(offset_EyesOpen_2_pre_cluster2_means.rename('Offset_EyesOpen_2_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_freq_EyesClosed_1_pre_cluster1_means.rename('AlphaFreq_EyesClosed_1_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_freq_EyesClosed_1_pre_cluster2_means.rename('AlphaFreq_EyesClosed_1_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_freq_EyesClosed_2_pre_cluster1_means.rename('AlphaFreq_EyesClosed_2_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_freq_EyesClosed_2_pre_cluster2_means.rename('AlphaFreq_EyesClosed_2_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_freq_EyesOpen_1_pre_cluster1_means.rename('AlphaFreq_EyesOpen_1_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_freq_EyesOpen_1_pre_cluster2_means.rename('AlphaFreq_EyesOpen_1_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_freq_EyesOpen_2_pre_cluster1_means.rename('AlphaFreq_EyesOpen_2_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_freq_EyesOpen_2_pre_cluster2_means.rename('AlphaFreq_EyesOpen_2_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_amp_EyesClosed_1_pre_cluster1_means.rename('AlphaAmp_EyesClosed_1_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_amp_EyesClosed_1_pre_cluster2_means.rename('AlphaAmp_EyesClosed_1_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_amp_EyesClosed_2_pre_cluster1_means.rename('AlphaAmp_EyesClosed_2_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_amp_EyesClosed_2_pre_cluster2_means.rename('AlphaAmp_EyesClosed_2_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_amp_EyesOpen_1_pre_cluster1_means.rename('AlphaAmp_EyesOpen_1_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_amp_EyesOpen_1_pre_cluster2_means.rename('AlphaAmp_EyesOpen_1_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_amp_EyesOpen_2_pre_cluster1_means.rename('AlphaAmp_EyesOpen_2_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_amp_EyesOpen_2_pre_cluster2_means.rename('AlphaAmp_EyesOpen_2_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_abs_amp_EyesClosed_1_pre_cluster1_means.rename('AbsAlpha_EyesClosed_1_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_abs_amp_EyesClosed_1_pre_cluster2_means.rename('AbsAlpha_EyesClosed_1_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_abs_amp_EyesClosed_2_pre_cluster1_means.rename('AbsAlpha_EyesClosed_2_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_abs_amp_EyesClosed_2_pre_cluster2_means.rename('AbsAlpha_EyesClosed_2_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_abs_amp_EyesOpen_1_pre_cluster1_means.rename('AbsAlpha_EyesOpen_1_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_abs_amp_EyesOpen_1_pre_cluster2_means.rename('AbsAlpha_EyesOpen_1_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_abs_amp_EyesOpen_2_pre_cluster1_means.rename('AbsAlpha_EyesOpen_2_pre_Cluster1'), left_index=True, right_index=True, how='left')
demographic_data_hlm = demographic_data_hlm.merge(alpha_abs_amp_EyesOpen_2_pre_cluster2_means.rename('AbsAlpha_EyesOpen_2_pre_Cluster2'), left_index=True, right_index=True, how='left')
demographic_data_hlm.to_csv(RESULTS_PATH / 'demographic_data_hlm_pre.csv')
# %% [markdown]
# ## 5. ICCs
# For these, we will first do them at the channel level, because we're not as worried about multiple comparisons as we are about reliability in general, followed by cluster level tests.
# %%
#Set to be ICC(2,1) - two way random effects, absolute agreement, single rater/measurement
icc_type = ["twoway", "agreement", "single"]
icc_chan_ec_exponent, icc_f_ec_exponent, icc_p_ec_exponent,icc_lower_ec_exponent,icc_uppwer_ec_exponent, icc_sample_size_ec_exponent = utils.icc_chans_max(np.array(exponent_df_EyesClosed_1_pre), np.array(exponent_df_EyesClosed_2_pre), icc_type)
icc_chan_eo_exponent, icc_f_eo_exponent, icc_p_eo_exponent, icc_lower_eo_exponent, icc_upper_ec_exponent, icc_sample_size_eo_exponent = utils.icc_chans_max(np.array(exponent_df_EyesOpen_1_pre), np.array(exponent_df_EyesOpen_2_pre), icc_type)
icc_chan_ec_offset, icc_f_ec_offset, icc_p_ec_offset, icc_ec_lower_offset,icc_ec_upper_offset, icc_sample_size_ec_offset = utils.icc_chans_max(np.array(offset_df_EyesClosed_1_pre), np.array(offset_df_EyesClosed_2_pre), icc_type)
icc_chan_eo_offset, icc_f_eo_offset, icc_p_eo_offset, icc_eo_lower_offset, icc_eo_upper_offset,icc_sample_size_eo_offset = utils.icc_chans_max(np.array(offset_df_EyesOpen_1_pre), np.array(offset_df_EyesOpen_2_pre), icc_type)
icc_chan_ec_af, icc_f_ec_af, icc_p_ec_af, icc_lower_ec_af, icc_upper_ec_af, icc_sample_size_ec_af = utils.icc_chans_max(np.array(alpha_freq_df_EyesClosed_1_pre), np.array(alpha_freq_df_EyesClosed_2_pre), icc_type)
icc_chan_eo_af, icc_f_eo_af, icc_p_eo_af, icc_lower_eo_af, icc_upper_eo_af, icc_sample_size_eo_af = utils.icc_chans_max(np.array(alpha_freq_df_EyesOpen_1_pre), np.array(alpha_freq_df_EyesOpen_2_pre), icc_type)
icc_chan_ec_aa, icc_f_ec_aa, icc_p_ec_aa, icc_lower_ec_aa, icc_upper_ec_aa, icc_sample_size_ec_aa = utils.icc_chans_max(np.array(alpha_amp_df_EyesClosed_1_pre), np.array(alpha_amp_df_EyesClosed_2_pre), icc_type)
icc_chan_eo_aa, icc_f_eo_aa, icc_p_eo_aa, icc_lower_eo_aa, icc_upper_eo_aa, icc_sample_size_eo_aa = utils.icc_chans_max(np.array(alpha_amp_df_EyesOpen_1_pre), np.array(alpha_amp_df_EyesOpen_2_pre), icc_type)
#Absolute alpha power ICCs
icc_chan_ec_abs_alpha, icc_f_ec_abs_alpha, icc_p_ec_abs_alpha, icc_lower_ec_abs_alpha, icc_upper_ec_abs_alpha, icc_sample_size_ec_abs_alpha = utils.icc_chans_max(np.array(abs_alpha_df_EyesClosed_1_pre), np.array(abs_alpha_df_EyesClosed_2_pre), icc_type)
icc_chan_eo_abs_alpha, icc_f_eo_abs_alpha, icc_p_eo_abs_alpha, icc_lower_eo_abs_alpha, icc_upper_eo_abs_alpha, icc_sample_size_eo_abs_alpha = utils.icc_chans_max(np.array(abs_alpha_df_EyesOpen_1_pre), np.array(abs_alpha_df_EyesOpen_2_pre), icc_type)
#ICCs at the cluster level
icc_type = ["twoway", "agreement", "single"]
icc_ec_exponent_cluster1, icc_f_ec_exponent_cluster1, icc_p_ec_exponent_cluster1, icc_lower_ec_exponent_cluster1, icc_upper_ec_exponent_cluster1, icc_sample_size_ec_exponent_cluster1 = utils.icc_chans_max(np.array(demographic_data_hlm[['Exponent_EyesClosed_1_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), np.array(demographic_data_hlm[['Exponent_EyesClosed_2_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), icc_type)
icc_ec_exponent_cluster2, icc_f_ec_exponent_cluster2, icc_p_ec_exponent_cluster2, icc_lower_ec_exponent_cluster2, icc_upper_ec_exponent_cluster2, icc_sample_size_ec_exponent_cluster2 = utils.icc_chans_max(np.array(demographic_data_hlm[['Exponent_EyesClosed_1_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), np.array(demographic_data_hlm[['Exponent_EyesClosed_2_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), icc_type)
icc_eo_exponent_cluster1, icc_f_eo_exponent_cluster1, icc_p_eo_exponent_cluster1, icc_lower_eo_exponent_cluster1, icc_upper_eo_exponent_cluster1, icc_sample_size_eo_exponent_cluster1 = utils.icc_chans_max(np.array(demographic_data_hlm[['Exponent_EyesOpen_1_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), np.array(demographic_data_hlm[['Exponent_EyesOpen_2_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), icc_type)
icc_eo_exponent_cluster2, icc_f_eo_exponent_cluster2, icc_p_eo_exponent_cluster2, icc_lower_eo_exponent_cluster2, icc_upper_eo_exponent_cluster2, icc_sample_size_eo_exponent_cluster2 = utils.icc_chans_max(np.array(demographic_data_hlm[['Exponent_EyesOpen_1_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), np.array(demographic_data_hlm[['Exponent_EyesOpen_2_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), icc_type)
icc_ec_offset_cluster1, icc_f_ec_offset_cluster1, icc_p_ec_offset_cluster1, icc_lower_ec_offset_cluster1, icc_upper_ec_offset_cluster1, icc_sample_size_ec_offset_cluster1 = utils.icc_chans_max(np.array(demographic_data_hlm[['Offset_EyesClosed_1_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), np.array(demographic_data_hlm[['Offset_EyesClosed_2_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), icc_type)
icc_ec_offset_cluster2, icc_f_ec_offset_cluster2, icc_p_ec_offset_cluster2, icc_lower_ec_offset_cluster2, icc_upper_ec_offset_cluster2, icc_sample_size_ec_offset_cluster2 = utils.icc_chans_max(np.array(demographic_data_hlm[['Offset_EyesClosed_1_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), np.array(demographic_data_hlm[['Offset_EyesClosed_2_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), icc_type)
icc_eo_offset_cluster1, icc_f_eo_offset_cluster1, icc_p_eo_offset_cluster1, icc_lower_eo_offset_cluster1, icc_upper_eo_offset_cluster1, icc_sample_size_eo_offset_cluster1 = utils.icc_chans_max(np.array(demographic_data_hlm[['Offset_EyesOpen_1_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), np.array(demographic_data_hlm[['Offset_EyesOpen_2_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), icc_type)
icc_eo_offset_cluster2, icc_f_eo_offset_cluster2, icc_p_eo_offset_cluster2, icc_lower_eo_offset_cluster2, icc_upper_eo_offset_cluster2, icc_sample_size_eo_offset_cluster2 = utils.icc_chans_max(np.array(demographic_data_hlm[['Offset_EyesOpen_1_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), np.array(demographic_data_hlm[['Offset_EyesOpen_2_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), icc_type)
icc_ec_af_cluster1, icc_f_ec_af_cluster1, icc_p_ec_af_cluster1, icc_lower_ec_af_cluster1, icc_upper_ec_af_cluster1, icc_sample_size_ec_af_cluster1 = utils.icc_chans_max(np.array(demographic_data_hlm[['AlphaFreq_EyesClosed_1_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_1_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), np.array(demographic_data_hlm[['AlphaFreq_EyesClosed_2_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_1_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), icc_type)
icc_ec_af_cluster2, icc_f_ec_af_cluster2, icc_p_ec_af_cluster2, icc_lower_ec_af_cluster2, icc_upper_ec_af_cluster2, icc_sample_size_ec_af_cluster2 = utils.icc_chans_max(np.array(demographic_data_hlm[['AlphaFreq_EyesClosed_1_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_1_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), np.array(demographic_data_hlm[['AlphaFreq_EyesClosed_2_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_1_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), icc_type)
icc_eo_af_cluster1, icc_f_eo_af_cluster1, icc_p_eo_af_cluster1, icc_lower_eo_af_cluster1, icc_upper_eo_af_cluster1, icc_sample_size_eo_af_cluster1 = utils.icc_chans_max(np.array(demographic_data_hlm[['AlphaFreq_EyesOpen_1_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_1_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), np.array(demographic_data_hlm[['AlphaFreq_EyesOpen_2_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_1_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), icc_type)
icc_eo_af_cluster2, icc_f_eo_af_cluster2, icc_p_eo_af_cluster2, icc_lower_eo_af_cluster2, icc_upper_eo_af_cluster2, icc_sample_size_eo_af_cluster2 = utils.icc_chans_max(np.array(demographic_data_hlm[['AlphaFreq_EyesOpen_1_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_1_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), np.array(demographic_data_hlm[['AlphaFreq_EyesOpen_2_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_1_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), icc_type)
icc_ec_aa_cluster1, icc_f_ec_aa_cluster1, icc_p_ec_aa_cluster1, icc_lower_ec_aa_cluster1, icc_upper_ec_aa_cluster1, icc_sample_size_ec_aa_cluster1 = utils.icc_chans_max(np.array(demographic_data_hlm[['AlphaAmp_EyesClosed_1_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_1_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), np.array(demographic_data_hlm[['AlphaAmp_EyesClosed_2_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_1_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), icc_type)
icc_ec_aa_cluster2, icc_f_ec_aa_cluster2, icc_p_ec_aa_cluster2, icc_lower_ec_aa_cluster2, icc_upper_ec_aa_cluster2, icc_sample_size_ec_aa_cluster2 = utils.icc_chans_max(np.array(demographic_data_hlm[['AlphaAmp_EyesClosed_1_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_1_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), np.array(demographic_data_hlm[['AlphaAmp_EyesClosed_2_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_1_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), icc_type)
icc_eo_aa_cluster1, icc_f_eo_aa_cluster1, icc_p_eo_aa_cluster1, icc_lower_eo_aa_cluster1, icc_upper_eo_aa_cluster1, icc_sample_size_eo_aa_cluster1 = utils.icc_chans_max(np.array(demographic_data_hlm[['AlphaAmp_EyesOpen_1_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_1_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), np.array(demographic_data_hlm[['AlphaAmp_EyesOpen_2_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_1_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), icc_type)
icc_eo_aa_cluster2, icc_f_eo_aa_cluster2, icc_p_eo_aa_cluster2, icc_lower_eo_aa_cluster2, icc_upper_eo_aa_cluster2, icc_sample_size_eo_aa_cluster2 = utils.icc_chans_max(np.array(demographic_data_hlm[['AlphaAmp_EyesOpen_1_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_1_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), np.array(demographic_data_hlm[['AlphaAmp_EyesOpen_2_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_1_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), icc_type)
icc_ec_abs_alpha_cluster1, icc_f_ec_abs_alpha_cluster1, icc_p_ec_abs_alpha_cluster1, icc_lower_ec_abs_alpha_cluster1, icc_upper_ec_abs_alpha_cluster1, icc_sample_size_ec_abs_alpha_cluster1 = utils.icc_chans_max(np.array(demographic_data_hlm[['AbsAlpha_EyesClosed_1_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_1_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), np.array(demographic_data_hlm[['AbsAlpha_EyesClosed_2_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_1_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), icc_type)
icc_eo_abs_alpha_cluster1, icc_f_eo_abs_alpha_cluster1, icc_p_eo_abs_alpha_cluster1, icc_lower_eo_abs_alpha_cluster1, icc_upper_eo_abs_alpha_cluster1, icc_sample_size_eo_abs_alpha_cluster1 = utils.icc_chans_max(np.array(demographic_data_hlm[['AbsAlpha_EyesOpen_1_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_1_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), np.array(demographic_data_hlm[['AbsAlpha_EyesOpen_2_pre_Cluster1']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_1_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), icc_type)
icc_ec_abs_alpha_cluster2, icc_f_ec_abs_alpha_cluster2, icc_p_ec_abs_alpha_cluster2, icc_lower_ec_abs_alpha_cluster2, icc_upper_ec_abs_alpha_cluster2, icc_sample_size_ec_abs_alpha_cluster2 = utils.icc_chans_max(np.array(demographic_data_hlm[['AbsAlpha_EyesClosed_1_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_1_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), np.array(demographic_data_hlm[['AbsAlpha_EyesClosed_2_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesClosed_1_pre'] & demographic_data_hlm['Good_Fits_EyesClosed_2_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_1_pre'] & demographic_data_hlm['Good_Alpha_EyesClosed_2_pre'] & demographic_data_hlm['Good_Ref_EyesClosed'] & demographic_data_hlm['Good_Retention_EyesClosed']]), icc_type)
icc_eo_abs_alpha_cluster2, icc_f_eo_abs_alpha_cluster2, icc_p_eo_abs_alpha_cluster2, icc_lower_eo_abs_alpha_cluster2, icc_upper_eo_abs_alpha_cluster2, icc_sample_size_eo_abs_alpha_cluster2 = utils.icc_chans_max(np.array(demographic_data_hlm[['AbsAlpha_EyesOpen_1_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_1_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), np.array(demographic_data_hlm[['AbsAlpha_EyesOpen_2_pre_Cluster2']][demographic_data_hlm['Good_Fits_EyesOpen_1_pre'] & demographic_data_hlm['Good_Fits_EyesOpen_2_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_1_pre'] & demographic_data_hlm['Good_Alpha_EyesOpen_2_pre'] & demographic_data_hlm['Good_Ref_EyesOpen'] & demographic_data_hlm['Good_Retention_EyesOpen']]), icc_type)
icc_cluster_results = pd.DataFrame({
'Measure': ['Exponent', 'Exponent', 'Exponent', 'Exponent', 'Offset', 'Offset', 'Offset', 'Offset', 'AlphaFreq', 'AlphaFreq', 'AlphaFreq', 'AlphaFreq', 'AlphaAmp', 'AlphaAmp', 'AlphaAmp', 'AlphaAmp','AbsAlpha', 'AbsAlpha', 'AbsAlpha', 'AbsAlpha'],
'Task': ['EyesClosed', 'EyesClosed', 'EyesOpen', 'EyesOpen', 'EyesClosed', 'EyesClosed', 'EyesOpen', 'EyesOpen', 'EyesClosed', 'EyesClosed', 'EyesOpen', 'EyesOpen', 'EyesClosed', 'EyesClosed', 'EyesOpen', 'EyesOpen','EyesClosed', 'EyesClosed', 'EyesOpen', 'EyesOpen'],
'Cluster': ['Cluster1', 'Cluster2', 'Cluster1', 'Cluster2', 'Cluster1', 'Cluster2', 'Cluster1', 'Cluster2', 'Cluster1', 'Cluster2', 'Cluster1', 'Cluster2', 'Cluster1', 'Cluster2', 'Cluster1', 'Cluster2','Cluster1', 'Cluster2', 'Cluster1', 'Cluster2'],
'ICC': [icc_ec_exponent_cluster1, icc_ec_exponent_cluster2, icc_eo_exponent_cluster1, icc_eo_exponent_cluster2,
icc_ec_offset_cluster1, icc_ec_offset_cluster2, icc_eo_offset_cluster1, icc_eo_offset_cluster2,
icc_ec_af_cluster1, icc_ec_af_cluster2, icc_eo_af_cluster1, icc_eo_af_cluster2,
icc_ec_aa_cluster1, icc_ec_aa_cluster2, icc_eo_aa_cluster1, icc_eo_aa_cluster2,
icc_ec_abs_alpha_cluster1, icc_ec_abs_alpha_cluster2, icc_eo_abs_alpha_cluster1, icc_eo_abs_alpha_cluster2],
'F-value': [icc_f_ec_exponent_cluster1, icc_f_ec_exponent_cluster2, icc_f_eo_exponent_cluster1, icc_f_eo_exponent_cluster2,
icc_f_ec_offset_cluster1, icc_f_ec_offset_cluster2, icc_f_eo_offset_cluster1, icc_f_eo_offset_cluster2,
icc_f_ec_af_cluster1, icc_f_ec_af_cluster2, icc_f_eo_af_cluster1, icc_f_eo_af_cluster2,
icc_f_ec_aa_cluster1, icc_f_ec_aa_cluster2, icc_f_eo_aa_cluster1, icc_f_eo_aa_cluster2,
icc_f_ec_abs_alpha_cluster1, icc_f_ec_abs_alpha_cluster2, icc_f_eo_abs_alpha_cluster1, icc_f_eo_abs_alpha_cluster2],
'p-value': [icc_p_ec_exponent_cluster1, icc_p_ec_exponent_cluster2, icc_p_eo_exponent_cluster1, icc_p_eo_exponent_cluster2,
icc_p_ec_offset_cluster1, icc_p_ec_offset_cluster2, icc_p_eo_offset_cluster1, icc_p_eo_offset_cluster2,
icc_p_ec_af_cluster1, icc_p_ec_af_cluster2, icc_p_eo_af_cluster1, icc_p_eo_af_cluster2,
icc_p_ec_aa_cluster1, icc_p_ec_aa_cluster2, icc_p_eo_aa_cluster1, icc_p_eo_aa_cluster2,
icc_p_ec_abs_alpha_cluster1, icc_p_ec_abs_alpha_cluster2, icc_p_eo_abs_alpha_cluster1, icc_p_eo_abs_alpha_cluster2],
'Lower CI': [icc_lower_ec_exponent_cluster1, icc_lower_ec_exponent_cluster2, icc_lower_eo_exponent_cluster1, icc_lower_eo_exponent_cluster2,
icc_lower_ec_offset_cluster1, icc_lower_ec_offset_cluster2, icc_lower_eo_offset_cluster1, icc_lower_eo_offset_cluster2,
icc_lower_ec_af_cluster1, icc_lower_ec_af_cluster2, icc_lower_eo_af_cluster1, icc_lower_eo_af_cluster2,
icc_lower_ec_aa_cluster1, icc_lower_ec_aa_cluster2, icc_lower_eo_aa_cluster1, icc_lower_eo_aa_cluster2,
icc_lower_ec_abs_alpha_cluster1, icc_lower_ec_abs_alpha_cluster2, icc_lower_eo_abs_alpha_cluster1, icc_lower_eo_abs_alpha_cluster2],
'Upper CI': [icc_upper_ec_exponent_cluster1, icc_upper_ec_exponent_cluster2, icc_upper_eo_exponent_cluster1, icc_upper_eo_exponent_cluster2,
icc_upper_ec_offset_cluster1, icc_upper_ec_offset_cluster2, icc_upper_eo_offset_cluster1, icc_upper_eo_offset_cluster2,
icc_upper_ec_af_cluster1, icc_upper_ec_af_cluster2, icc_upper_eo_af_cluster1, icc_upper_eo_af_cluster2,
icc_upper_ec_aa_cluster1, icc_upper_ec_aa_cluster2, icc_upper_eo_aa_cluster1, icc_upper_eo_aa_cluster2,
icc_upper_ec_abs_alpha_cluster1, icc_upper_ec_abs_alpha_cluster2, icc_upper_eo_abs_alpha_cluster1, icc_upper_eo_abs_alpha_cluster2],
'N': [icc_sample_size_ec_exponent_cluster1, icc_sample_size_ec_exponent_cluster2, icc_sample_size_eo_exponent_cluster1, icc_sample_size_eo_exponent_cluster2,
icc_sample_size_ec_offset_cluster1, icc_sample_size_ec_offset_cluster2, icc_sample_size_eo_offset_cluster1, icc_sample_size_eo_offset_cluster2,
icc_sample_size_ec_af_cluster1, icc_sample_size_ec_af_cluster2, icc_sample_size_eo_af_cluster1, icc_sample_size_eo_af_cluster2,
icc_sample_size_ec_aa_cluster1, icc_sample_size_ec_aa_cluster2, icc_sample_size_eo_aa_cluster1, icc_sample_size_eo_aa_cluster2,
icc_sample_size_ec_abs_alpha_cluster1, icc_sample_size_ec_abs_alpha_cluster2, icc_sample_size_eo_abs_alpha_cluster1, icc_sample_size_eo_abs_alpha_cluster2]
})
icc_cluster_results['ICC'] = icc_cluster_results['ICC'].apply(lambda x: x[0] if isinstance(x, (list, np.ndarray)) else x)
icc_cluster_results['F-value'] = icc_cluster_results['F-value'].apply(lambda x: x[0] if isinstance(x, (list, np.ndarray)) else x)
icc_cluster_results['p-value'] = icc_cluster_results['p-value'].apply(lambda x: x[0] if isinstance(x, (list, np.ndarray)) else x)
icc_cluster_results['Lower CI'] = icc_cluster_results['Lower CI'].apply(lambda x: x[0] if isinstance(x, (list, np.ndarray)) else x)
icc_cluster_results['Upper CI'] = icc_cluster_results['Upper CI'].apply(lambda x: x[0] if isinstance(x, (list, np.ndarray)) else x)
icc_cluster_results['N'] = icc_cluster_results['N'].apply(lambda x: x[0] if isinstance(x, (list, np.ndarray)) else x)
icc_cluster_results.to_csv(RESULTS_PATH / 'icc_cluster_results_pre.csv', index=False)
# %%
#Plot the ICCs for each measure and task, as a topography
import matplotlib.pyplot as plt
import mne
mapping = 'Blues'
contours = 0
image_interp = 'cubic'
size = 5
vlim = (0, 1)
fig, axs = plt.subplots(2, 4)
im,cm = mne.viz.plot_topomap(icc_chan_ec_exponent, chan_info, axes=axs[0,0], show=False, contours=contours, image_interp=image_interp, size=size, vlim=vlim)
im,cm = mne.viz.plot_topomap(icc_chan_eo_exponent, chan_info, axes=axs[0,1], show=False, contours=contours, image_interp=image_interp, size=size, vlim=vlim)
im,cm = mne.viz.plot_topomap(icc_chan_ec_offset, chan_info, axes=axs[0,2], show=False, contours=contours, image_interp=image_interp, size=size, vlim=vlim)
im,cm = mne.viz.plot_topomap(icc_chan_eo_offset, chan_info, axes=axs[0,3], show=False, contours=contours, image_interp=image_interp, size=size, vlim=vlim)
im,cm = mne.viz.plot_topomap(icc_chan_ec_af, chan_info, axes=axs[1,0], show=False, contours=contours, image_interp=image_interp, size=size, vlim=vlim)
im,cm = mne.viz.plot_topomap(icc_chan_eo_af, chan_info, axes=axs[1,1], show=False, contours=contours, image_interp=image_interp, size=size, vlim=vlim)
im,cm = mne.viz.plot_topomap(icc_chan_ec_aa, chan_info, axes=axs[1,2], show=False, contours=contours, image_interp=image_interp, size=size, vlim=vlim)
im,cm = mne.viz.plot_topomap(icc_chan_eo_aa, chan_info, axes=axs[1,3], show=False, contours=contours, image_interp=image_interp, size=size, vlim=vlim)
# Plot colorbar
ax_x_start = 0.12
ax_y_start = 0.1
ax_width = 0.8
ax_height = 0.02
cbar_ax = fig.add_axes([ax_x_start, ax_y_start, ax_width, ax_height])
clb = fig.colorbar(im, cax=cbar_ax, orientation='horizontal')
clb.set_label('ICC', fontsize=12)
clb.ax.tick_params(labelsize=12)
axs[0,0].set_title('Exponent EC')
axs[0,1].set_title('Exponent EO')
axs[0,2].set_title('Offset EC')
axs[0,3].set_title('Offset EO')
axs[1,0].set_title('IAPF EC')
axs[1,1].set_title('IAPF EO')
axs[1,2].set_title('\u03B1 PW EC')
axs[1,3].set_title('\u03B1 PW EO')
for i in range(0, 4):
axs[0, i].images[0].set_cmap(mapping)
axs[1, i].images[0].set_cmap(mapping)
fig.set_size_inches(8, 5)
fig.subplots_adjust(wspace=0.02, hspace=0.02)
plt.savefig(FIGURES_PATH / 'ICCs_pre.png', dpi=600)
plt.close()
# %%
#Here, we just plot the ICCs for absolute alpha.
import matplotlib.pyplot as plt
import mne
mapping = 'Blues'
contours = 0
image_interp = 'cubic'
size = 5
vlim = (0, 1)
fig, axs = plt.subplots(1, 2)
im,cm = mne.viz.plot_topomap(icc_chan_ec_abs_alpha, chan_info, axes=axs[0], show=False, contours=contours, image_interp=image_interp, size=size, vlim=vlim)
im,cm = mne.viz.plot_topomap(icc_chan_eo_abs_alpha, chan_info, axes=axs[1], show=False, contours=contours, image_interp=image_interp, size=size, vlim=vlim)
# Plot colorbar
ax_x_start = 0.12
ax_y_start = 0.1
ax_width = 0.8
ax_height = 0.02
cbar_ax = fig.add_axes([ax_x_start, ax_y_start, ax_width, ax_height])
clb = fig.colorbar(im, cax=cbar_ax, orientation='horizontal')
clb.set_label('ICC', fontsize=12)
clb.ax.tick_params(labelsize=12)
axs[0].set_title('Absolute \u03B1 Power EC')
axs[1].set_title('Absolute \u03B1 Power EO')
for i in range(0, 1):
axs[0].images[0].set_cmap(mapping)
axs[1].images[0].set_cmap(mapping)
fig.set_size_inches(5, 3)
fig.subplots_adjust(wspace=0.02, hspace=0.02)
plt.savefig(FIGURES_PATH / 'ICCs_abs_alpha_pre.png', dpi=600)
plt.close()
# %%
# Plot the ICCs for each measure and task, as a topography (binned/flat shades)
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import numpy as np
import mne
contours = 0
image_interp = 'cubic'
size = 5
vlim = (0, 1)
bounds = [0.0, 0.4, 0.59, 0.74, 1.0]
tick_locs = [0.2, 0.495, 0.665, 0.875]
tick_labels = ['Poor', 'Fair', 'Good', 'Excellent']
cmap = plt.get_cmap('Blues', 5)
cmap = mcolors.ListedColormap(cmap(np.arange(5)))
norm = mcolors.BoundaryNorm(bounds, cmap.N)
fig, axs = plt.subplots(2, 4)
def _tp(data, ax):
try:
im, cm = mne.viz.plot_topomap(
data, chan_info, axes=ax, show=False,
contours=contours, image_interp=image_interp,
size=size, vlim=vlim, cmap=cmap, cnorm=norm
)
except TypeError:
im, cm = mne.viz.plot_topomap(
data, chan_info, axes=ax, show=False,
contours=contours, image_interp=image_interp,
size=size, vlim=vlim, cmap=cmap
)
im.set_norm(norm)
return im
im = _tp(icc_chan_ec_exponent, axs[0,0])
_tp(icc_chan_eo_exponent, axs[0,1])
_tp(icc_chan_ec_offset, axs[0,2])
_tp(icc_chan_eo_offset, axs[0,3])
_tp(icc_chan_ec_af, axs[1,0])
_tp(icc_chan_eo_af, axs[1,1])
_tp(icc_chan_ec_aa, axs[1,2])
_tp(icc_chan_eo_aa, axs[1,3])
# Colorbar (categorical)
ax_x_start, ax_y_start, ax_width, ax_height = 0.12, 0.10, 0.80, 0.02
cbar_ax = fig.add_axes([ax_x_start, ax_y_start, ax_width, ax_height])
clb = fig.colorbar(im, cax=cbar_ax, orientation='horizontal',
boundaries=bounds, ticks=tick_locs)
clb.set_label('ICC (binned)', fontsize=12)
clb.ax.set_xticklabels(tick_labels)
clb.ax.tick_params(labelsize=12)
# Titles
axs[0,0].set_title('Exponent EC')
axs[0,1].set_title('Exponent EO')
axs[0,2].set_title('Offset EC')
axs[0,3].set_title('Offset EO')
axs[1,0].set_title('IAPF EC')
axs[1,1].set_title('IAPF EO')
axs[1,2].set_title('\u03B1 PW EC')
axs[1,3].set_title('\u03B1 PW EO')
# Size/spacing
fig.set_size_inches(8, 5)
fig.subplots_adjust(wspace=0.02, hspace=0.02)
plt.savefig(FIGURES_PATH / 'ICCs_pre_discrete.png', dpi=600)
plt.close()
# %%
# Plot the ICCs for absolute alpha, as a topography (binned/flat shades)
import matplotlib.pyplot as plt
import mne
mapping = 'Blues'
contours = 0
image_interp = 'nearest'
size = 5
vlim = (0, 1)
fig, axs = plt.subplots(1, 2)
im = _tp(icc_chan_ec_abs_alpha, axs[0])
_tp(icc_chan_eo_abs_alpha, axs[1])
# Colorbar (categorical)
ax_x_start, ax_y_start, ax_width, ax_height = 0.12, 0.10, 0.80, 0.02
cbar_ax = fig.add_axes([ax_x_start, ax_y_start, ax_width, ax_height])
clb = fig.colorbar(im, cax=cbar_ax, orientation='horizontal',
boundaries=bounds, ticks=tick_locs)
clb.set_label('ICC (binned)', fontsize=12)
clb.ax.set_xticklabels(tick_labels)
clb.ax.tick_params(labelsize=12)
axs[0].set_title('Absolute \u03B1 Power EC')
axs[1].set_title('Absolute \u03B1 Power EO')
fig.set_size_inches(5, 3)
fig.subplots_adjust(wspace=0.02, hspace=0.02)
plt.savefig(FIGURES_PATH / 'ICCs_abs_alpha_pre_discrete.png', dpi=600)
plt.close()
# %% [markdown]
# ## 6. Final Sample Details, Baseline, Correlations and HLMs
# These analyses are all handled at the cluster level. Aside from some data-wrangling, the actual analyses are in the IfAdo_analyses rmarkdown notebook.
# %%
#If we don't have demographic_data_hlm, reload it from the csv
if 'demographic_data_hlm' not in locals():
demographic_data_hlm = pd.read_csv(RESULTS_PATH / 'demographic_data_hlm_pre.csv', index_col=0)
df_hlm = demographic_data_hlm.reset_index()
# %%
# Create a filtered dataframe with participants who have good fits
df_hlm_filtered = df_hlm.copy()
# For participants with session2 = 'yes', require good fits in all 4 pre tasks
# For participants with session2 = 'no', require good fits in session 1 pre tasks
condition_session2_yes = (
(df_hlm_filtered['session2'] == 'yes') &
(df_hlm_filtered['Good_Fits_EyesClosed_1_pre'] == True) &
(df_hlm_filtered['Good_Fits_EyesClosed_2_pre'] == True) &
(df_hlm_filtered['Good_Fits_EyesOpen_1_pre'] == True) &
(df_hlm_filtered['Good_Fits_EyesOpen_2_pre'] == True) &
(df_hlm_filtered['Good_Retention_EyesClosed'] == True) &
(df_hlm_filtered['Good_Retention_EyesOpen'] == True) &
(df_hlm_filtered['Good_Ref_EyesClosed'] == True) &
(df_hlm_filtered['Good_Ref_EyesOpen'] == True)
)
condition_session2_no = (
(df_hlm_filtered['session2'] == 'no') &
(df_hlm_filtered['Good_Fits_EyesClosed_1_pre'] == True) &
(df_hlm_filtered['Good_Fits_EyesOpen_1_pre'] == True)
)
df_hlm_good_fits = df_hlm_filtered[condition_session2_yes | condition_session2_no].copy()
print(f"Original dataframe: {len(df_hlm_filtered)} rows, {df_hlm_filtered['participant_id'].nunique()} participants")
print(f"Filtered dataframe: {len(df_hlm_good_fits)} rows, {df_hlm_good_fits['participant_id'].nunique()} participants")
session2_yes_count = df_hlm_good_fits[df_hlm_good_fits['session2'] == 'yes']['participant_id'].nunique()
session2_no_count = df_hlm_good_fits[df_hlm_good_fits['session2'] == 'no']['participant_id'].nunique()
print(f"Participants with session2 = 'yes': {session2_yes_count}")
print(f"Participants with session2 = 'no': {session2_no_count}")
print(f"Total participants with good fits: {session2_yes_count + session2_no_count}")
# %%
import re
import pandas as pd
import numpy as np
def make_hlm_df(measure, condition, cluster=None, df=df_hlm, require_session2=True, require_good_fits=True, require_alpha_fits = False, dropna=False):
df = df_hlm
m = str(measure).strip().lower()
measure_map = {
'exponent': 'Exponent',
'offset': 'Offset',
'alpha_freq': 'AlphaFreq', 'alphafreq': 'AlphaFreq', 'alpha-frequency': 'AlphaFreq',
'alpha_amp': 'AlphaAmp', 'alphaamp': 'AlphaAmp', 'alpha-amplitude': 'AlphaAmp', 'absolute-alpha': 'AbsAlpha',
'fit': 'Fit', 'error': 'Error'
}
measure_key = measure_map.get(m, None)
if measure_key is None:
measure_key = str(measure).strip().replace(' ', '').title()
c = str(condition).strip().lower()
if 'closed' in c:
cond_key = 'EyesClosed'
elif 'open' in c:
cond_key = 'EyesOpen'
else:
cond_key = str(condition).strip().replace(' ', '').title()
good1 = f"Good_Fits_{cond_key}_1_pre"
good2 = f"Good_Fits_{cond_key}_2_pre"
good3 = f"Good_Retention_{cond_key}"
good4 = f"Good_Ref_{cond_key}"
for rc in (good1, good2, good3, good4):
if rc not in df.columns:
df[rc] = False
if require_alpha_fits:
good1_alpha = f"Good_Alpha_{cond_key}_1_pre"
good2_alpha = f"Good_Alpha_{cond_key}_2_pre"
if good1_alpha not in df.columns:
df[good1_alpha] = False
if good2_alpha not in df.columns:
df[good2_alpha] = False
# find all cluster-specific columns for the measure+condition across timepoints
# expected pattern: <Measure>_<Condition>_<time>_pre_Cluster<NUM>
pat = re.compile(rf"^{re.escape(measure_key)}_{re.escape(cond_key)}_[12]_pre_Cluster(\d+)$", flags=re.IGNORECASE)
matched_cols = [col for col in df.columns if pat.match(col)]
if cluster is not None:
if isinstance(cluster, int):
want = f"Cluster{cluster}"
else:
cs = str(cluster).strip()
want = cs if cs.lower().startswith('cluster') else f"Cluster{cs}"
matched_cols = [c for c in matched_cols if c.endswith(want)]
if not matched_cols:
col_t1 = f"{measure_key}_{cond_key}_1_pre_Cluster{cluster if cluster is not None else ''}".rstrip('_')
col_t2 = f"{measure_key}_{cond_key}_2_pre_Cluster{cluster if cluster is not None else ''}".rstrip('_')
for rc in (col_t1, col_t2):
if rc not in df.columns:
df[rc] = np.nan
matched_cols = [c for c in (col_t1, col_t2) if c in df.columns]
required_ids = ['participant_id', 'age', 'sex', 'session2', good1, good2, good3, good4]
if require_alpha_fits:
required_ids += [good1_alpha, good2_alpha]
for rc in required_ids:
if rc not in df.columns:
if rc in (good1, good2, good3, good4) + (() if not require_alpha_fits else (good1_alpha, good2_alpha)):
df[rc] = False
else:
df[rc] = np.nan
# melt across all matched cluster columns
id_vars = ['participant_id', 'age', 'sex', 'session2', good1, good2, good3, good4]
if require_alpha_fits:
id_vars += [good1_alpha, good2_alpha]
melted = df[id_vars + matched_cols].melt(
id_vars=id_vars,
value_vars=matched_cols,
var_name='orig_measure_col',
value_name='value'
)
# extract time (1 or 2) from the original column name
melted['time'] = melted['orig_measure_col'].str.extract(r'_(1|2)_pre_', expand=False).astype(float)
# extract cluster number and make it nullable integer
melted['cluster'] = pd.to_numeric(melted['orig_measure_col'].str.extract(r'Cluster(\d+)', expand=False), errors='coerce').astype('Int64')
# center age
melted['age'] = pd.to_numeric(melted['age'], errors='coerce')
melted['age_c'] = melted['age'] - melted['age'].mean()
# set measure name to canonical measure_key
melted['measure'] = measure_key
# optional filtering
if require_session2 and require_good_fits:
melted = melted[
(melted['session2'] == 'yes') &
(melted[good1] == True) &
(melted[good2] == True) &
(melted[good3] == True) & (melted[good4] == True)
].copy()
elif require_session2:
melted = melted[melted['session2'] == 'yes'].copy()
elif require_good_fits:
melted = melted[(melted[good1] == True) & (melted[good2] == True) & (melted[good3] == True)].copy()
if require_alpha_fits:
melted = melted[
(melted[good1_alpha] == True) &
(melted[good2_alpha] == True)
].copy()
if dropna:
melted = melted.dropna(subset=['value'])
# drop the helper original column before returning
melted = melted.drop(columns=['orig_measure_col'])
return melted
# %%
#Create data frames for each measure, task, and cluster, using make_hlm_df function. Note that we don't need to specify require session 2, or good fits as they're defaults
#Exponents
exponents_both_clusters_eyesclosed_hlm = make_hlm_df('Exponent', 'EyesClosed', cluster=None, dropna=True)
exponents_both_clusters_eyesopen_hlm = make_hlm_df('Exponent', 'EyesOpen', cluster=None, dropna=True)
exponent_cluster1_eyesclosed_hlm = make_hlm_df('Exponent', 'EyesClosed', 'Cluster1',dropna=True)
exponent_cluster2_eyesclosed_hlm = make_hlm_df('Exponent', 'EyesClosed', 'Cluster2', dropna=True)
exponent_cluster1_eyesopen_hlm = make_hlm_df('Exponent', 'EyesOpen', 'Cluster1', dropna=True)
exponent_cluster2_eyesopen_hlm = make_hlm_df('Exponent', 'EyesOpen', 'Cluster2', dropna=True)
#Offsets
offset_cluster1_eyesclosed_hlm = make_hlm_df('Offset', 'EyesClosed', 'Cluster1', dropna=True)
offset_cluster2_eyesclosed_hlm = make_hlm_df('Offset', 'EyesClosed', 'Cluster2', dropna=True)
offset_cluster1_eyesopen_hlm = make_hlm_df('Offset', 'EyesOpen', 'Cluster1', dropna=True)
offset_cluster2_eyesopen_hlm = make_hlm_df('Offset', 'EyesOpen', 'Cluster2', dropna=True)
#alpha frequency
alpha_freq_eyesclosed_hlm = make_hlm_df('AlphaFreq', 'EyesClosed', cluster=None, require_alpha_fits=True, dropna=True)
alpha_freq_eyesopen_hlm = make_hlm_df('AlphaFreq', 'EyesOpen', cluster=None, require_alpha_fits=True, dropna=True)
alpha_freq_eyesclosed_cluster1_hlm = make_hlm_df('AlphaFreq', 'EyesClosed', 'Cluster1', require_alpha_fits=True, dropna=True)
alpha_freq_eyesclosed_cluster2_hlm = make_hlm_df('AlphaFreq', 'EyesClosed', 'Cluster2', require_alpha_fits=True, dropna=True)
alpha_freq_eyesopen_cluster1_hlm = make_hlm_df('AlphaFreq', 'EyesOpen', 'Cluster1', require_alpha_fits=True, dropna=True)
alpha_freq_eyesopen_cluster2_hlm = make_hlm_df('AlphaFreq', 'EyesOpen', 'Cluster2', require_alpha_fits=True, dropna=True)
#alpha amplitude
alpha_amp_eyesclosed_hlm = make_hlm_df('AlphaAmp', 'EyesClosed', cluster=None, require_alpha_fits=True, dropna=True)
alpha_amp_eyesopen_hlm = make_hlm_df('AlphaAmp', 'EyesOpen', cluster=None, require_alpha_fits=True, dropna=True)
alpha_amp_eyesclosed_cluster1_hlm = make_hlm_df('AlphaAmp', 'EyesClosed', 'Cluster1', require_alpha_fits=True, dropna=True)
alpha_amp_eyesclosed_cluster2_hlm = make_hlm_df('AlphaAmp', 'EyesClosed', 'Cluster2', require_alpha_fits=True, dropna=True)
alpha_amp_eyesopen_cluster1_hlm = make_hlm_df('AlphaAmp', 'EyesOpen', 'Cluster1', require_alpha_fits=True, dropna=True)
alpha_amp_eyesopen_cluster2_hlm = make_hlm_df('AlphaAmp', 'EyesOpen', 'Cluster2', require_alpha_fits=True, dropna=True)
#Absolute alpha
abs_alpha_eyesclosed_hlm = make_hlm_df('AbsAlpha', 'EyesClosed', cluster=None, require_alpha_fits=True, dropna=True)
abs_alpha_eyesopen_hlm = make_hlm_df('AbsAlpha', 'EyesOpen', cluster=None, require_alpha_fits=True, dropna=True)
abs_alpha_eyesclosed_cluster1_hlm = make_hlm_df('AbsAlpha', 'EyesClosed', 'Cluster1', require_alpha_fits=True, dropna=True)
abs_alpha_eyesclosed_cluster2_hlm = make_hlm_df('AbsAlpha', 'EyesClosed', 'Cluster2', require_alpha_fits=True, dropna=True)
abs_alpha_eyesopen_cluster1_hlm = make_hlm_df('AbsAlpha', 'EyesOpen', 'Cluster1', require_alpha_fits=True, dropna=True)
abs_alpha_eyesopen_cluster2_hlm = make_hlm_df('AbsAlpha', 'EyesOpen', 'Cluster2', require_alpha_fits=True, dropna=True)
hlm_dfs = {
'exponent_cluster1_eyesclosed': exponent_cluster1_eyesclosed_hlm,
'exponent_cluster2_eyesclosed': exponent_cluster2_eyesclosed_hlm,
'exponent_cluster1_eyesopen': exponent_cluster1_eyesopen_hlm,
'exponent_cluster2_eyesopen': exponent_cluster2_eyesopen_hlm,
'offset_cluster1_eyesclosed': offset_cluster1_eyesclosed_hlm,
'offset_cluster2_eyesclosed': offset_cluster2_eyesclosed_hlm,
'offset_cluster1_eyesopen': offset_cluster1_eyesopen_hlm,
'offset_cluster2_eyesopen': offset_cluster2_eyesopen_hlm,
'alpha_freq_cluster1_eyesclosed': alpha_freq_eyesclosed_cluster1_hlm,
'alpha_freq_cluster2_eyesclosed': alpha_freq_eyesclosed_cluster2_hlm,
'alpha_freq_cluster1_eyesopen': alpha_freq_eyesopen_cluster1_hlm,
'alpha_freq_cluster2_eyesopen': alpha_freq_eyesopen_cluster2_hlm,
'alpha_amp_cluster1_eyesclosed': alpha_amp_eyesclosed_cluster1_hlm,
'alpha_amp_cluster2_eyesclosed': alpha_amp_eyesclosed_cluster2_hlm,
'alpha_amp_cluster1_eyesopen': alpha_amp_eyesopen_cluster1_hlm,
'alpha_amp_cluster2_eyesopen': alpha_amp_eyesopen_cluster2_hlm,
'abs_alpha_cluster1_eyesclosed': abs_alpha_eyesclosed_cluster1_hlm,
'abs_alpha_cluster2_eyesclosed': abs_alpha_eyesclosed_cluster2_hlm,
'abs_alpha_cluster1_eyesopen': abs_alpha_eyesopen_cluster1_hlm,
'abs_alpha_cluster2_eyesopen': abs_alpha_eyesopen_cluster2_hlm
}
hlm_dfs_all = {
'exponents_both_clusters_eyesclosed': exponents_both_clusters_eyesclosed_hlm,
'exponents_both_clusters_eyesopen': exponents_both_clusters_eyesopen_hlm,
'offsets_both_clusters_eyesclosed': make_hlm_df('Offset', 'EyesClosed', cluster=None, dropna=True),
'offsets_both_clusters_eyesopen': make_hlm_df('Offset', 'EyesOpen', cluster=None, dropna=True),
'alpha_freq_both_clusters_eyesclosed': alpha_freq_eyesclosed_hlm,
'alpha_freq_both_clusters_eyesopen': alpha_freq_eyesopen_hlm,
'alpha_amp_both_clusters_eyesclosed': alpha_amp_eyesclosed_hlm,
'alpha_amp_both_clusters_eyesopen': alpha_amp_eyesopen_hlm,
'abs_alpha_both_clusters_eyesclosed': abs_alpha_eyesclosed_hlm,
'abs_alpha_both_clusters_eyesopen': abs_alpha_eyesopen_hlm
}
# %%
# Sample size diagnostics - count participants for each analysis type
from contextlib import redirect_stdout
#Save the printed info to a text file
with open('sample_size_demographics.txt', 'w') as f:
with redirect_stdout(f):
print("Original Sample Size and Demographics")
print("=" * 60)
# 1. Overall sample sizes
print("\n1. Sample Size")
print("-" * 30)
print(f"Total participants in dataset: {len(demographic_data)}")
print(f"Participants with session2 = 'yes': {(demographic_data['session2'] == 'yes').sum()}")
print(f"Participants with session2 = 'no': {(demographic_data['session2'] == 'no').sum()}")
print(f"Age: M={demographic_data['age'].mean():.2f}, SD={demographic_data['age'].std():.2f}, Range={demographic_data['age'].min():.0f}-{demographic_data['age'].max():.0f}")
print(f"Sex: {(demographic_data['sex'] == 'F').sum()} Female, {(demographic_data['sex'] == 'M').sum()} Male")
# 2. Specparam fits sample sizes
print("\n2. Specparam Fits - overall, not means")
print("-" * 30)
for task in ['EyesClosed', 'EyesOpen']:
for session in ['1', '2']:
col_name = f'Good_Fits_{task}_{session}_pre'
if col_name in demographic_data.columns:
count = demographic_data[col_name].sum()
print(f"Good fits - {task} Session {session}: {count}")
# 2.1 Mean and SD of R^2 for each task and session, for those with good fits Mean_Fits_EyesClosed_1_pre
print("\n2.1 Mean and SD of R^2 for each task and session (good fits only)")
print("-" * 50)
for task in ['EyesClosed', 'EyesOpen']:
for session in ['1', '2']:
col_good = f'Good_Fits_{task}_{session}_pre'
col_r2 = f'Mean_Fits_{task}_{session}_pre'
if col_good in demographic_data.columns and col_r2 in demographic_data.columns:
good_r2 = demographic_data.loc[demographic_data[col_good] == True, col_r2]
mean_r2 = good_r2.mean()
sd_r2 = good_r2.std()
print(f"{task} Session {session} - Mean R^2: {mean_r2:.4f}, SD R^2: {sd_r2:.4f}")
#2.2 Mean and SD of error for each task and session, for those with good fits
print("\n2.2 Mean and SD of Error for each task and session (good fits only)")
print("-" * 50)
for task in ['EyesClosed', 'EyesOpen']:
for session in ['1', '2']:
col_good = f'Good_Fits_{task}_{session}_pre'
col_error = f'Mean_Error_{task}_{session}_pre'
if col_good in demographic_data.columns and col_error in demographic_data.columns:
good_error = demographic_data.loc[demographic_data[col_good] == True, col_error]
mean_error = good_error.mean()
sd_error = good_error.std()
print(f"{task} Session {session} - Mean Error: {mean_error:.4f}, SD Error: {sd_error:.4f}")
# Count participants with good fits in all 4 pre tasks (our main analysis sample)
all_good_fits = (
(demographic_data['Good_Fits_EyesClosed_1_pre'] == True) &
(demographic_data['Good_Fits_EyesClosed_2_pre'] == True) &
(demographic_data['Good_Fits_EyesOpen_1_pre'] == True) &
(demographic_data['Good_Fits_EyesOpen_2_pre'] == True) &
(demographic_data['session2'] == 'yes')
)
# 3. Alpha peak detection sample sizes
print("\n3. Alpha peak detection - overall, not means")
print("-" * 30)
for task in ['EyesClosed', 'EyesOpen']:
for session in ['1', '2']:
col_name = f'Good_Alpha_{task}_{session}_pre'
if col_name in demographic_data.columns:
count = demographic_data[col_name].sum()
print(f"Good alpha peaks - {task} Session {session}: {count}")
all_good_alpha = (
(demographic_data['Good_Alpha_EyesClosed_1_pre'] == True) &
(demographic_data['Good_Alpha_EyesClosed_2_pre'] == True) &
(demographic_data['Good_Alpha_EyesOpen_1_pre'] == True) &
(demographic_data['Good_Alpha_EyesOpen_2_pre'] == True) &
(demographic_data['session2'] == 'yes')
)
# 4. Task-specific breakdowns for session2=yes participants
print("\n4. Task specific sample sizes (session2 = 'yes' only)")
print("-" * 50)
session2_yes = demographic_data[demographic_data['session2'] == 'yes'].copy()
for task in ['EyesClosed', 'EyesOpen']:
print(f"\n{task}:")
# Both timepoints good
both_good_fits = (
(session2_yes[f'Good_Fits_{task}_1_pre'] == True) &
(session2_yes[f'Good_Fits_{task}_2_pre'] == True)
)
both_good_alpha = (
(session2_yes[f'Good_Alpha_{task}_1_pre'] == True) &
(session2_yes[f'Good_Alpha_{task}_2_pre'] == True)
)
print(f" Good fits both timepoints: {both_good_fits.sum()}")
print(f" Good alpha both timepoints: {both_good_alpha.sum()}")
print(f" Both criteria both timepoints: {(both_good_fits & both_good_alpha).sum()}")
#Put chan stuff here.
for task in ['EyesClosed', 'EyesOpen']:
print(f"\n{task} - Channels removed during preprocessing (session2 = 'yes' only):")
for session in ['1', '2']:
col_name = f'Channels_Removed_task-{task}_ses-{session}_acq-pre'
if col_name in session2_yes.columns:
removed_counts = session2_yes[col_name].dropna()
removed_counts = removed_counts[removed_counts != 99]
mean_removed = removed_counts.mean()
sd_removed = removed_counts.std()
max_removed = removed_counts.max()
print(f" Session {session}: Mean channels removed: {mean_removed:.2f}, SD: {sd_removed:.2f}, Max: {max_removed}")
# 5. Demographics for the samples used in the HLMs. Note that these are all session2.
print("\n5. Demographics for analysis samples (session2 = 'yes' only), with fit/error info")
print("-" * 40)
for task in ['EyesClosed', 'EyesOpen']:
print(f"\n{task}:")
task_sample = session2_yes[session2_yes[f'Good_Fits_{task}_1_pre'] & session2_yes[f'Good_Fits_{task}_2_pre'] & session2_yes[f'Good_Ref_{task}'] & session2_yes[f'Good_Retention_{task}']].copy()
alpha_task_sample = session2_yes[session2_yes[f'Good_Fits_{task}_1_pre'] & session2_yes[f'Good_Fits_{task}_2_pre'] & session2_yes[f'Good_Alpha_{task}_1_pre'] & session2_yes[f'Good_Alpha_{task}_2_pre'] & session2_yes[f'Good_Ref_{task}'] & session2_yes[f'Good_Retention_{task}']].copy()
if len(task_sample) > 0:
print(f" Participants in exponent/offset HLMs (n={len(task_sample)}):")
print(f" Age: M={task_sample['age'].mean():.2f}, SD={task_sample['age'].std():.2f}, Range={task_sample['age'].min():.0f}-{task_sample['age'].max():.0f}")
print(f" Sex: {(task_sample['sex'] == 'F').sum()} Female, {(task_sample['sex'] == 'M').sum()} Male")
#mean fits and error
print(f" Mean R^2 Session 1: {task_sample[f'Mean_Fits_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].mean():.2f}, SD: {task_sample[f'Mean_Fits_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].std():.2f}")
print(f" Mean R^2 Session 2: {task_sample[f'Mean_Fits_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].mean():.2f}, SD: {task_sample[f'Mean_Fits_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].std():.2f}")
print(f" Fit Range Session 1: {task_sample[f'Mean_Fits_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].min():.2f} - {task_sample[f'Mean_Fits_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].max():.2f}")
print(f" Fit Range Session 2: {task_sample[f'Mean_Fits_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].min():.2f} - {task_sample[f'Mean_Fits_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].max():.2f}")
print(f" Mean Error Session 1: {task_sample[f'Mean_Error_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].mean():.2f}, SD: {task_sample[f'Mean_Error_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].std():.2f}")
print(f" Mean Error Session 2: {task_sample[f'Mean_Error_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].mean():.2f}, SD: {task_sample[f'Mean_Error_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].std():.2f}")
print(f" Error Range Session 1: {task_sample[f'Mean_Error_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre']].min():.2f} - {task_sample[f'Mean_Error_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre']].max():.2f}")
print(f" Error Range Session 2: {task_sample[f'Mean_Error_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].min():.2f} - {task_sample[f'Mean_Error_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].max():.2f}")
#And using the interpolated means
print(f" Mean R^2 Session 1 (interp): {task_sample[f'Mean_Fits_interp_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].mean():.2f}, SD: {task_sample[f'Mean_Fits_interp_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].std():.2f}")
print(f" Mean R^2 Session 2 (interp): {task_sample[f'Mean_Fits_interp_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].mean():.2f}, SD: {task_sample[f'Mean_Fits_interp_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].std():.2f}")
print(f" Fit Range Session 1 (interp): {task_sample[f'Mean_Fits_interp_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].min():.3f} - {task_sample[f'Mean_Fits_interp_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].max():.2f}")
print(f" Fit Range Session 2 (interp): {task_sample[f'Mean_Fits_interp_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].min():.2f} - {task_sample[f'Mean_Fits_interp_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].max():.2f}")
print(f" Mean Error Session 1 (interp): {task_sample[f'Mean_Error_interp_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].mean():.2f}, SD: {task_sample[f'Mean_Error_interp_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].std():.2f}")
print(f" Mean Error Session 2 (interp): {task_sample[f'Mean_Error_interp_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].mean():.2f}, SD: {task_sample[f'Mean_Error_interp_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].std():.2f}")
print(f" Error Range Session 1 (interp): {task_sample[f'Mean_Error_interp_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].min():.2f} - {task_sample[f'Mean_Error_interp_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].max():.2f}")
print(f" Error Range Session 2 (interp): {task_sample[f'Mean_Error_interp_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].min():.2f} - {task_sample[f'Mean_Error_interp_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].max():.2f}")
#And Epoch counts
print(f" Original Epochs Session 1: {task_sample[f'EpochsOriginal_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].mean():.2f}, SD: {task_sample[f'EpochsOriginal_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].std():.2f}")
print(f" Original Epochs Session 2: {task_sample[f'EpochsOriginal_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].mean():.2f}, SD: {task_sample[f'EpochsOriginal_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].std():.2f}")
print(f" Cleaned Epochs Session 1: {task_sample[f'EpochsRetained_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].mean():.2f}, SD: {task_sample[f'EpochsRetained_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].std():.2f}")
print(f" Cleaned Epochs Session 2: {task_sample[f'EpochsRetained_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].mean():.2f}, SD: {task_sample[f'EpochsRetained_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].std():.2f}")
#Min and Max Epochs
print(f" Original Epochs Range Session 1: {task_sample[f'EpochsOriginal_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].min():.0f} - {task_sample[f'EpochsOriginal_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].max():.0f}")
print(f" Original Epochs Range Session 2: {task_sample[f'EpochsOriginal_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].min():.0f} - {task_sample[f'EpochsOriginal_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].max():.0f}")
print(f" Cleaned Epochs Range Session 1: {task_sample[f'EpochsRetained_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].min():.0f} - {task_sample[f'EpochsRetained_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].max():.0f}")
print(f" Cleaned Epochs Range Session 2: {task_sample[f'EpochsRetained_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].min():.0f} - {task_sample[f'EpochsRetained_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].max():.0f}")
#Proportion of epochs retained
print(f" Proportion of Epochs Retained Session 1: {task_sample[f'EpochsProportion_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].mean()*100:.2f}, SD: {task_sample[f'EpochsProportion_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].std()*100:.2f}")
print(f" Proportion of Epochs Retained Session 2: {task_sample[f'EpochsProportion_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].mean()*100:.2f}, SD: {task_sample[f'EpochsProportion_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].std()*100:.2f}")
#min and max proportion
print(f" Proportion of Epochs Retained Range Session 1: {task_sample[f'EpochsProportion_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].min()*100:.2f} - {task_sample[f'EpochsProportion_{task}_1_pre'][task_sample[f'Good_Fits_{task}_1_pre'] == True].max()*100:.2f}")
print(f" Proportion of Epochs Retained Range Session 2: {task_sample[f'EpochsProportion_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].min()*100:.2f} - {task_sample[f'EpochsProportion_{task}_2_pre'][task_sample[f'Good_Fits_{task}_2_pre'] == True].max()*100:.2f}")
if len(alpha_task_sample) > 0:
print(f" Participants in alpha HLMs (n={len(alpha_task_sample)}):")
print(f" Age: M={alpha_task_sample['age'].mean():.2f}, SD={alpha_task_sample['age'].std():.2f}, Range={alpha_task_sample['age'].min():.0f}-{alpha_task_sample['age'].max():.0f}")
print(f" Sex: {(alpha_task_sample['sex'] == 'F').sum()} Female, {(alpha_task_sample['sex'] == 'M').sum()} Male")
print(f" Mean R^2 Session 1: {alpha_task_sample[f'Mean_Fits_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].mean():.2f}, SD: {alpha_task_sample[f'Mean_Fits_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].std():.2f}")
print(f" Mean R^2 Session 2: {alpha_task_sample[f'Mean_Fits_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].mean():.2f}, SD: {alpha_task_sample[f'Mean_Fits_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].std():.2f}")
print(f" Fit Range Session 1: {alpha_task_sample[f'Mean_Fits_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].min():.2f} - {alpha_task_sample[f'Mean_Fits_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].max():.2f}")
print(f" Fit Range Session 2: {alpha_task_sample[f'Mean_Fits_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].min():.2f} - {alpha_task_sample[f'Mean_Fits_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].max():.2f}")
print(f" Mean Error Session 1: {alpha_task_sample[f'Mean_Error_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].mean():.2f}, SD: {alpha_task_sample[f'Mean_Error_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].std():.2f}")
print(f" Mean Error Session 2: {alpha_task_sample[f'Mean_Error_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].mean():.2f}, SD: {alpha_task_sample[f'Mean_Error_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].std():.2f}")
print(f" Error Range Session 1: {alpha_task_sample[f'Mean_Error_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].min():.2f} - {alpha_task_sample[f'Mean_Error_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].max():.2f}")
print(f" Error Range Session 2: {alpha_task_sample[f'Mean_Error_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].min():.2f} - {alpha_task_sample[f'Mean_Error_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].max():.2f}")
print(f" Mean R^2 Session 1 (interp): {alpha_task_sample[f'Mean_Fits_interp_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].mean():.2f}, SD: {alpha_task_sample[f'Mean_Fits_interp_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].std():.2f}")
print(f" Mean R^2 Session 2 (interp): {alpha_task_sample[f'Mean_Fits_interp_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].mean():.2f}, SD: {alpha_task_sample[f'Mean_Fits_interp_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].std():.2f}")
print(f" Fit Range Session 1 (interp): {alpha_task_sample[f'Mean_Fits_interp_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].min():.2f} - {alpha_task_sample[f'Mean_Fits_interp_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].max():.2f}")
print(f" Fit Range Session 2 (interp): {alpha_task_sample[f'Mean_Fits_interp_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].min():.2f} - {alpha_task_sample[f'Mean_Fits_interp_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].max():.2f}")
print(f" Mean Error Session 1 (interp): {alpha_task_sample[f'Mean_Error_interp_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].mean():.2f}, SD: {alpha_task_sample[f'Mean_Error_interp_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].std():.2f}")
print(f" Mean Error Session 2 (interp): {alpha_task_sample[f'Mean_Error_interp_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].mean():.2f}, SD: {alpha_task_sample[f'Mean_Error_interp_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].std():.2f}")
print(f" Error Range Session 1 (interp): {alpha_task_sample[f'Mean_Error_interp_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].min():.2f} - {alpha_task_sample[f'Mean_Error_interp_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].max():.2f}")
print(f" Error Range Session 2 (interp): {alpha_task_sample[f'Mean_Error_interp_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].min():.2f} - {alpha_task_sample[f'Mean_Error_interp_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].max():.2f}")
#And Epoch counts
print(f" Original Epochs Session 1: {alpha_task_sample[f'EpochsOriginal_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].mean():.2f}, SD: {alpha_task_sample[f'EpochsOriginal_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].std():.2f}")
print(f" Original Epochs Session 2: {alpha_task_sample[f'EpochsOriginal_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].mean():.2f}, SD: {alpha_task_sample[f'EpochsOriginal_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].std():.2f}")
print(f" Cleaned Epochs Session 1: {alpha_task_sample[f'EpochsRetained_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].mean():.2f}, SD: {alpha_task_sample[f'EpochsRetained_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].std():.2f}")
print(f" Cleaned Epochs Session 2: {alpha_task_sample[f'EpochsRetained_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].mean():.2f}, SD: {alpha_task_sample[f'EpochsRetained_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].std():.2f}")
#Min and Max Epochs
print(f" Original Epochs Range Session 1: {alpha_task_sample[f'EpochsOriginal_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].min():.0f} - {alpha_task_sample[f'EpochsOriginal_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].max():.0f}")
print(f" Original Epochs Range Session 2: {alpha_task_sample[f'EpochsOriginal_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].min():.0f} - {alpha_task_sample[f'EpochsOriginal_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].max():.0f}")
print(f" Cleaned Epochs Range Session 1: {alpha_task_sample[f'EpochsRetained_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].min():.0f} - {alpha_task_sample[f'EpochsRetained_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].max():.0f}")
print(f" Cleaned Epochs Range Session 2: {alpha_task_sample[f'EpochsRetained_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].min():.0f} - {alpha_task_sample[f'EpochsRetained_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].max():.0f}")
#Proportion of epochs retained
print(f" Proportion of Epochs Retained Session 1: {alpha_task_sample[f'EpochsProportion_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].mean()*100:.2f}, SD: {alpha_task_sample[f'EpochsProportion_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].std()*100:.2f}")
print(f" Proportion of Epochs Retained Session 2: {alpha_task_sample[f'EpochsProportion_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].mean()*100:.2f}, SD: {alpha_task_sample[f'EpochsProportion_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].std()*100:.2f}")
#min and max proportion
print(f" Proportion of Epochs Retained Range Session 1: {alpha_task_sample[f'EpochsProportion_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].min()*100:.2f} - {alpha_task_sample[f'EpochsProportion_{task}_1_pre'][alpha_task_sample[f'Good_Fits_{task}_1_pre']].max()*100:.2f}")
print(f" Proportion of Epochs Retained Range Session 2: {alpha_task_sample[f'EpochsProportion_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].min()*100:.2f} - {alpha_task_sample[f'EpochsProportion_{task}_2_pre'][alpha_task_sample[f'Good_Fits_{task}_2_pre']].max()*100:.2f}")
pass
# %% [markdown]
# ### 6.2 HLMs predicting FOOOF parameters from time, age and sex are implemented in R.
# %%
#Save each of the hlm dataframes to a csv file
for name, hlm_df in hlm_dfs.items():
hlm_df.to_csv(RESULTS_PATH / f'{name}_for_hlm.csv', index=False)

IfAdo_preprocessing.ipynb at commit e056acf, under MIT · at the source

Overview

Authors: Polina Politanskaia¹, Jacinta Bywater¹, Anna J Finley², Hannah A D Keage³, Nicholas J Kelley⁴, Daniel J McKeown¹, Victor R Schinazi¹, Douglas J Angus¹

ORCID iDs: Daniel J McKeown, Douglas J Angus

Institutions: Bond University (Australia); North Dakota State University (United States); Adelaide University (Australia); The University of Adelaide (Australia); University of Southampton (United Kingdom)

Journal: Cerebral cortex (New York, N.Y. : 1991), volume 36, issue 7, article bhag113

Dates: received 3 March 2026; accepted 28 June 2026; published online 10 August 2026; in print July 2026

Type: Research article · Language: English

License: CC BY-NC

Identifiers: DOI 10.1093/cercor/bhag113 · PMID 42574751 · PMCID PMC13456336 · OpenAlex W7202125338

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: EEG (modality), human (organism), methods / tools (subfield)

Methods: Spectral & time-frequency, Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning, Preprocessing, fMRI & imaging

Keywords: EEG, aperiodic, alpha, ageing, reliability

MeSH: Aging*, Brain*, Electroencephalography*, Rest*, Adult, Aged, Alpha Rhythm, Female, Follow-Up Studies, Humans, Male, Middle Aged, Reproducibility of Results, Young Adult (* major topic)

Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex

Funding: Bond University Postgraduate Research Scholarship; National Research Foundation Singapore; Campus for Research Excellence and Technological Enterprise

Citations: cited by 1 paper (Europe PMC); 65 references in the paper

Abstract

Several aspects of parameterized neural activity, including the aperiodic exponent and individual peak alpha frequency, have emerged as promising biomarkers for ageing, pathology, and cognitive decline. Their potential clinical application is tempered by a lack of evidence on long-term temporal stability. Existing investigations have largely relied on cross-sectional designs or considered stability for up to 90 days. Here, we examined five-year reliability, stability, and age-associated changes in periodic and aperiodic neural activity using electroencephalography in adults aged 20-70 years. Resting-state EEG was recorded in two sessions, approximately five years apart. We extracted the aperiodic exponent, aperiodic offset, peak alpha power, and individual alpha peak frequency and examined test-retest reliability at both the channel and cluster levels. All parameters demonstrated fair to excellent test-retest reliability (intraclass correlations = 0.51-0.88). Linear mixed models revealed that individual peak alpha frequency decreased, the aperiodic exponent flattened, and parameterized alpha power remained unchanged. There were no interactions between time and age. Our findings suggest that parameterized activity is reliable over long timeframes, and demonstrates changes consistent with ageing-related processes. Spectral parameterization may provide a means of characterizing within-person neurophysiological changes across adulthood. Future research should explore the utility of identifying deviations that may indicate pathology.

Reproduced under the paper's license (CC BY-NC), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 4 matches between paragraphs and lines of code.

MindSpaceLab/Aperiodic_Test_Retest_5year

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: e056acf5244f1b0becfc82591736031431f59ab3, 16 June 2026

Languages: Python (4), R (1), Jupyter (1)

Size: 56 files, 6 scripts

Software Heritage: not archived

Found in: “Data availability”

Holds: README, license file, environment (pyproject.toml), 2 notebooks

Not found: CITATION.cff, tests, continuous integration, documentation

Tools: MNE-Python (4 files), specparam (formerly FOOOF) (4 files), ICLabel (3 files), Matplotlib (3 files), NumPy (3 files), pandas (3 files), MNE-BIDS (2 files), PyPREP (2 files), scikit-learn (2 files), SciPy (2 files), broom (1 file), cowplot (1 file), easystats (1 file), emmeans (1 file), ggplot2 (1 file), ggpubr (1 file), lme4 (1 file), lmerTest (1 file), NeuroDSP (1 file), psych (1 file), seaborn (1 file), tidyverse (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

8 files

notebooks/IfAdo_analyses.Rmd, R, 2,077 lines
notebooks/IfAdo_preprocessing.ipynb, Jupyter, 1,986 lines, 2 matches
src/__init__.py, Python, 3 lines
src/preprocessing.py, Python, 312 lines, 2 matches
src/utils.py, Python, 1,086 lines
src/viz.py, Python, 47 lines
LICENSE, License, 21 lines
README.md, Text, 20 lines

At the source: github.com/MindSpaceLab/Aperiodic_Test_Retest_5year

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
6 scripts, each with its path and the digest of its content;
4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

doi:10.18112/openneuro.ds005385.v1.0.2, at OpenNeuro; found in “Data availability”

Data availability

All code used for all analyses and plots is publicly available at https://github.com/MindSpaceLab/Aperiodic_Test_Retest_5year. Raw data are available at https://doi.org/10.18112/openneuro.ds005385.v1.0.2.

Reproduced under the paper's license (CC BY-NC), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 8 authors, 5 keywords, 14 MeSH terms, 3 funders, 63 references.

Cite

This paper

Politanskaia, P., Bywater, J., Finley, A. J., Keage, H. A. D., Kelley, N. J., McKeown, D. J., Schinazi, V. R., & Angus, D. J. (2026). Long-term reliability and stability of parameterized resting state EEG: evidence from a five-year follow-up. Cerebral cortex (New York, N.Y. : 1991), 36(7), bhag113. https://doi.org/10.1093/cercor/bhag113

BibTeX

@article{politanskaia2026long, author = {Politanskaia, Polina and Bywater, Jacinta and Finley, Anna J and Keage, Hannah A D and Kelley, Nicholas J and McKeown, Daniel J and Schinazi, Victor R and Angus, Douglas J}, title = {{Long-term reliability and stability of parameterized resting state EEG: evidence from a five-year follow-up}}, journal = {Cerebral cortex (New York, N.Y. : 1991)}, year = {2026}, month = jul, volume = {36}, number = {7}, pages = {bhag113}, publisher = {Oxford University Press}, issn = {1047-3211}, doi = {10.1093/cercor/bhag113}, url = {https://doi.org/10.1093/cercor/bhag113}, pmid = {42574751}, pmcid = {PMC13456336} }

RIS

TY - JOUR AU - Politanskaia, Polina AU - Bywater, Jacinta AU - Finley, Anna J AU - Keage, Hannah A D AU - Kelley, Nicholas J AU - McKeown, Daniel J AU - Schinazi, Victor R AU - Angus, Douglas J TI - Long-term reliability and stability of parameterized resting state EEG: evidence from a five-year follow-up T2 - Cerebral cortex (New York, N.Y. : 1991) J2 - Cereb Cortex PY - 2026 DA - 2026/07/01 VL - 36 IS - 7 SP - bhag113 SN - 1047-3211 PB - Oxford University Press DO - 10.1093/cercor/bhag113 UR - https://doi.org/10.1093/cercor/bhag113 LA - en ER -

CSL-JSON

{ "id": "10.1093/cercor/bhag113", "type": "article-journal", "title": "Long-term reliability and stability of parameterized resting state EEG: evidence from a five-year follow-up", "container-title": "Cerebral cortex (New York, N.Y. : 1991)", "author": [ { "family": "Politanskaia", "given": "Polina" }, { "family": "Bywater", "given": "Jacinta" }, { "family": "Finley", "given": "Anna J" }, { "family": "Keage", "given": "Hannah A D" }, { "family": "Kelley", "given": "Nicholas J" }, { "family": "McKeown", "given": "Daniel J" }, { "family": "Schinazi", "given": "Victor R" }, { "family": "Angus", "given": "Douglas J" } ], "container-title-short": "Cereb Cortex", "volume": "36", "issue": "7", "page": "bhag113", "DOI": "10.1093/cercor/bhag113", "PMID": "42574751", "PMCID": "PMC13456336", "ISSN": "1047-3211", "publisher": "Oxford University Press", "URL": "https://doi.org/10.1093/cercor/bhag113", "language": "en", "issued": { "date-parts": [ [ 2026, 7, 1 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.