<jupyter_start><jupyter_text>Group Analysis - FOOOFed EEG Analysis: TaskExamine the results of parameterizing spectra for task EEG data, and predicting behavioural data. Throughout these analyses, we compare the `YNG` and `OLD` subject groups.This notebook does the analyses and creates the plots for Figure 6.<jupyter_code>%matplotlib inline
%config InlineBackend.figure_format='retina'
import warnings
from copy import deepcopy
from os.path import join as pjoin

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import spearmanr, ttest_1samp, normaltest

from statsmodels.stats.anova import anova_lm
from statsmodels.stats.diagnostic import compare_cox, compare_j

from fooof.plts import plot_spectra
from fooof.sim.gen import gen_aperiodic, gen_power_spectrum
# Import custom code for this analysis
import sys
sys.path.append('../code')

from plts import *
from utils import *
from data_mgmt import *
from analysis import *<jupyter_output><empty_output><jupyter_text>Settings<jupyter_code># Import general settings
from settings import YNG_INDS, OLD_INDS
from settings import RESULTS_PATH
# Set folder for FOOOF results
folder = 'FOOOF'

# Model settings - which group to run
group = 'ALL' # 'ALL', 'YNG', 'OLD'

print('Number of YNG subjects: ', len(YNG_INDS))
print('Number of OLD subjects: ', len(OLD_INDS))
# Data settings
srate = 512
tmin, tmax = -0.85, 1.1
times = np.arange(tmin, tmax, 1/srate)
seg_times = [(-0.85, -0.35), (0.1, 0.6), (0.5, 1.0)]
n_subjs = 31
dict_structure = {'YNG' : {}, 'OLD' : {}, 'ALL' : {}}
load_side = 'Contra'

# Set the time segments indices to compute the difference measures with
i1, i2 = 2, 0
# Wether to save out plots or not
SAVE_FIG = False<jupyter_output><empty_output><jupyter_text>Check dropped trials<jupyter_code># Load dropped trials & components
dropped_trials = np.load(pjoin(RESULTS_PATH, 'Group', 'dropped_trials.npy'))
dropped_components = np.load(pjoin(RESULTS_PATH, 'Group', 'dropped_components.npy'))

# Check dropped trials for each subject
print('SubNum \t\t # Dropped Trials \t # Dropped Components')
for ind, trials, components in zip(range(n_subjs), dropped_trials, dropped_components):
    temp_trials = trials[trials < 999.]
    temp_comps = components[components < 999.]
    #print(ind, '\t\t', len(temp_trials), '\t\t\t', len(temp_comps))<jupyter_output>SubNum 		 # Dropped Trials 	 # Dropped Components<jupyter_text>Group FOOOFing - Trial Averaged DataData objects are 3d arrays, with the shape `[n_loads, n_subjs, n_times]`.  Load Data<jupyter_code># Load behavioural data
behav_data = pd.read_csv(pjoin(RESULTS_PATH, 'Behav', 'neural_aging_data_behaviour.csv'))
# Convert data types
behav_data['Age'] = behav_data['Age'].astype('str')
behav_data['Load'] = behav_data['Load'].astype('str')
# Calculate average behaviour across loads
avg_behav = behav_data.groupby('SubjID').mean()
# Load and extract FOOOF data
all_offsets, all_exps = load_fooof_task_ap(RESULTS_PATH, load_side, folder)
all_alphas_cf = load_fooof_task_pe(RESULTS_PATH, load_side, 0, folder)
all_alphas_pw = load_fooof_task_pe(RESULTS_PATH, load_side, 1, folder)
all_alphas_bw = load_fooof_task_pe(RESULTS_PATH, load_side, 2, folder)
all_r2s, all_errs = load_fooof_task_md(RESULTS_PATH, load_side, folder)
# Load canonical alpha analysis
canalph_group = np.load(pjoin(RESULTS_PATH, 'Group', 'canonical_group.npy'))

# Load individualized frequency canonical alpha analysis
canalph_icf_group = np.load(pjoin(RESULTS_PATH, 'Group', 'canonical_icf_group.npy'))
# Average across analytic alpha measures to get canonical alpha measure
seg_masks = []
for seg in seg_times:
    seg_masks.append(np.logical_and(times >= seg[0], times <= seg[1]))

