MachineLearning.py

# -*- coding: utf-8 -*-
"""
Updated Oct 18 2022

@author: Qianliang Li (glia@dtu.dk)

This script contains the code for the machine learning analysis

It should be run after Preprocessing.py and FeatureEstimation.py
"""

# Set working directory
import os
wkdir = "/home/glia/EEG/Final_scripts/"
os.chdir(wkdir)

# Load all libraries from the Preamble
from Preamble import *

# Load the questionnaire data
final_qdf = pd.read_csv("final_qdf.csv", sep=",", na_values=' ')

# Define cases as >= 44 total PCL
# Type: numpy array with subject id
cases = np.array(final_qdf["Subject_ID"][final_qdf["PCL_total"]>=44])
n_groups = 2
Groups = ["CTRL", "PTSD"]

# Make function to load EEG features
def load_features_df():
    # Load all features
    power_df = pd.read_pickle(Feature_savepath+"Power_df.pkl")
    fTBR_data_df = pd.read_pickle(Feature_savepath+"fTBR_df.pkl")
    asymmetry_df = pd.read_pickle(Feature_savepath+"asymmetry_df.pkl")
    PAF_data_df = pd.read_pickle(Feature_savepath+"PAF_data_FOOOF_df.pkl")
    PAF_data_df_global = pd.read_pickle(Feature_savepath+"PAF_data_FOOOF_global_df.pkl")
    OOF_data_df = pd.read_pickle(Feature_savepath+"OOF_data_FOOOF_df.pkl")
    con_data_df = pd.read_pickle(Feature_savepath+"con_data_source_drop_interpol_df.pkl")
    pec_data_df = pd.read_pickle(Feature_savepath+"pec_data_drop_interpol_ch_df.pkl")
    microstate_transition_data_df = pd.read_pickle(Feature_savepath+"microstate_transition_data_df.pkl")
    microstate_time_df = pd.read_pickle(Feature_savepath+"microstate_time_df.pkl")
    microstate_entropy_df = pd.read_pickle(Feature_savepath+"microstate_entropy_df.pkl")
    GC_data_df = pd.read_pickle(Feature_savepath+"GC_data_source_drop_interpol_df.pkl")
    H_data_df = pd.read_pickle(Feature_savepath+"H_data_df.pkl")
    H_data_df_global = pd.read_pickle(Feature_savepath+"H_data_global_df.pkl")
    
    # List of features
    EEG_features_name_list = [['Power'],
                              ['Frontal Theta Beta Ratio',
                              'Asymmetry'],
                              ['Peak Alpha Frequency',
                              'Global Peak Alpha Frequency'],
                              ["1/f exponent"],
                              ['imcoh'],
                              ['wpli'],
                              ['Power Envelope Correlation'],
                              ['Microstate Transition',
                              'Microstate Ratio Time',
                              'Microstate Entropy'],
                              ["Granger Causality"],
                              ['DFA Exponent',
                              'Global DFA Exponent']]
    
    # Arrange them to fit one 2D dataframe
    # Make function to add measurement column for indexing
    def add_measurement_column(df, measurement = "Text"):
        dummy_variable = [measurement]*df.shape[0]
        df.insert(1, "Measurement", dummy_variable)
        return df
    # Make function to convert column tuple to string
    def convertTupleHeader(header):
        header = list(header)
        str = "_".join(header)
        return str
    
    # Prepare overall dataframe
    EEG_features_df = pd.DataFrame(Subject_id, columns = ["Subject_ID"])
    
    # Add power spectral densities
    power_df = add_measurement_column(power_df, "Power")
    temp_df = power_df.pivot_table(index="Subject_ID",columns=["Measurement",
                                        "Quant_status", "Eye_status",

                                        "Freq_band", "Channel"], dropna=False,
                                   values="PSD").reset_index(drop=True)
    temp_df.columns = [convertTupleHeader(temp_df.columns[i]) for i in range(len(temp_df.columns))]
    
    EEG_features_df = pd.concat([EEG_features_df,temp_df], axis=1)
    
    # Add frontal theta/beta ratio
    # fTBR_data_df = add_measurement_column(fTBR_data_df, "Frontal Theta Beta Ratio")
    temp_df = fTBR_data_df.pivot_table(index="Subject_ID",columns=["Measurement",
                                        "Eye_status"], dropna=False,
                                       values="TBR").reset_index(drop=True)
    temp_df.columns = [convertTupleHeader(temp_df.columns[i]) for i in range(len(temp_df.columns))]
    
    EEG_features_df = pd.concat([EEG_features_df,temp_df], axis=1)
    
    # Add power asymmetry
    asymmetry_df = add_measurement_column(asymmetry_df, "Asymmetry")
    temp_df = asymmetry_df.pivot_table(index="Subject_ID",columns=["Measurement",
                                        "Eye_status", "Freq_band", "ROI"], dropna=False,
                                        values="Asymmetry_score").reset_index(drop=True)
    temp_df.columns = [convertTupleHeader(temp_df.columns[i]) for i in range(len(temp_df.columns))]
    
    EEG_features_df = pd.concat([EEG_features_df,temp_df], axis=1)
    
    # Add peak alpha frequency
    PAF_data_df = add_measurement_column(PAF_data_df, "Peak Alpha Frequency")
    temp_df = PAF_data_df.pivot_table(index="Subject_ID",columns=["Measurement",
                                        "Eye_status", "Channel"], dropna=False,
                                        values="Value").reset_index(drop=True)
    # NaN values are interpolated with means across channels for each condition
    eye_status = list(final_epochs[0].event_id.keys())
    for ee in eye_status:
        temp = temp_df.loc[:,("Peak Alpha Frequency",ee)] # get data
        temp = temp.T.fillna(temp.mean(axis=1)).T # fill (transpose used because fillna is axis=0)
        temp_df.loc[:,("Peak Alpha Frequency",ee)] = temp.to_numpy()
    # If there are still NaN the values are interpolated across channels and condition
    temp_df = temp_df.T.fillna(temp_df.mean(axis=1)).T
    
    temp_df.columns = [convertTupleHeader(temp_df.columns[i]) for i in range(len(temp_df.columns))]
    
    EEG_features_df = pd.concat([EEG_features_df,temp_df], axis=1)
    
    # Add 1/f exponent
    OOF_data_df = add_measurement_column(OOF_data_df, "1/f exponent")
    temp_df = OOF_data_df.pivot_table(index="Subject_ID",columns=["Measurement",
                                        "Eye_status", "Channel"], dropna=False,
                                        values="Value").reset_index(drop=True)
    # NaN values are interpolated with means across channels for each condition
    eye_status = list(final_epochs[0].event_id.keys())
    for ee in eye_status:
        temp = temp_df.loc[:,("1/f exponent",ee)] # get data
        temp = temp.T.fillna(temp.mean(axis=1)).T # fill (transpose used because fillna is axis=0)
        temp_df.loc[:,("1/f exponent",ee)] = temp.to_numpy()
    # If there are still NaN the values are interpolated across channels and condition
    temp_df = temp_df.T.fillna(temp_df.mean(axis=1)).T
    
    temp_df.columns = [convertTupleHeader(temp_df.columns[i]) for i in range(len(temp_df.columns))]
    
    EEG_features_df = pd.concat([EEG_features_df,temp_df], axis=1)
    
    # Add global peak alpha frequency
    #PAF_data_df_global = add_measurement_column(PAF_data_df_global, "Global_Peak_Alpha_Frequency") # already exists
    temp_df = PAF_data_df_global.pivot_table(index="Subject_ID",columns=["Measurement",
                                        "Eye_status"], dropna=False,
                                        values="Value").reset_index(drop=True)
    # NaN values are interpolated across eye condition
    temp_df = temp_df.T.fillna(temp_df.mean(axis=1)).T
    
    temp_df.columns = [convertTupleHeader(temp_df.columns[i]) for i in range(len(temp_df.columns))]
    
    EEG_features_df = pd.concat([EEG_features_df,temp_df], axis=1)
    
    # Add connectivity measurements
    #con_data_df = add_measurement_column(con_data_df, "Connectivity") # already exists
    temp_df = con_data_df.pivot_table(index="Subject_ID",columns=["Con_measurement",
                                        "Eye_status", "chx", "chy", "Freq_band"], dropna=True,
                                        values="Value").reset_index(drop=True)
    # Drop coh and plv, which are more susceptible to volume conduction
    temp_df = temp_df.drop("coh",axis=1)
    temp_df = temp_df.drop("plv",axis=1)
    
    temp_df.columns = [convertTupleHeader(temp_df.columns[i]) for i in range(len(temp_df.columns))]
    
    EEG_features_df = pd.concat([EEG_features_df,temp_df], axis=1)
    
    # Add orthogonalized power enveloped correlations
    pec_data_df = add_measurement_column(pec_data_df, "Power Envelope Correlation")
    temp_df = pec_data_df.pivot_table(index="Subject_ID",columns=["Measurement",
                                        "Eye_status", "chx", "chy", "Freq_band"], dropna=True,
                                        values="Value").reset_index(drop=True)
    temp_df.columns = [convertTupleHeader(temp_df.columns[i]) for i in range(len(temp_df.columns))]
    
    EEG_features_df = pd.concat([EEG_features_df,temp_df], axis=1)
    
    # Add microstate transition probabilities
    microstate_transition_data_df = add_measurement_column(microstate_transition_data_df, "Microstate Transition")
    temp_df = microstate_transition_data_df.pivot_table(index="Subject_ID",columns=["Measurement",
                                        "Eye_status", "Transition"], dropna=False,
                                        values="Value").reset_index(drop=True)
    temp_df.columns = [convertTupleHeader(temp_df.columns[i]) for i in range(len(temp_df.columns))]
    
    EEG_features_df = pd.concat([EEG_features_df,temp_df], axis=1)
    
    # Add microstate time covered
    microstate_time_df = add_measurement_column(microstate_time_df, "Microstate Ratio Time")
    # Convert microstate to str before using it as column
    microstate_time_df = microstate_time_df.astype({"Microstate": str})
    
    temp_df = microstate_time_df.pivot_table(index="Subject_ID",columns=["Measurement",
                                        "Eye_status", "Microstate"], dropna=False,
                                        values="Value").reset_index(drop=True)
    temp_df.columns = [convertTupleHeader(temp_df.columns[i]) for i in range(len(temp_df.columns))]
    
    EEG_features_df = pd.concat([EEG_features_df,temp_df], axis=1)
    
    # Add microstate entropy
    #microstate_entropy_df = add_measurement_column(microstate_entropy_df, "Microstate_Entropy") # already exists