canalpha = np.zeros_like(all_alphas_pw)
canal_icf = np.zeros_like(all_alphas_pw)
for subi, (canalph_data, canalph_icf_data) in enumerate(zip(canalph_group, canalph_icf_group)):
    for lodi in range(3):
        for segi, mask in enumerate(seg_masks):
            canalpha[lodi, subi, segi] = np.mean(canalph_data[lodi, mask])
            canal_icf[lodi, subi, segi] = np.mean(canalph_icf_data[lodi, mask])<jupyter_output><empty_output><jupyter_text>Data Checks<jupyter_code># Check number of missing FOOOFed alphas
print('Number of missing FOOOFed alphas: \t', np.sum(np.isnan(all_alphas_pw)))
# Check if there are any NaN values
nans = np.isnan(all_alphas_pw)
print('Total number of NaN values:\t\t', np.sum(nans))<jupyter_output>Total number of NaN values:		 0<jupyter_text>Data Management<jupyter_code># Make a data dictionary - each with shape [n_conds, n_times]
data_dict = deepcopy(dict_structure)
diff_data_dict = deepcopy(dict_structure)
behav_dict = deepcopy(dict_structure)
# Collect the data and labels into lists for checking through
labels = ['offset', 'exponent',
          'alpha_cf', 'alpha_pw', 'alpha_bw',
          'canalpha', 'canal_icf']
datas = [all_offsets, all_exps,
         all_alphas_cf, all_alphas_pw, all_alphas_bw,
         canalpha, canal_icf]
# Set up data & diff_data dicts
for label, data in zip(labels, datas):
    data_dict['YNG'][label], data_dict['OLD'][label] = reshape_data(data)
    data_dict['ALL'][label] = np.concatenate([data_dict['YNG'][label],
                                              data_dict['OLD'][label]])
    
    diff_data_dict['YNG'][label] = calc_diff(data_dict['YNG'][label], i1, i2)
    diff_data_dict['OLD'][label] = calc_diff(data_dict['OLD'][label], i1, i2)
    diff_data_dict['ALL'][label] = np.concatenate([diff_data_dict['YNG'][label],
                                                   diff_data_dict['OLD'][label]])
# Set up the behavioural data dict
for label in ["d'", "Load", 'CDA']:
    behav_dict['ALL'][label] = behav_data[label].values
    behav_dict['YNG'][label] = behav_data[behav_data['Age'] == '1'][label].values
    behav_dict['OLD'][label] = behav_data[behav_data['Age'] == '2'][label].values<jupyter_output><empty_output><jupyter_text>Data CheckingCheck the FOOOF fit measures, checking whether there is a systematic difference between the groups that could explain the differences.<jupyter_code># Check if the fit metrics vary by group
for label, data in zip(['Error', 'R^2'], [all_errs, all_r2s]):
    print(label)
    m1, m2 = np.nanmean(data[:, YNG_INDS, :]), np.nanmean(data[:, OLD_INDS, :])
    print('\tMeans: \t\t  {:1.4f} \t {:1.4f} \t diff: {:1.4f}'.format(m1, m2, np.abs(m1-m2)))
    print_stat('\tt-test', *ttest_ind(data[:, YNG_INDS, :].flatten(),
                                      data[:, OLD_INDS, :].flatten()))
# Check the correlation between FOOOF fit error & R^2
print_stat('Err-R^2 Corr', *spearmanr(all_errs.flatten(), all_r2s.flatten()))
# Check if the fit metrics correlate with measured FOOOF parameters
print('\t\t  Error \t   R^2')
for label, data in zip(labels, datas):
    print_stat(label,
               spearmanr(all_errs.flatten(), data.flatten())[0],
               spearmanr(all_r2s.flatten(), data.flatten())[0])<jupyter_output>Error 	   R^2
offset: 	  0.0828 	 0.6035
exponent: 	 -0.1176 	 0.6869
alpha_cf: 	  0.1621 	-0.0953
alpha_pw: 	  0.7591 	 0.1712
alpha_bw: 	  0.2681 	 0.2580
canalpha: 	  0.2465 	-0.0506
canal_icf: 	  0.2475 	-0.0955<jupyter_text>Descriptive Data Checks<jupyter_code># Print out mean values, per group
print('Average Values per measure, per group: \n')
print('\t\t   YNG \t\t  OLD')
for label in labels:
    print_stat(label,
               np.nanmean(diff_data_dict['YNG'][label]),
               np.nanmean(diff_data_dict['OLD'][label]))
# Print out standard deviations, per group
print('Standard deviation per measure, per group: \n')
print('\t\t   YNG \t\t  OLD')
for label in labels:
    print_stat(label,
               np.nanstd(diff_data_dict['YNG'][label]),
               np.nanstd(diff_data_dict['OLD'][label]))
# Check for normality, per measure, per group
print('Normal test per measure, per group: \n')
print('YNG\t\t   stat \t  p-val')
for label in labels:
    print_stat(label, *normaltest(diff_data_dict['YNG'][label]))
print('\nOLD\t\t   stat \t  p-val')
for label in labels:
    print_stat(label, *normaltest(diff_data_dict['OLD'][label]))
# Print out tests for group differences
print('T-Tests for Differences Between Groups (YNG vs OLD), per Metric: \n')
print('\t\t   t-val \t  p-val')
for label in labels:
    print_stat(label, *nan_ttest(diff_data_dict['YNG'][label],
                                 diff_data_dict['OLD'][label]))
# Check correlations of physiology to behaviour, per group
print('Correlation (r-vals) of each metric with behaviour, per Group: \n')
print('\t\t   YNG \t\t  OLD')
for label in labels:
    print_stat(label,
    nan_corr(diff_data_dict['YNG'][label], behav_dict['YNG']["d'"])[0],
    nan_corr(diff_data_dict['OLD'][label], behav_dict['OLD']["d'"])[0])<jupyter_output>Correlation (r-vals) of each metric with behaviour, per Group: 

		   YNG 		  OLD
offset: 	  0.1435 	 0.0341
exponent: 	  0.2729 	-0.1506
alpha_cf: 	  0.0938 	-0.1061
alpha_pw: 	 -0.1167 	 0.3955
alpha_bw: 	 -0.1312 	 0.1600
canalpha: 	 -0.2216 	 0.2717
canal_icf: 	 -0.2401 	 0.3341<jupyter_text>Fit Models to Predict BehaviourPredict behaviour output from evoked responses of alpha and aperiodic parameters.<jupyter_code># Set up model definitions
models_defs = {
    'base_model' : 'behav ~ load',
    'cda_model' : 'behav ~ load + cda',
    'canal_model' : 'behav ~ load + al_pw_dif',
    'canicf_model' : 'behav ~ load + al_pw_dif',
    'f_al_model' : 'behav ~ load  + al_pw_dif',
    'f_al_p_model' : 'behav ~ load + al_cf_dif + al_pw_dif + al_bw_dif',
    'f_ap_model' : 'behav ~ load + off_dif + exp_dif',
    'f_ap_p_model' : 'behav ~ load + off_dif + exp_dif + al_pw_dif', 
}
# Set up model data
models_data = {
    'base_model' : lambda group : {'behav' : behav_dict[group]["d'"],
                                   'load' : behav_dict[group]['Load']},
    
    'cda_model' : lambda group : {'behav' : behav_dict[group]["d'"],
                                  'load' : behav_dict[group]['Load'],
                                  'cda' : behav_dict[group]['CDA']}, 
    
    'canal_model' : lambda group : {'behav' : behav_dict[group]["d'"],
                                    'load' : behav_dict[group]['Load'],
                                    'al_pw_dif' : diff_data_dict[group]['canalpha']},
    
    'canicf_model' : lambda group : {'behav' : behav_dict[group]["d'"],
                                        'load' : behav_dict[group]['Load'],
                                        'al_pw_dif' : diff_data_dict[group]['canal_icf']},
    
    'f_al_model' : lambda group : {'behav' : behav_dict[group]["d'"],
                                   'load' : behav_dict[group]['Load'],
                                   'al_pw_dif' : diff_data_dict[group]['alpha_pw']},
    
    'f_al_p_model' : lambda group : {'behav' : behav_dict[group]["d'"],
                                     'load' : behav_dict[group]['Load'],
                                     'al_cf_dif' : diff_data_dict[group]['alpha_cf'],
                                     'al_pw_dif' : diff_data_dict[group]['alpha_pw'],
                                     'al_bw_dif' : diff_data_dict[group]['alpha_bw']},
    
    'f_ap_model' : lambda group :  {'behav' : behav_dict[group]["d'"],
                                    'load' : behav_dict[group]['Load'],
                                    'off_dif' : diff_data_dict[group]['offset'],
                                    'exp_dif' : diff_data_dict[group]['exponent']},
    
    'f_ap_p_model' : lambda group :  {'behav' : behav_dict[group]["d'"],
                                      'load' : behav_dict[group]['Load'],
                                      'al_pw_dif' : diff_data_dict[group]['alpha_pw'],
                                      'off_dif' : diff_data_dict[group]['offset'],
                                      'exp_dif' : diff_data_dict[group]['exponent']},

}
# Fit models, across groups
models_results = deepcopy(dict_structure)
for group in models_results.keys():
    for model_name, model_def in models_defs.items():
        models_results[group][model_name] = run_model(model_def,
                                                      models_data[model_name](group),
                                                      print_model=False)
# Check the model performances - per model def per group
print('\t\t\t  R^2 Adj \t  p-val')
for group_label, group_models in models_results.items():
    print(group_label)
    for model_label, model_results in group_models.items():
        print_stat('\t' + model_label,
                   model_results.rsquared_adj,
                   model_results.f_pvalue)
# Check a particular model
group = 'OLD'
model = 'f_ap_model'
print(models_results[group][model].summary())<jupyter_output>OLS Regression Results                            
==============================================================================
Dep. Variable:                  behav   R-squared:                       0.355
Model:                            OLS   Adj. R-squared:                  0.285
Method:                 Least Squares   F-statistic:                     5.095
Date:                Sat, 17 Oct 2020   Prob (F-statistic):            0.00226
Time:                        01:57:12   Log-Likelihood:                -50.625
No. Observations:                  42   AIC:                             111.3
Df Residuals:                      37   BIC:                             119.9
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|[...]<jupyter_text>Model ComparisonsExplicitly test for differences between different model fits.  Comparing Nested Models`statsmodels` offers three tests for nested models: f test, lagrange multiplier, likelihood ratio. Note that these three can be called from a results object, as `compare_x_test` with `f`, `lm` and `lr` as `x`. F-test can also be run with `anova_lm`.  Comparing Non-Nested ModelsStatmodels offers two tests for non-nested model: cox test & j testThey are better described in the R implementations:- cox_test: http://math.furman.edu/~dcs/courses/math47/R/library/lmtest/html/coxtest.html- j_test: http://math.furman.edu/~dcs/courses/math47/R/library/lmtest/html/jtest.html<jupyter_code># Compare nested models: alpha models vs base models
for group in ['OLD', 'YNG']:
    print('\nGroup: ', group, '\n')
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        print('Canonical alpha vs. Base Model')
        print(anova_lm(models_results[group]['base_model'],
                       models_results[group]['canal_model']))
        print('\n')
        print('FOOOF alpha vs. Base Model')
        print(anova_lm(models_results[group]['base_model'],
                       models_results[group]['f_al_model']))
        print('\n')
        print('FOOOF aperiodic vs. Base Model')
        print(anova_lm(models_results[group]['base_model'],
                       models_results[group]['f_ap_model']))
# Compare different alpha models
print('Canonical alpha vs. FOOOFed Alpha')
for group in ['OLD', 'YNG']:
    print('Group: ', group)
    print_stat('\tAlpha-Model Compare',
               *compare_cox(models_results[group]['canal_model'],
                            models_results[group]['f_al_model']))
# Compare if the + models are any better
for group in ['OLD', 'YNG']:
    print('\nGroup: ', group, '\n')
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        print('FOOOF Alpha + vs. FOOOF Alpha')
        print(anova_lm(models_results[group]['f_al_model'],
                       models_results[group]['f_al_p_model']))
        print('\n')
        print('FOOOF AP + vs. FOOOF AP')
        print(anova_lm(models_results[group]['f_ap_model'],
                       models_results[group]['f_ap_p_model']))<jupyter_output>Group:  OLD 

FOOOF Alpha + vs. FOOOF Alpha
   df_resid        ssr  df_diff   ss_diff         F    Pr(>F)
0      38.0  26.429804      0.0       NaN       NaN       NaN
1      36.0  26.202530      2.0  0.227274  0.156127  0.856027


FOOOF AP + vs. FOOOF AP
   df_resid        ssr  df_diff  ss_diff         F   Pr(>F)
0      37.0  27.399601      0.0      NaN       NaN      NaN
1      36.0  24.018710      1.0  3.38089  5.067385  0.03058

Group:  YNG 

FOOOF Alpha + vs. FOOOF Alpha
   df_resid        ssr  df_diff   ss_diff         F    Pr(>F)
0      47.0  34.338756      0.0       NaN       NaN       NaN
1      45.0  32.655867      2.0  1.682889  1.159516  0.322831


FOOOF AP + vs. FOOOF AP
   df_resid        ssr  df_diff   ss_diff         F    Pr(>F)
0      46.0  28.706282      0.0       NaN       NaN       NaN
1      45.0  28.637508      1.0  0.068773  0.108068  0.743879<jupyter_text>ConclusionsIn general, we can see that predicting behaviour from parameterized spectral features appears to be promising.Key points:- In the cases in which alpha is predictive the parameterized alpha tends to do better- There is also some predictive power of the aperiodic components  Check Which Parameters ChangeNext we will check which parameters show a significant task related difference.In the cell below, per group, per measure, the average difference of the measure (from during the trial vs baseline) is printed, as well as a 1 sample t-test, which tests if this number is significantly different from 0 (no change).<jupyter_code>for age in ['YNG', 'OLD']:
    print(age)
    for label in ['offset', 'exponent', 'alpha_cf', 'alpha_pw', 'alpha_bw']:
        temp = diff_data_dict[age][label]
        print('\n', label)
        print('\t avg diff: \t {:+1.3f}'.format(np.mean(temp)))
        print('\t effect size: \t {:+1.3f}'.format(cohens_d_1samp(temp)))
        print_stat('\t t-test', *ttest_1samp(temp, 0))<jupyter_output>YNG

 offset
	 avg diff: 	 -0.057
	 effect size: 	 -0.533
	 t-test: 	 -3.8094 	 0.0004

 exponent
	 avg diff: 	 -0.035
	 effect size: 	 -0.353
	 t-test: 	 -2.5212 	 0.0149

 alpha_cf
	 avg diff: 	 +0.029
	 effect size: 	 +0.061
	 t-test: 	  0.4381 	 0.6632

 alpha_pw
	 avg diff: 	 +0.099
	 effect size: 	 +0.750
	 t-test: 	  5.3537 	 0.0000

 alpha_bw
	 avg diff: 	 -0.081
	 effect size: 	 -0.232
	 t-test: 	 -1.6569 	 0.1038
OLD

 offset
	 avg diff: 	 -0.051
	 effect size: 	 -0.444
	 t-test: 	 -2.8806 	 0.0063

 exponent
	 avg diff: 	 -0.007
	 effect size: 	 -0.059
	 t-test: 	 -0.3799 	 0.7060

 alpha_cf
	 avg diff: 	 +0.300
	 effect size: 	 +0.350
	 t-test: 	  2.2651 	 0.0289

 alpha_pw
	 avg diff: 	 -0.002
	 effect size: 	 -0.019
	 t-test: 	 -0.1250 	 0.9011

 alpha_bw
	 avg diff: 	 +0.082
	 effect size: 	 +0.131
	 t-test: 	  0.8492 	 0.4007<jupyter_text>Plot the difference of FOOOFs Now we recreate power spectra that visualize these changes, per group.To do so, we will use the values shown above to reconstruct how spectra change from baseline to task engaged.<jupyter_code># Generate spectra that show the task related change in spectral activity
for group, inds in zip(['YNG', 'OLD'], [YNG_INDS, OLD_INDS]):

    time_offs = all_offsets[:, inds, :].mean(0).mean(0)
    time_exps = all_exps[:, inds, :].mean(0).mean(0)
    time_al_cf = all_alphas_cf[:, inds, :].mean(0).mean(0)
    time_al_pw = all_alphas_pw[:, inds, :].mean(0).mean(0)
    time_al_bw = all_alphas_bw[:, inds, :].mean(0).mean(0)

    fs, base_spectrum = gen_power_spectrum([3, 30],
                                           [time_offs[0], time_exps[0]],
                                           [time_al_cf[0], time_al_pw[0],
                                            time_al_bw[0]/2], nlv=0)
    fs, task_spectrum = gen_power_spectrum([3, 30],
                                           [time_offs[2], time_exps[2]],
                                           [time_al_cf[2], time_al_pw[2],
                                            time_al_bw[2]/2], nlv=0)

    plot_spectra(fs, [base_spectrum, task_spectrum], True, True,
                 labels=['Baseline', 'Delay'], figsize=(7, 6), lw=3.5)
    plt.title(group, {'fontsize' : 16, 'fontweight' : 'bold'})

    if SAVE_FIG:
        plt.tight_layout()
        plt.savefig('../figures/' + 'EEG-TaskModels-' + group + '.pdf')<jupyter_output><empty_output><jupyter_text>Correlations Between FeaturesCheck the correlation structure within and between FOOOF and canonical measures.<jupyter_code># Settings
group = 'ALL'
cur_data = data_dict
#cur_data = diff_data_dict
# Calculate and collect correlations between features
corrs = np.zeros(shape=[len(datas), len(datas)])
print('Correlations between features \t\t r-val \t\t p-val')
for i1, label1 in enumerate(labels):
    for i2, label2 in enumerate(labels):
        corr = nan_corr(cur_data[group][label1].flatten(),
                        cur_data[group][label2].flatten())
        if i1 != i2:
            corrs[i1, i2] = corr[0]
        else:
            corrs[i1, i2] = 0
        if i1 < i2:
            print('  {:10s} &   {} : \t\t{:+1.3f} \t\t{:1.4f}'.format(label1, label2, *corr))
# Plot the correlation matrix between features
plt.imshow(corrs, cmap='bwr', vmin=-1, vmax=1)
plt.gca().set_xticklabels([''] + labels, rotation=45);
plt.gca().set_yticklabels([''] + labels, rotation=45);
plt.colorbar();

if SAVE_FIG:
    plt.tight_layout()
    plt.savefig('../figures/EEG-ParamCorrs.pdf')<jupyter_output><empty_output><jupyter_text>Explore Relation Between Canonical Alpha and Spectral ParametersIn this section, we will explore how the canonical alpha measure relates to FOOOF measures.The idea of the FOOOF model is that the measured power at a particular frequency, say 10 Hz or alpha, reflects a combination of aperiodic power and any band-specific oscillations. If so, the total power at a given frequency, as is traditionally computed, should be predictable by a combination of FOOOF features.To check this, we will further explore the relationship between canonically measured alpha power, and FOOOF features for the aperiodic component (offset and exponent) and the FOOOF measure of alpha peak-specific power.<jupyter_code># Calculate correlation between canonical and FOOOF alpha
print_stat('C-Alpha vs F-Alpha     ', *nan_corr(np.array(all_alphas_pw.flatten()),
                                                np.array(canalpha.flatten())))
print_stat('C-ICF-Alpha vs F-Alpha ', *nan_corr(np.array(all_alphas_pw.flatten()),
                                                np.array(canal_icf.flatten())))
# Generate the power @ alpha frequency given the aperiodic component, from the FOOOF fits
ap_alpha = []
for cf, off, exp in zip(all_alphas_cf.flatten(), all_offsets.flatten(), all_exps.flatten()):
    ap_alpha.append(gen_aperiodic(np.array([10]), [off, exp])[0])
ap_alpha = np.array(ap_alpha)

# Calculate the total power at 10 Hz (or about) from the combination of aperiodic & alpha peak
foo_total = ap_alpha + all_alphas_pw.flatten()
# Calculate correlation between canonical alpha and aperiodic component @ 10 Hz
print_stat('C-Alpha vs 10Hz-AP', *nan_corr(np.array(ap_alpha),
                                           np.array(canalpha.flatten())))
# Calculate correlation between the canonical alpha and the FOOOF model total @ 10 Hz
print_stat('C-Alpha vs Total-FOOOF', *nan_corr(np.array(foo_total.flatten()),
                                               np.array(canalpha.flatten())))
print_stat('C-ICF-Alpha vs Total-FOOOF', *nan_corr(np.array(foo_total.flatten()),
                                                   np.array(canal_icf.flatten())))<jupyter_output>C-Alpha vs Total-FOOOF: 	  0.2628 	 0.0000
C-ICF-Alpha vs Total-FOOOF: 	  0.2545 	 0.0000<jupyter_text>To put it all together, let's see to what extent we can predict the canonically measured alpha from FOOOF features.<jupyter_code># Predict canonical alpha
#   Note: if this is updated to predict ICF instead, results are approximately the same
data = pd.DataFrame()

group = 'ALL'
data['exp'] = diff_data_dict[group]['exponent']
data['off'] = diff_data_dict[group]['offset']
data['fal'] = diff_data_dict[group]['alpha_pw']
data['canal'] = diff_data_dict[group]['canalpha']

mod = run_model("canal ~ exp + off + fal", data, True)
# Check the predictions from individual components
m1 = run_model("canal ~ exp", data, False)
m2 = run_model("canal ~ off", data, False)
m3 = run_model("canal ~ fal", data, False)

print('Canal - EXP: \t{:1.4f}'.format(m1.rsquared_adj))
print('Canal - OFF: \t{:1.4f}'.format(m2.rsquared_adj))
print('Canal - FAL: \t{:1.4f}'.format(m3.rsquared_adj))<jupyter_output>Canal - EXP: 	0.1709
Canal - OFF: 	0.0290
Canal - FAL: 	0.8006
StarcoderdataJupyterStructuredCleanDedup