# -*- 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
microstate_entropy_df["Measurement"] = ["Microstate Entropy"]*microstate_entropy_df.shape[0]
temp_df = microstate_entropy_df.pivot_table(index="Subject_ID",columns=["Measurement",
"Eye_status"], 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 Granger Causality # fixed so it is always chx --> chy
GC_data_df = add_measurement_column(GC_data_df, "Granger Causality")
temp_df = GC_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))]
# Check if any GC was not calculated properly
assert temp_df.shape[1] == (31*31-31)*5*2 # expected number of features with 31 sources
EEG_features_df = pd.concat([EEG_features_df,temp_df], axis=1)
# Add DFA exponent
H_data_df = add_measurement_column(H_data_df, "DFA Exponent")
temp_df = H_data_df.pivot_table(index="Subject_ID",columns=["Measurement",
"Eye_status", "Channel", "Freq_band"], 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 Global Hurst exponent
H_data_df_global = H_data_df_global.drop("Measurement",axis=1)
H_data_df_global = add_measurement_column(H_data_df_global, "Global DFA Exponent") # already exists
temp_df = H_data_df_global.pivot_table(index="Subject_ID",columns=["Measurement",
"Eye_status", "Freq_band"], 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)
return EEG_features_df, EEG_features_name_list
# %% ML predction of PTSD status using all features
# mRMR -> mRMR -> SVM/RF/LogReg with L2
# 10 repetitions of 10 fold outer and 10 fold inner two-layer CV
EEG_features_df, EEG_features_name_list = load_features_df()
# Add group status
Group_status = np.array([0]*EEG_features_df.shape[0]) # CTRL = 0
Group_status[np.array([i in cases for i in EEG_features_df["Subject_ID"]])] = 1 # PTSD = 1
EEG_features_df.insert(1, "Group_status", Group_status)
# Subject info columns
Subject_info_cols = list(EEG_features_df.columns[0:2])
n_subject_info_cols = len(Subject_info_cols)
n_discrete_cols = 2
# To ensure proper stratification into train/test set I will stratify using group status and study status
# A variable that encodes for this is created
n_studies = 3
study_group_status = EEG_features_df["Group_status"].copy()
for i in range(n_studies):
# Get study index
study_idx = (EEG_features_df["Subject_ID"]>=(i+1)*100000)&(EEG_features_df["Subject_ID"]<(i+2)*100000)
# Assign label
study_group_status[(study_idx)&(EEG_features_df["Group_status"]==0)] = 2*i # CTRL
study_group_status[(study_idx)&(EEG_features_df["Group_status"]==1)] = 2*i+1 # PTSD
# Target variable
Target = ["Group_status"]
Target_col = EEG_features_df.iloc[:,1:].columns.isin(Target)
Target_col_idx = np.where(Target_col == True)[0]
# Make 3 models and save them to use enseemble in the end
CLF_models = ["SVM", "LogReg", "RF"]
n_models = len(CLF_models)
# Repeat the classification to see if I just got a lucky seed
n_repetitions = 10
k_out = 10
accuracy_arr = np.zeros((n_repetitions,k_out,n_models,3))
model_par = []
final_models = []
final_features = []
final_y_preds = []
np.random.seed(42)
# Prepare the splits beforehand to make sure the repetitions are not the same
Rep_Outer_CV = []
Rep_Outer_CV_test = [] # a list with only the test indices
for rep in range(n_repetitions):
skf = StratifiedKFold(n_splits=k_out, shuffle=True, random_state=(rep+1)*123) # using 10% equals around 21 test subjects
# I am also using converting it to an iterable list instad of the generator to promote reuse
Outer_CV = []
Outer_CV_test = []
for train_index, test_index in skf.split(EEG_features_df,study_group_status):
Outer_CV.append([train_index,test_index])
Outer_CV_test.append(test_index)
Rep_Outer_CV.append(Outer_CV)
Rep_Outer_CV_test.append(Outer_CV_test)
# Check none of the repetitions are the same by using only the test sets
# The list is first flattened
Rep_Outer_CV_test_flat = [item for sublist in Rep_Outer_CV_test for item in sublist]
# All elements are converted to strings
# This makes it easier to look for uniques, and the indices are already in numerical order
Rep_Outer_CV_test_flat_str = ["".join([str(x) for x in ele])for ele in Rep_Outer_CV_test_flat]
def allUnique(x):
seen = set()
return not any(i in seen or seen.add(i) for i in x)
assert allUnique(Rep_Outer_CV_test_flat_str)
# Get current time
c_time1 = time.localtime()
c_time1 = time.strftime("%a %d %b %Y %H:%M:%S", c_time1)
print(c_time1)
for rep in range(n_repetitions):
# The outer fold CV has already been saved as lists
Outer_CV = Rep_Outer_CV[rep]
# Pre-allocate memory
model_par0 = []
final_models0 = []
final_features0 = []
final_y_preds0 = []
Outer_counter = 0
for train_index, test_index in Outer_CV:
test_df = EEG_features_df.iloc[test_index]
train_df = EEG_features_df.iloc[train_index]
# Training data will be standardized
standardizer = preprocessing.StandardScaler().fit(train_df.iloc[:,n_discrete_cols:])
train_df_standard = train_df.copy()
train_df_standard.iloc[:,n_discrete_cols:] = standardizer.transform(train_df_standard.iloc[:,n_discrete_cols:])
# Test data will also be standardized but using mean and std from training data
test_df_standard = test_df.copy()
test_df_standard.iloc[:,n_discrete_cols:] = standardizer.transform(test_df_standard.iloc[:,n_discrete_cols:])
# Get the training data
X_train = train_df_standard.copy().drop(Subject_info_cols, axis=1)
y_train = train_df_standard[Target]
# Get test data
X_test = test_df_standard.copy().drop(Subject_info_cols, axis=1)
y_test = test_df_standard[Target]
# Prepare initial filtering of feature types to alleviate imbalance in number of features
# Use a variable that is optimized using inner CV
n_feat_mRMR1 = [20,30,40,50]
n_feat_mRMR2 = [30,40,50,60]
# Max features from mRMR for each eye status
max_mRMR_features = n_feat_mRMR1[-1]
max_mRMR_features2 = n_feat_mRMR2[-1]
eye_status = ["Eyes Closed", "Eyes Open"]
n_eye_status = len(eye_status)
mRMR_features = []
feat_counter = []
# Perform mRMR for each feature type
for fset in range(len(EEG_features_name_list)):
temp_feat = EEG_features_name_list[fset]
other_feats = EEG_features_name_list[:fset]+EEG_features_name_list[fset+1:]
other_feats = [item for sublist in other_feats for item in sublist] # make the list flatten
for e in range(n_eye_status):
ee = eye_status[e]
temp_ee = "{}_{}".format(temp_feat,ee)
# Retrieve the dataset for each feature
col_idx = np.zeros(len(X_train.columns), dtype=bool)
for fsub in range(len(temp_feat)):
temp_feat0 = temp_feat[fsub]
col_idx0 = X_train.columns.str.contains(temp_feat0+"_")
col_idx = np.logical_or(col_idx,col_idx0) # append all trues
temp_X_train = X_train.loc[:,col_idx]
# Check if any of the other features were wrongly chosen
for fcheck in range(len(other_feats)):
if any(temp_X_train.columns.str.contains(other_feats[fcheck]+"_")==True):
temp_X_train = temp_X_train.loc[:,np.invert(temp_X_train.columns.str.contains(other_feats[fcheck]+"_"))]
if temp_X_train.size == 0: # if no columns are left, e.g. imcoh when removing coh, then add features again
temp_X_train = X_train.loc[:,X_train.columns.str.contains(temp_feat0+"_")]
# Index for eye status
temp_X_train = temp_X_train.loc[:,temp_X_train.columns.str.contains(ee)]
# Save number of original features fed into mRMR
feat_counter.append([temp_feat,temp_X_train.shape[1]])
# Do not use mRMR if there are fewer than max_mRMR_features
if temp_X_train.shape[1] <= max_mRMR_features:
filter_features = temp_X_train.columns
else:
# mRMR
filter_feat_selector = mRMR_feature_select(temp_X_train.to_numpy(),y_train.to_numpy(),
num_features_to_select=max_mRMR_features,
K_MAX=1000,n_jobs=-1,verbose=False)
# Get selected features
filter_features = temp_X_train.columns[filter_feat_selector]
# Save features selected
mRMR_features.append(list(filter_features))
# Check that no features were excluded when checking for wrongly chosen
total_features = 0
for feat_c in range(len(feat_counter)):
total_features += feat_counter[feat_c][1]
assert total_features == X_train.shape[1]
# Part 2 with second mRMR and classifiers in loop
k_fold = 10
skf2 = StratifiedKFold(n_splits=k_fold, shuffle=True, random_state=(rep+1)*123) # using 10% equals around 21 test subjects
# I am also using converting it to an iterable list instad of the generator to promote reuse
Inner_CV = []
for train_index2, test_index2 in skf2.split(train_df,study_group_status.iloc[train_index]):
Inner_CV.append([train_index2,test_index2])
# SVM with L2-norm (it is by default squared)
# Prepare hyper parameters
exponent = np.linspace(-3,1,9)
exponent = np.round(exponent,5) # sometimes linspace are not so exact
C_parameter_SVM = np.power(np.array([10]*len(exponent)),exponent)
kernels = ["linear"]
# rbf overfit easier, whereas linear empirically works better in high D data
# Sequential Feature Selection for Logistic Regression
# In-built model selection CV
# The L2 is the inverse of C for LogReg
exponent = np.linspace(-3,1,9)
exponent = np.round(exponent,5) # sometimes linspace are not so exact
C_parameter_LogReg = np.power(np.array([10]*len(exponent)),exponent)
# Random forest classifier
# Prepare hyper parameters
trees = np.array([10, 100, 500, 1000])
depth = np.linspace(1,2,2) # using more depth leads to a lot of overfitting
# Prepare arrays for validation errors
val_err_svm = np.zeros((len(n_feat_mRMR1),len(n_feat_mRMR2),len(C_parameter_SVM),len(kernels)))
val_feat_svm = np.zeros((len(n_feat_mRMR1),len(n_feat_mRMR2),len(C_parameter_SVM),len(kernels),max_mRMR_features2))
val_feat_svm.fill(np.nan)
val_err_logreg = np.zeros((len(n_feat_mRMR1),len(n_feat_mRMR2),len(C_parameter_LogReg)))
val_feat_logreg = np.zeros((len(n_feat_mRMR1),len(n_feat_mRMR2),len(C_parameter_LogReg),max_mRMR_features2))
val_feat_logreg.fill(np.nan)
val_err_rf = np.zeros((len(n_feat_mRMR1),len(n_feat_mRMR2),len(trees),len(depth)))
val_feat_rf = np.zeros((len(n_feat_mRMR1),len(n_feat_mRMR2),len(trees),len(depth),max_mRMR_features2))
val_feat_rf.fill(np.nan)
val_feat_import_rf = np.zeros((len(n_feat_mRMR1),len(n_feat_mRMR2),len(trees),len(depth),max_mRMR_features2))
val_feat_import_rf.fill(np.nan)
for n1 in range(len(n_feat_mRMR1)):
# Get the features from first mRMR
mRMR_features1 = [item for sublist in mRMR_features for item in sublist[0:n_feat_mRMR1[n1]]]
# print("{} features are left after first mRMR filtering".format(len(mRMR_features1)))
X_train_mRMR = X_train[mRMR_features1].copy()
# Second round of mRMR on the reduced list
filter_feat_selector = mRMR_feature_select(X_train_mRMR.to_numpy(),y_train.to_numpy(),
num_features_to_select=max_mRMR_features2,
K_MAX=1000,n_jobs=-1,verbose=False)
for n2 in range(len(n_feat_mRMR2)):
# Get selected features from mRMR2
filter_features = X_train_mRMR.columns[filter_feat_selector[0:n_feat_mRMR2[n2]]]
X_train_mRMR2 = X_train_mRMR[filter_features].copy()
# SVM with recursive feature elemination
# Stratified CV to find best regularization strength and number of features
# CV is built-in RFECV
min_features_to_select = 1
# Using normal for loop as RFECV has inbuilt multiprocessing
for C in range(len(C_parameter_SVM)):
for K in range(len(kernels)):
# Define the model
svc = SVC(C=C_parameter_SVM[C], kernel=kernels[K], tol=1e-3, cache_size=4000)
# Perform recurive feature elimination with in-built CV
rfecv = RFECV(estimator=svc, n_jobs=-1, scoring="balanced_accuracy",
cv=Inner_CV,
min_features_to_select=min_features_to_select)
rfecv.fit(X_train_mRMR2,y_train.to_numpy().ravel())
# Get CV score
err = rfecv.grid_scores_[rfecv.n_features_-min_features_to_select]
# Save results for hyperparameter optimization
val_err_svm[n1,n2,C,K] = err
# Save the features chosen by RFE based on index from pre mRMR
rfe_feat_idx = np.where(X_train.columns.isin(filter_features[rfecv.support_]))[0]
val_feat_svm[n1,n2,C,K,0:len(rfe_feat_idx)] = rfe_feat_idx
# print("Finished SVM run {} out of {}".format(C+1,np.prod([len(C_parameter_SVM),len(kernels)])))
# Logistic regression with sequential forward selection
# In-bult CV
k_features = np.arange(1,n_feat_mRMR2[n2]+1,1) # try up to all features
inner_cv_scores = np.zeros((len(C_parameter_LogReg),len(k_features)))
for C in range(len(C_parameter_LogReg)):
LogReg = LogisticRegression(penalty="l2", C=C_parameter_LogReg[C],
max_iter = 50000)
sfs = SFS(LogReg, k_features = (k_features[0],k_features[-1]),
forward = True, scoring = "balanced_accuracy",
verbose = 0, floating = False,
cv = Inner_CV, n_jobs = -1)
sfs = sfs.fit(X_train_mRMR2, y_train.to_numpy().ravel())
# Save CV scores for each SFS step
for feat in range(len(k_features)):
inner_cv_scores[C,feat] = sfs.get_metric_dict()[k_features[feat]]["avg_score"]
# Find the best number of features
# I am rounding to make it more smooth and disregard small improvements
K = np.where(np.round(inner_cv_scores[C,:],2)==np.max(np.round(inner_cv_scores[C,:],2)))[0]
if len(K) > 1:
K = K[np.where(K == np.min(K))[0]]
err = inner_cv_scores[C,K]
# Save validation error and features
val_err_logreg[n1,n2,C] = err
# Save the features chosen by RFE based on index from pre mRMR
sfs_feat_idx = np.where(X_train.columns.isin(sfs.subsets_[int(K)+1]["feature_names"]))[0] # K+1 since it starts with 1 feature
val_feat_logreg[n1,n2,C,0:len(sfs_feat_idx)] = sfs_feat_idx
# print("Finished LogReg run {} out of {}".format(C+1,len(C_parameter_LogReg)))
# Random forest
model_test_err = np.zeros((k_fold,len(trees),len(depth)))
RF_feat_import = np.zeros((k_fold,len(trees),len(depth),X_train_mRMR2.shape[1]))
counter = 0
for cv in range(k_fold):
# Retrieve CV indices
train_index_rf = Inner_CV[cv][0]
test_index_rf = Inner_CV[cv][1]
# Retrieve datasets
X_train2 = X_train_mRMR2.iloc[train_index_rf]; X_test2 = X_train_mRMR2.iloc[test_index_rf]
y_train2 = y_train.to_numpy().ravel()[train_index_rf]; y_test2 = y_train.to_numpy().ravel()[test_index_rf]
for t in range(len(trees)):
for d in range(len(depth)):
RF = RandomForestClassifier(n_estimators=trees[t], max_features="sqrt", max_depth=depth[d],
n_jobs=-1, random_state=None)
RF.fit(X_train2.to_numpy(),y_train2)
# Validation error
RF_y = RF.predict(X_test2.to_numpy())
err = balanced_accuracy_score(y_test2, RF_y)
# Save to array
model_test_err[cv,t,d] = err
# Save feature importance array
RF_feat_import[cv,t,d] = RF.feature_importances_
counter += 1
# print("Finished RF run {} out of {}".format(counter, k_fold))
# Average the errors over the CV folds
model_test_err_mean = np.mean(model_test_err, axis=0)
val_err_rf[n1,n2,:,:] = model_test_err_mean
# Average the feature importances over the CV folds
RF_feat_import = np.mean(RF_feat_import, axis=0)
val_feat_import_rf[n1,n2,:,:,0:n_feat_mRMR2[n2]] = RF_feat_import
# Save the features used by the RF based on index from pre mRMR
rf_feat_idx = np.where(X_train.columns.isin(filter_features))[0]
val_feat_rf[n1,n2,:,:,0:len(rf_feat_idx)] = rf_feat_idx
print("Finished {}% of total run".format((n1+1)*1/len(n_feat_mRMR1)*100))
# Choose the optimal parameters
### SVM
n1, n2, C, K = np.where(val_err_svm==np.max(val_err_svm))
if len(C) > 1:
print("There are multiple SVM runs with the same validation error")
print("Choosing run with fewest features to alleviate overfitting")
rfe_feat_len = []
for i2 in range(len(C)):
n1_temp = int(n1[i2])
n2_temp = int(n2[i2])
C_temp = int(C[i2])
K_temp = int(K[i2])
temp_feats = val_feat_svm[n1_temp,n2_temp,C_temp,K_temp][~np.isnan(val_feat_svm[n1_temp,n2_temp,C_temp,K_temp])].astype(int)
rfe_feat_len.append(len(temp_feats))
rfe_feat_min = np.where(rfe_feat_len==np.min(rfe_feat_len))[0]
if len(rfe_feat_min) > 1:
print("Multiple SVM runs with same number of fewest features")
print("Choosing run with lowest C (highest regularization) to alleviate overfitting")
C_min = np.argmin(C[rfe_feat_min])
n1, n2, C, K = [int(n1[rfe_feat_min][C_min]), int(n2[rfe_feat_min][C_min]),
int(C[rfe_feat_min][C_min]), int(K[rfe_feat_min][C_min])]
else:
n1, n2, C, K = [int(n1[int(rfe_feat_min)]), int(n2[int(rfe_feat_min)]),
int(C[int(rfe_feat_min)]), int(K[int(rfe_feat_min)])]
else:
n1, n2, C, K = [int(n1), int(n2), int(C), int(K)]
mRMR_chosen1 = n_feat_mRMR1[n1]
mRMR_chosen2 = n_feat_mRMR2[n2]
C_chosen = C_parameter_SVM[C]
K_chosen = kernels[K]
# Save the best validation erro
val_error = val_err_svm[n1, n2, C, K]
# Get the subsets chosen
chosen_feats = val_feat_svm[n1,n2,C,K][~np.isnan(val_feat_svm[n1,n2,C,K])].astype(int)
rfe_features = list(X_train.columns[chosen_feats])
n_final_feat = len(rfe_features)
x_train_mRMR2_rfe = X_train[rfe_features]
# Fit on all training data
model = SVC(C=C_chosen, kernel=K_chosen, tol=1e-3, cache_size=4000)
model = model.fit(x_train_mRMR2_rfe,y_train.to_numpy().ravel())
# Get training output
model_train_y = model.predict(x_train_mRMR2_rfe)
# Get training error
train_error = balanced_accuracy_score(y_train, model_train_y)
# Get prediction of test data
y_pred = model.predict(X_test[rfe_features])
# Use model to predict class on test data
test_error = balanced_accuracy_score(y_test, y_pred)
# Save or prepare to save
accuracy_arr[rep,Outer_counter,0,:] = [train_error,val_error,test_error]
SVM_model_par = [mRMR_chosen1,mRMR_chosen2,C_chosen,n_final_feat,K_chosen]
SVM_model = model
SVM_y_pred = [[model_train_y],[y_pred]]
### LogReg with SFS
n1, n2, C = np.where(val_err_logreg==np.max(val_err_logreg))
if len(C) > 1:
print("There are multiple LogReg runs with the same validation error")
print("Choosing run with fewest features to alleviate overfitting")
sfs_feat_len = []
for i2 in range(len(C)):
n1_temp = int(n1[i2])
n2_temp = int(n2[i2])
C_temp = int(C[i2])
temp_feats = val_feat_logreg[n1_temp,n2_temp,C_temp][~np.isnan(val_feat_logreg[n1_temp,n2_temp,C_temp])].astype(int)
sfs_feat_len.append(len(temp_feats))
sfs_feat_min = np.where(sfs_feat_len==np.min(sfs_feat_len))[0]
if len(sfs_feat_min) > 1:
print("Multiple LogReg runs with same number of fewest features")
print("Choosing run with lowest C (highest regularization) to alleviate overfitting")
C_min = np.argmin(C[sfs_feat_min])
n1, n2, C = [int(n1[sfs_feat_min][C_min]), int(n2[sfs_feat_min][C_min]),
int(C[sfs_feat_min][C_min])]
else:
n1, n2, C = [int(n1[int(sfs_feat_min)]), int(n2[int(sfs_feat_min)]),
int(C[int(sfs_feat_min)])]
else:
n1, n2, C = [int(n1), int(n2), int(C)]
mRMR_chosen1 = n_feat_mRMR1[n1]
mRMR_chosen2 = n_feat_mRMR2[n2]
C_chosen = C_parameter_LogReg[C]
# Save the best validation erro
val_error = val_err_logreg[n1, n2, C]
# Get the subsets chosen
chosen_feats = val_feat_logreg[n1,n2,C][~np.isnan(val_feat_logreg[n1,n2,C])].astype(int)
sfs_features = list(X_train.columns[chosen_feats])
n_final_feat = len(sfs_features)
x_train_mRMR2_sfs = X_train[sfs_features]
# Fit on all training data
model = LogisticRegression(penalty="l2", C=C_chosen, max_iter = 50000)
model = model.fit(x_train_mRMR2_sfs,y_train.to_numpy().ravel())
# Get training output
model_train_y = model.predict(x_train_mRMR2_sfs)
# Get training error
train_error = balanced_accuracy_score(y_train, model_train_y)
# Get prediction of test data
y_pred = model.predict(X_test[sfs_features])
# Use model to predict class on test data
test_error = balanced_accuracy_score(y_test, y_pred)
# Save or prepare to save
accuracy_arr[rep,Outer_counter,1,:] = [train_error,val_error,test_error]
LogReg_model_par = [mRMR_chosen1,mRMR_chosen2,C_chosen,n_final_feat]
LogReg_model = model
LogReg_y_pred = [[model_train_y],[y_pred]]
### RF
n1, n2, t, d = np.where(val_err_rf==np.max(val_err_rf))
if len(d) > 1:
print("There are multiple RF runs with the same validation error")
print("Choosing run with lowest depth to alleviate overfitting")
d_min = np.where(d==np.min(d))[0]
if len(d_min) > 1:
print("Multiple RF runs with same number number of depths")
print("Choosing run with lowest trees to alleviate overfitting")
t_min = np.argmin(t[d_min]) # argmin just takes the first
# If there are multiple with same parameters and validation error it is most likely the same
n1, n2, t, d = [int(n1[d_min][t_min]), int(n2[d_min][t_min]),
int(t[d_min][t_min]), int(d[d_min][t_min])]
else:
n1, n2, t, d = [int(n1[int(d_min)]), int(n2[int(d_min)]),
int(t[int(d_min)]), int(d[int(d_min)])]
else:
n1, n2, t, d = [int(n1), int(n2), int(t), int(d)]
mRMR_chosen1 = n_feat_mRMR1[n1]
mRMR_chosen2 = n_feat_mRMR2[n2]
t_chosen = trees[t]
d_chosen = depth[d]
# Save the best validation erro
val_error = val_err_rf[n1, n2, t, d]
# Get the chosen features and feature importances
chosen_feats = val_feat_rf[n1,n2,t,d][~np.isnan(val_feat_rf[n1,n2,t,d])].astype(int)
rf_features = list(X_train.columns[chosen_feats])
n_final_feat = len(rf_features)
x_train_mRMR2_rf = X_train[rf_features]
rf_feat_importances = val_feat_import_rf[n1,n2,t,d][~np.isnan(val_feat_import_rf[n1,n2,t,d])]
# Fit on all training data
model = RandomForestClassifier(n_estimators=t_chosen, max_features="sqrt", max_depth=d_chosen,
n_jobs=-1, random_state=None)
model = model.fit(x_train_mRMR2_rf,y_train.to_numpy().ravel())
# Get training output
model_train_y = model.predict(x_train_mRMR2_rf)
# Get training error
train_error = balanced_accuracy_score(y_train, model_train_y)
# Get prediction of test data
y_pred = model.predict(X_test[rf_features])
# Use model to predict class on test data
test_error = balanced_accuracy_score(y_test, y_pred)
# Save or prepare to save
accuracy_arr[rep,Outer_counter,2,:] = [train_error,val_error,test_error]
RF_model_par = [mRMR_chosen1,mRMR_chosen2,t_chosen,d_chosen]
RF_model = model
RF_y_pred = [[model_train_y],[y_pred]]
# Save the rest
model_par0.append([SVM_model_par,LogReg_model_par,RF_model_par])
final_models0.append([SVM_model,LogReg_model,RF_model])
final_features0.append([rfe_features, sfs_features, [rf_features,rf_feat_importances]])
final_y_preds0.append([SVM_y_pred,LogReg_y_pred,RF_y_pred])
# Move counter
Outer_counter += 1
print("Finished outer fold {} out of {} for rep: {}".format(Outer_counter,k_out,rep+1))
# Save results from all outer folds
model_par.append(model_par0)
final_models.append(final_models0)
final_features.append(final_features0)
final_y_preds.append(final_y_preds0)
# Save results to file
Rep_mRMR2_SVM_LogReg_RF = [accuracy_arr, model_par, final_models, final_features, final_y_preds]
# Run with variable feat in first and second mRMR
with open(Model_savepath+"Rep10_10x10CV_mRMR2_SVM_LogReg_RF_Group_Status_results_010222.pkl", "wb") as filehandle:
pickle.dump(Rep_mRMR2_SVM_LogReg_RF, filehandle)
# Get current time
c_time2 = time.localtime()
c_time2 = time.strftime("%a %d %b %Y %H:%M:%S", c_time2)
print("Started", c_time1, "\nCurrent Time",c_time2)
# Print total progress
print("Finished outer fold repetition {} out of {}".format(rep+1,n_repetitions))
# %% ML predction of PTSD status using each feature separately
# mRMR -> mRMR -> SVM/RF/LogReg with L2
# 10 repetitions of 10 fold outer and 10 fold inner two-layer CV
# With concurrent processes for faster computation
EEG_features_df, EEG_features_name_list = load_features_df()
# Add group status
Group_status = np.array([0]*EEG_features_df.shape[0]) # CTRL = 0
Group_status[np.array([i in cases for i in EEG_features_df["Subject_ID"]])] = 1 # PTSD = 1
EEG_features_df.insert(1, "Group_status", Group_status)
# Subject info columns
Subject_info_cols = list(EEG_features_df.columns[0:2])
n_subject_info_cols = len(Subject_info_cols)
n_discrete_cols = 2
# To ensure proper stratification into train/test set I will stratify using group status and study status
# A variable that encodes for this is created
n_studies = 3
study_group_status = EEG_features_df["Group_status"].copy()
for i in range(n_studies):
# Get study index
study_idx = (EEG_features_df["Subject_ID"]>=(i+1)*100000)&(EEG_features_df["Subject_ID"]<(i+2)*100000)
# Assign label
study_group_status[(study_idx)&(EEG_features_df["Group_status"]==0)] = 2*i # CTRL
study_group_status[(study_idx)&(EEG_features_df["Group_status"]==1)] = 2*i+1 # PTSD
# Target variable
Target = ["Group_status"]
Target_col = EEG_features_df.iloc[:,1:].columns.isin(Target)
Target_col_idx = np.where(Target_col == True)[0]
# Make 3 models and save them to use enseemble in the end
CLF_models = ["SVM", "LogReg", "RF"]
n_models = len(CLF_models)
# Repeat the classification to see if I just got a lucky seed
n_repetitions = 10
k_out = 10
# accuracy_arr = np.zeros((n_repetitions,k_out,len(EEG_features_name_list),n_models,3))
Ind_feat_clf = [0]*n_repetitions
np.random.seed(42)
# Prepare the splits beforehand to make sure the repetitions are not the same
Rep_Outer_CV = []
Rep_Outer_CV_test = [] # a list with only the test indices
for rep in range(n_repetitions):
skf = StratifiedKFold(n_splits=k_out, shuffle=True, random_state=(rep+1)*123) # using 10% equals around 21 test subjects
# I am also using converting it to an iterable list instad of the generator to promote reuse
Outer_CV = []
Outer_CV_test = []
for train_index, test_index in skf.split(EEG_features_df,study_group_status):
Outer_CV.append([train_index,test_index])
Outer_CV_test.append(test_index)
Rep_Outer_CV.append(Outer_CV)
Rep_Outer_CV_test.append(Outer_CV_test)
# Check none of the repetitions are the same by using only the test sets
# The list is first flattened
Rep_Outer_CV_test_flat = [item for sublist in Rep_Outer_CV_test for item in sublist]
# All elements are converted to strings
# This makes it easier to look for uniques, and the indices are already in numerical order
Rep_Outer_CV_test_flat_str = ["".join([str(x) for x in ele])for ele in Rep_Outer_CV_test_flat]
def allUnique(x):
seen = set()
return not any(i in seen or seen.add(i) for i in x)
assert allUnique(Rep_Outer_CV_test_flat_str)
# Get current time
c_time1 = time.localtime()
c_time1 = time.strftime("%a %d %b %Y %H:%M:%S", c_time1)
print(c_time1)
def Each_feat_10rep_10x10CV(rep):
# The outer fold CV has already been saved as lists
Outer_CV = Rep_Outer_CV[rep]
# Pre-allocate memory
Ind_feat_clf0 = []
accuracy_arr0 = np.zeros((k_out,len(EEG_features_name_list),n_models,3))
Outer_counter = 0
for train_index, test_index in Outer_CV:
test_df = EEG_features_df.iloc[test_index]
train_df = EEG_features_df.iloc[train_index]
# Training data will be standardized
standardizer = preprocessing.StandardScaler().fit(train_df.iloc[:,n_discrete_cols:])
train_df_standard = train_df.copy()
train_df_standard.iloc[:,n_discrete_cols:] = standardizer.transform(train_df_standard.iloc[:,n_discrete_cols:])
# Test data will also be standardized but using mean and std from training data
test_df_standard = test_df.copy()
test_df_standard.iloc[:,n_discrete_cols:] = standardizer.transform(test_df_standard.iloc[:,n_discrete_cols:])
# Get the training data
X_train = train_df_standard.copy().drop(Subject_info_cols, axis=1)
y_train = train_df_standard[Target]
# Get test data
X_test = test_df_standard.copy().drop(Subject_info_cols, axis=1)
y_test = test_df_standard[Target]
# Part 2 with second mRMR and classifiers in loop
k_fold = 10
skf2 = StratifiedKFold(n_splits=k_fold, shuffle=True, random_state=(rep+1)*123) # using 10% equals around 21 test subjects
# I am also using converting it to an iterable list instad of the generator to promote reuse
Inner_CV = []
for train_index2, test_index2 in skf2.split(train_df,study_group_status.iloc[train_index]):
Inner_CV.append([train_index2,test_index2])
# Prepare initial filtering of feature types to alleviate imbalance in number of features
# Use a variable that is optimized using inner CV
n_feat_mRMR1 = [20,30,40,50]
# Max features from mRMR for each eye status
max_mRMR_features = n_feat_mRMR1[-1]
eye_status = ["Eyes Closed", "Eyes Open"]
n_eye_status = len(eye_status)
feat_counter = []
# Perform mRMR for each feature type followed by RFE SVM/SFS LogReg/RF
Ind_feat_clf00 = []
for fset in range(len(EEG_features_name_list)):
temp_feat = EEG_features_name_list[fset]
other_feats = EEG_features_name_list[:fset]+EEG_features_name_list[fset+1:]
other_feats = [item for sublist in other_feats for item in sublist] # make the list flatten
# Retrieve the dataset for each feature
col_idx = np.zeros(len(X_train.columns), dtype=bool)
for fsub in range(len(temp_feat)):
temp_feat0 = temp_feat[fsub]
col_idx0 = X_train.columns.str.contains(temp_feat0+"_")
col_idx = np.logical_or(col_idx,col_idx0) # append all trues
temp_X_train = X_train.loc[:,col_idx]
# Check if any of the other features were wrongly chosen
for fcheck in range(len(other_feats)):
if any(temp_X_train.columns.str.contains(other_feats[fcheck]+"_")==True):
temp_X_train = temp_X_train.loc[:,np.invert(temp_X_train.columns.str.contains(other_feats[fcheck]))]
if temp_X_train.size == 0: # if no columns are left, e.g. imcoh when removing coh, then add features again
temp_X_train = X_train.loc[:,X_train.columns.str.contains(temp_feat0+"_")]
# Save number of original features fed into mRMR
feat_counter.append([temp_feat,temp_X_train.shape[1]])
# If there are no features selected, then skip the rest of the loop
if temp_X_train.shape[1] == 0:
continue
# Do not use mRMR if there are fewer than n_features
if temp_X_train.shape[1] <= max_mRMR_features:
filter_features = temp_X_train.columns
else:
# mRMR
filter_feat_selector = mRMR_feature_select(temp_X_train.to_numpy(),y_train.to_numpy(),
num_features_to_select=max_mRMR_features,
K_MAX=1000,n_jobs=1,verbose=False)
# Get selected features
filter_features = temp_X_train.columns[filter_feat_selector]
# SVM with L2-norm (it is by default squared)
# Prepare hyper parameters
exponent = np.linspace(-3,1,9)
exponent = np.round(exponent,5) # sometimes linspace are not so exact
C_parameter_SVM = np.power(np.array([10]*len(exponent)),exponent)
# rbf kernel overfit easier, whereas linear empirically works better in high D data
# Sequential Feature Selection for Logistic Regression
# In-built model selection CV
# The L2 is the inverse of C for LogReg
exponent = np.linspace(-3,1,9)
exponent = np.round(exponent,5) # sometimes linspace are not so exact
C_parameter_LogReg = np.power(np.array([10]*len(exponent)),exponent)
# Random forest classifier
# Prepare hyper parameters
trees = np.array([10, 100, 500, 1000])
depth = np.linspace(1,2,2) # using more depth leads to a lot of overfitting
# Prepare arrays for validation errors
val_err_svm = np.zeros((len(n_feat_mRMR1),len(C_parameter_SVM)))
val_feat_svm = np.zeros((len(n_feat_mRMR1),len(C_parameter_SVM),max_mRMR_features))
val_feat_svm.fill(np.nan)
val_err_logreg = np.zeros((len(n_feat_mRMR1),len(C_parameter_LogReg)))
val_feat_logreg = np.zeros((len(n_feat_mRMR1),len(C_parameter_LogReg),max_mRMR_features))
val_feat_logreg.fill(np.nan)
val_err_rf = np.zeros((len(n_feat_mRMR1),len(trees),len(depth)))
val_feat_rf = np.zeros((len(n_feat_mRMR1),len(trees),len(depth),max_mRMR_features))
val_feat_rf.fill(np.nan)
val_feat_import_rf = np.zeros((len(n_feat_mRMR1),len(trees),len(depth),max_mRMR_features))
val_feat_import_rf.fill(np.nan)
min_features_to_select = 1
for n1 in range(len(n_feat_mRMR1)):
mRMR_features1 = filter_features[0:n_feat_mRMR1[n1]]
# Use the selected filter features
temp_X_train_mRMR = temp_X_train[mRMR_features1]
# SVM with recursive feature elemination
# Stratified CV to find best regularization strength and number of features
# CV is built-in RFECV
for C in range(len(C_parameter_SVM)):
# Define the model
svc = SVC(C=C_parameter_SVM[C], kernel="linear", tol=1e-3, cache_size=4000)
# Perform recurive feature elimination with in-built CV
rfecv = RFECV(estimator=svc, n_jobs=1, scoring="balanced_accuracy",
cv=Inner_CV,
min_features_to_select=min_features_to_select)
rfecv.fit(temp_X_train_mRMR,y_train.to_numpy().ravel())
# Get CV score
err = rfecv.grid_scores_[rfecv.n_features_-min_features_to_select]
# Save results for hyperparameter optimization
val_err_svm[n1,C] = err
# Save the features chosen by RFE based on index from pre mRMR
rfe_feat_idx = np.where(X_train.columns.isin(mRMR_features1[rfecv.support_]))[0]
val_feat_svm[n1,C,0:len(rfe_feat_idx)] = rfe_feat_idx
# print("Finished SVM run {} out of {}".format(C+1,len(C_parameter_SVM)))
# Logistic regression with sequential forward selection
# In-bult CV
k_max = np.min([temp_X_train_mRMR.shape[1],n_feat_mRMR1[n1]])
k_features = np.arange(1,k_max+1,1) # try up to all features
inner_cv_scores = np.zeros((len(C_parameter_LogReg),len(k_features)))
for C in range(len(C_parameter_LogReg)):
LogReg = LogisticRegression(penalty="l2", C=C_parameter_LogReg[C],
max_iter = 50000)
sfs = SFS(LogReg, k_features = (k_features[0],k_features[-1]),
forward = True, scoring = "balanced_accuracy",
verbose = 0, floating = False,
cv = Inner_CV, n_jobs = 1)
sfs = sfs.fit(temp_X_train_mRMR, y_train.to_numpy().ravel())
# Save CV scores for each SFS step
for feat in range(len(k_features)):
inner_cv_scores[C,feat] = sfs.get_metric_dict()[k_features[feat]]["avg_score"]
# Find the best number of features
# I am rounding to make it more smooth and disregard small improvements
K = np.where(np.round(inner_cv_scores[C,:],2)==np.max(np.round(inner_cv_scores[C,:],2)))[0]
if len(K) > 1:
K = K[np.where(K == np.min(K))[0]]
err = inner_cv_scores[C,K]
# Save validation error and features
val_err_logreg[n1,C] = err
# Save the features chosen by RFE based on index from pre mRMR
sfs_feat_idx = np.where(X_train.columns.isin(sfs.subsets_[int(K)+1]["feature_names"]))[0] # K+1 since it starts with 1 feature
val_feat_logreg[n1,C,0:len(sfs_feat_idx)] = sfs_feat_idx
# print("Finished LogReg run {} out of {}".format(C+1,len(C_parameter_LogReg)))
# Random forest
model_test_err = np.zeros((k_fold,len(trees),len(depth)))
RF_feat_import = np.zeros((k_fold,len(trees),len(depth),temp_X_train_mRMR.shape[1]))
counter = 0
for cv in range(k_fold):
# Retrieve CV indices
train_index_rf = Inner_CV[cv][0]
test_index_rf = Inner_CV[cv][1]
# Retrieve datasets
X_train2 = temp_X_train_mRMR.iloc[train_index_rf]; X_test2 = temp_X_train_mRMR.iloc[test_index_rf]
y_train2 = y_train.to_numpy().ravel()[train_index_rf]; y_test2 = y_train.to_numpy().ravel()[test_index_rf]
for t in range(len(trees)):
for d in range(len(depth)):
RF = RandomForestClassifier(n_estimators=trees[t], max_features="sqrt", max_depth=depth[d],
n_jobs=1, random_state=None)
RF.fit(X_train2.to_numpy(),y_train2)
# Validation error
RF_y = RF.predict(X_test2.to_numpy())
err = balanced_accuracy_score(y_test2,RF_y)
# Save to array
model_test_err[cv,t,d] = err
# Save feature importance array
RF_feat_import[cv,t,d] = RF.feature_importances_
counter += 1
# print("Finished RF run {} out of {}".format(counter, k_fold))
# Average the errors over the CV folds
model_test_err_mean = np.mean(model_test_err, axis=0)
val_err_rf[n1,:,:] = model_test_err_mean
# Average the feature importances over the CV folds
RF_feat_import = np.mean(RF_feat_import, axis=0)
val_feat_import_rf[n1,:,:,0:RF_feat_import.shape[2]] = RF_feat_import
# Save the features used by the RF based on index from pre mRMR
rf_feat_idx = np.where(X_train.columns.isin(mRMR_features1))[0]
val_feat_rf[n1,:,:,0:len(rf_feat_idx)] = rf_feat_idx
# Print progress
current_progress = (n1+1)/(len(n_feat_mRMR1))*100
print("Finished {}% of inner fold optimization for feat: {}".format(current_progress,temp_feat[0]))
# Choose the optimal parameters
### SVM
n1, C = np.where(val_err_svm==np.max(val_err_svm))
if len(C) > 1:
print("There are multiple SVM runs with the same validation error")
print("Choosing run with fewest features to alleviate overfitting")
rfe_feat_len = []
for i2 in range(len(C)):
n1_temp = int(n1[i2])
C_temp = int(C[i2])
temp_feats = val_feat_svm[n1_temp,C_temp][~np.isnan(val_feat_svm[n1_temp,C_temp])].astype(int)
rfe_feat_len.append(len(temp_feats))
rfe_feat_min = np.where(rfe_feat_len==np.min(rfe_feat_len))[0]
if len(rfe_feat_min) > 1:
print("Multiple SVM runs with same number of fewest features")
print("Choosing run with lowest C (highest regularization) to alleviate overfitting")
C_min = np.argmin(C[rfe_feat_min])
n1, C = [int(n1[rfe_feat_min][C_min]),
int(C[rfe_feat_min][C_min])]
else:
n1, C = [int(n1[int(rfe_feat_min)]),
int(C[int(rfe_feat_min)])]
else:
n1, C = [int(n1), int(C)]
mRMR_chosen1 = n_feat_mRMR1[n1]
C_chosen = C_parameter_SVM[C]
# Save the best validation error
val_error = val_err_svm[n1, C]
# Get the subsets chosen
chosen_feats = val_feat_svm[n1,C][~np.isnan(val_feat_svm[n1,C])].astype(int)
rfe_features = list(X_train.columns[chosen_feats])
n_final_feat = len(rfe_features)
x_train_mRMR2_rfe = X_train[rfe_features]
# Fit on all training data
model = SVC(C=C_chosen, kernel="linear", tol=1e-3, cache_size=4000)
model = model.fit(x_train_mRMR2_rfe,y_train.to_numpy().ravel())
# Get training output
model_train_y = model.predict(x_train_mRMR2_rfe)
# Get training error
train_error = balanced_accuracy_score(y_train, model_train_y)
# Get prediction of test data
y_pred = model.predict(X_test[rfe_features])
# Use model to predict class on test data
test_error = balanced_accuracy_score(y_test, y_pred)
# Save or prepare to save
accuracy_arr0[Outer_counter,fset,0,:] = [train_error,val_error,test_error]
SVM_model_par = [mRMR_chosen1,C_chosen,n_final_feat]
SVM_model = model
SVM_y_pred = [[model_train_y],[y_pred]]
### LogReg with SFS
n1, C = np.where(val_err_logreg==np.max(val_err_logreg))
if len(C) > 1:
print("There are multiple LogReg runs with the same validation error")
print("Choosing run with fewest features to alleviate overfitting")
sfs_feat_len = []
for i2 in range(len(C)):
n1_temp = int(n1[i2])
C_temp = int(C[i2])
temp_feats = val_feat_logreg[n1_temp,C_temp][~np.isnan(val_feat_logreg[n1_temp,C_temp])].astype(int)
sfs_feat_len.append(len(temp_feats))
sfs_feat_min = np.where(sfs_feat_len==np.min(sfs_feat_len))[0]
if len(sfs_feat_min) > 1:
print("Multiple LogReg runs with same number of fewest features")
print("Choosing run with lowest C (highest regularization) to alleviate overfitting")
C_min = np.argmin(C[sfs_feat_min])
n1, C = [int(n1[sfs_feat_min][C_min]),
int(C[sfs_feat_min][C_min])]
else:
n1, C = [int(n1[int(sfs_feat_min)]),
int(C[int(sfs_feat_min)])]
else:
n1, C = [int(n1), int(C)]
mRMR_chosen1 = n_feat_mRMR1[n1]
C_chosen = C_parameter_LogReg[C]
# Save the best validation erro
val_error = val_err_logreg[n1, C]
# Get the subsets chosen
chosen_feats = val_feat_logreg[n1,C][~np.isnan(val_feat_logreg[n1,C])].astype(int)
sfs_features = list(X_train.columns[chosen_feats])
n_final_feat = len(sfs_features)
x_train_mRMR2_sfs = X_train[sfs_features]
# Fit on all training data
model = LogisticRegression(penalty="l2", C=C_chosen, max_iter = 50000)
model = model.fit(x_train_mRMR2_sfs,y_train.to_numpy().ravel())
# Get training output
model_train_y = model.predict(x_train_mRMR2_sfs)
# Get training error
train_error = balanced_accuracy_score(y_train, model_train_y)
# Get prediction of test data
y_pred = model.predict(X_test[sfs_features])
# Use model to predict class on test data
test_error = balanced_accuracy_score(y_test, y_pred)
# Save or prepare to save
accuracy_arr0[Outer_counter,fset,1,:] = [train_error,val_error,test_error]
LogReg_model_par = [mRMR_chosen1,C_chosen,n_final_feat]
LogReg_model = model
LogReg_y_pred = [[model_train_y],[y_pred]]
### RF
n1, t, d = np.where(val_err_rf==np.max(val_err_rf))
if len(d) > 1:
print("There are multiple RF runs with the same validation error")
print("Choosing run with lowest depth to alleviate overfitting")
d_min = np.where(d==np.min(d))[0]
if len(d_min) > 1:
print("Multiple RF runs with same number number of depths")
print("Choosing run with lowest trees to alleviate overfitting")
t_min = np.argmin(t[d_min]) # argmin just takes the first
# If there are multiple with same parameters and validation error it is most likely the same
n1, t, d = [int(n1[d_min][t_min]),
int(t[d_min][t_min]), int(d[d_min][t_min])]
else:
n1, t, d = [int(n1[int(d_min)]),
int(t[int(d_min)]), int(d[int(d_min)])]
else:
n1, t, d = [int(n1), int(t), int(d)]
mRMR_chosen1 = n_feat_mRMR1[n1]
t_chosen = trees[t]
d_chosen = depth[d]
# Save the best validation error
val_error = val_err_rf[n1, t, d]
# Get the chosen features and feature importances
chosen_feats = val_feat_rf[n1,t,d][~np.isnan(val_feat_rf[n1,t,d])].astype(int)
rf_features = list(X_train.columns[chosen_feats])
n_final_feat = len(rf_features)
x_train_mRMR2_rf = X_train[rf_features]
rf_feat_importances = val_feat_import_rf[n1,t,d][~np.isnan(val_feat_import_rf[n1,t,d])]
# Fit on all training data
model = RandomForestClassifier(n_estimators=t_chosen, max_features="sqrt", max_depth=d_chosen,
n_jobs=1, random_state=None)
model = model.fit(x_train_mRMR2_rf,y_train.to_numpy().ravel())
# Get training output
model_train_y = model.predict(x_train_mRMR2_rf)
# Get training error
train_error = balanced_accuracy_score(y_train, model_train_y)
# Get prediction of test data
y_pred = model.predict(X_test[rf_features])
# Use model to predict class on test data
test_error = balanced_accuracy_score(y_test, y_pred)
# Save or prepare to save
accuracy_arr0[Outer_counter,fset,2,:] = [train_error,val_error,test_error]
RF_model_par = [mRMR_chosen1,t_chosen,d_chosen]
RF_model = model
RF_y_pred = [[model_train_y],[y_pred]]
# Save the rest
model_par00 = [SVM_model_par,LogReg_model_par,RF_model_par]
final_models00 = [SVM_model,LogReg_model,RF_model]
final_features00 = [rfe_features, sfs_features, [rf_features,rf_feat_importances]]
final_y_preds00 = [SVM_y_pred,LogReg_y_pred,RF_y_pred]
data_splits00 = [X_train,y_train,X_test,y_test,standardizer.mean_,standardizer.scale_]
# Save the results for the specific feature
res = [temp_feat, model_par00, final_models00, final_features00, final_y_preds00, data_splits00]
Ind_feat_clf00.append(res)
print("Finished outer fold {} out of {} for rep: {} for feat {}".format(Outer_counter+1,k_out,rep+1,temp_feat))
# Check that no features were excluded when checking for wrongly chosen
total_features = 0
for feat_c in range(len(feat_counter)):
total_features += feat_counter[feat_c][1]
assert total_features == X_train.shape[1]
# Save results over outer folds
Ind_feat_clf0.append(Ind_feat_clf00)
print("Finished outer fold {} out of {} for rep: {} for all features".format(Outer_counter+1,k_out,rep+1))
# Move counter
Outer_counter += 1
out = [accuracy_arr0, Ind_feat_clf0] # [Rep][Outer CV][Feature][Variable]
# Print total progress
print("Finished outer fold repetition {} out of {}".format(rep+1,n_repetitions))
# Get current time
c_time2 = time.localtime()
c_time2 = time.strftime("%a %d %b %Y %H:%M:%S", c_time2)
print("Started", c_time1, "\nCurrent Time",c_time2)
return rep, out
Ind_feat_clf = [0]*n_repetitions
with concurrent.futures.ProcessPoolExecutor(max_workers=10) as executor:
for rep, result in executor.map(Each_feat_10rep_10x10CV, range(n_repetitions)): # Function and arguments
Ind_feat_clf[rep] = result # [acc or parameters], for parameters: [Rep][Outer CV][Feature][Variable]
# Save results to file and overwrite when new results arrive
with open(Model_savepath+"Each_feat_Rep10_10x10CV_mRMR_SVM_LogReg_RF_results_070222.pkl", "wb") as filehandle:
pickle.dump(Ind_feat_clf, filehandle)
# %% Sparse clustering of all EEG features to look for subtypes
EEG_features_df, EEG_features_name_list = load_features_df()
# Add group status
Group_status = np.array([0]*EEG_features_df.shape[0]) # CTRL = 0
Group_status[np.array([i in cases for i in EEG_features_df["Subject_ID"]])] = 1 # PTSD = 1
EEG_features_df.insert(1, "Group_status", Group_status)
# Only use PTSD patient group
EEG_features_df2 = EEG_features_df.loc[EEG_features_df["Group_status"]==1,:]
Subject_info_cols = ["Subject_ID","Group_status"]
# Standardize the values
X = np.array(EEG_features_df2.copy().drop(Subject_info_cols, axis=1))
standardizer = preprocessing.StandardScaler().fit(X)
X = standardizer.transform(X)
# Use gridsearch and permutations to estimate gap statistic and use it to
# determine number of clusters and sparsity s
# I will use 100 permutations and test 2 to 6 clusters as Zhang 2020
max_clusters = 6
n_sparsity_feat = 20
perm_res = []
Timestamps = []
# Get current time
c_time1 = time.localtime()
c_time1 = time.strftime("%a %d %b %Y %H:%M:%S", c_time1)
print(c_time1)
Timestamps.append(c_time1)
for k in range(1,max_clusters):
# Cannot permute with 1 cluster
n_clusters = k+1
perm = pysparcl.cluster.permute_modified(X, k=n_clusters, verbose=True,
nvals=n_sparsity_feat, nperms=100)
perm_res.append(perm)
# Get current time
c_time2 = time.localtime()
c_time2 = time.strftime("%a %d %b %Y %H:%M:%S", c_time2)
print("Started", c_time1, "\nCurrent Time",c_time2)
Timestamps.append(c_time2)
# Save the results
with open(Feature_savepath+"All_feats_no_coh_plv_kmeans_perm.pkl", "wb") as file:
pickle.dump(perm_res, file)
# # Load
# with open(Feature_savepath+"All_feats_no_coh_plv_kmeans_perm.pkl", "rb") as file:
# perm_res = pickle.load(file)
# Convert results to array
perm_res_arr = np.zeros((len(perm_res)*n_sparsity_feat,4))
for i in range(len(perm_res)):
_, gaps, sdgaps, wbounds, _ = perm_res[i].values()
for i2 in range(n_sparsity_feat):
perm_res_arr[n_sparsity_feat*i+i2,0] = i+2 # cluster size
perm_res_arr[n_sparsity_feat*i+i2,1] = gaps[i2] # gap statistic
perm_res_arr[n_sparsity_feat*i+i2,2] = sdgaps[i2] # gap statistic std
perm_res_arr[n_sparsity_feat*i+i2,3] = wbounds[i2] # sparsity feature s
# For each sparsity s, determine best k using one-standard-error criterion
# Meaning the cluster and sparsity is chosen for the smallest value of k for a fixed s
# that fulfill Gap(k) >= Gap(k+1)-std(k+1)
def one_standard_deviation_search(gaps, std):
best_gaps = np.argmax(gaps)
current_gaps = gaps[best_gaps]
current_std = std[best_gaps]
current_gaps_idx = best_gaps
while (gaps[current_gaps_idx-1] >= current_gaps - current_std):
if current_gaps_idx == 0:
break
else:
current_gaps_idx -= 1
current_gaps = gaps[current_gaps_idx]
current_std = std[current_gaps_idx]
out = current_gaps, current_std, current_gaps_idx
return out
best_ks = np.zeros((n_sparsity_feat, 2))
all_s = np.unique(perm_res_arr[:,3])
plt.figure(figsize=(12,12))
for i2 in range(n_sparsity_feat):
current_s = all_s[i2]
gaps = perm_res_arr[perm_res_arr[:,3] == current_s,1]
std = perm_res_arr[perm_res_arr[:,3] == current_s,2]
_, _, idx = one_standard_deviation_search(gaps, std)
# Save to array
best_ks[i2,0] = current_s
best_ks[i2,1] = perm_res_arr[perm_res_arr[:,3] == current_s,0][idx]
# Plot gap
plt.errorbar(perm_res_arr[perm_res_arr[:,3] == current_s,0].astype("int"),
gaps, yerr=std, capsize=5, label = np.round(current_s,3))
plt.title("Gap statistic for different fixed s")
plt.legend(loc=1)
plt.xlabel("Number of clusters")
plt.ylabel("Gap statistic")
best_k = int(scipy.stats.mode(best_ks[:,1])[0])
# Determine s using fixed k as lowest s within 1 std of max gap statistic
# According to Witten & Tibshirani, 2010
best_gaps_idx = np.argmax(perm_res_arr[perm_res_arr[:,0] == best_k,1])
best_gaps = perm_res_arr[perm_res_arr[:,0] == best_k,1][best_gaps_idx]
best_gaps_std = perm_res_arr[perm_res_arr[:,0] == best_k,2][best_gaps_idx]
one_std_crit = perm_res_arr[perm_res_arr[:,0] == best_k,1]>=best_gaps-best_gaps_std
best_s = np.array([perm_res_arr[perm_res_arr[:,0] == best_k,3][one_std_crit][0]])
# Perform clustering with k clusters
sparcl = pysparcl.cluster.kmeans(X, k=best_k, wbounds=best_s)[0]
# # Save the results
# with open(Feature_savepath+"All_feats_sparse_kmeans.pkl", "wb") as file:
# pickle.dump(sparcl, file)
with open(Feature_savepath+"All_feats_sparse_kmeans.pkl", "rb") as file:
sparcl = pickle.load(file)
# %% Sparse kmeans -> mRMR -> SVM/RF/LogReg with L2
# Using sparse kmeans selected features, only eyes closed as they were primarily chosen
# Prediction of subtype 1
EEG_features_df, EEG_features_name_list = load_features_df()
with open(Feature_savepath+"All_feats_sparse_kmeans.pkl", "rb") as file:
sparcl = pickle.load(file)
# Use concatenated 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']]
# Add group status
Group_status = np.array([0]*EEG_features_df.shape[0]) # CTRL = 0
Group_status[np.array([i in cases for i in EEG_features_df["Subject_ID"]])] = 1 # PTSD = 1
EEG_features_df.insert(1, "Group_status", Group_status)
# Only keep PTSD subtype 1 by dropping subtype 2
subtype1 = np.array(Subject_id)[PTSD_idx][sparcl["cs"]==0]
subtype2 = np.array(Subject_id)[PTSD_idx][sparcl["cs"]==1]
EEG_features_df = EEG_features_df.set_index("Subject_ID")
EEG_features_df = EEG_features_df.drop(subtype2)
EEG_features_df = EEG_features_df.reset_index()
# Check it was dropped correctly
assert all(subtype1 == EEG_features_df.loc[EEG_features_df["Group_status"] == 1,"Subject_ID"])
# Subject info columns
Subject_info_cols = list(EEG_features_df.columns[0:2])
n_subject_info_cols = len(Subject_info_cols)
n_discrete_cols = 2
# Get features from sparse kmeans
nonzero_idx = sparcl["ws"].nonzero()
sparcl_features = EEG_features_df.copy().drop(Subject_info_cols, axis=1).columns[nonzero_idx]
sum(sparcl_features.str.contains("Eyes Open"))/len(sparcl_features) # less than 3% are EO
# Only use eyes closed (2483 features)
sparcl_features = sparcl_features[sparcl_features.str.contains("Eyes Closed")]
EEG_features_df = pd.concat([EEG_features_df[Subject_info_cols],EEG_features_df[sparcl_features]],axis=1)
# To ensure proper stratification into train/test set I will stratify using group status and study status
# A variable that encodes for this is created
n_studies = 3
study_group_status = EEG_features_df["Group_status"].copy()
for i in range(n_studies):
# Get study index
study_idx = (EEG_features_df["Subject_ID"]>=(i+1)*100000)&(EEG_features_df["Subject_ID"]<(i+2)*100000)
# Assign label
study_group_status[(study_idx)&(EEG_features_df["Group_status"]==0)] = 2*i # CTRL
study_group_status[(study_idx)&(EEG_features_df["Group_status"]==1)] = 2*i+1 # PTSD
# Target variable
Target = ["Group_status"]
Target_col = EEG_features_df.iloc[:,1:].columns.isin(Target)
Target_col_idx = np.where(Target_col == True)[0]
# Make 3 models and save them to use enseemble in the end
CLF_models = ["SVM", "LogReg", "RF"]
n_models = len(CLF_models)
# Repeat the classification to see if I just got a lucky seed
n_repetitions = 10
k_out = 10
accuracy_arr = np.zeros((n_repetitions,k_out,n_models,3))
model_par = []
final_models = []
final_features = []
final_y_preds = []
np.random.seed(42)
# Prepare the splits beforehand to make sure the repetitions are not the same
Rep_Outer_CV = []
Rep_Outer_CV_test = [] # a list with only the test indices
for rep in range(n_repetitions):
skf = StratifiedKFold(n_splits=k_out, shuffle=True, random_state=(rep+1)*123) # using 10% equals around 21 test subjects
# I am also using converting it to an iterable list instad of the generator to promote reuse
Outer_CV = []
Outer_CV_test = []
for train_index, test_index in skf.split(EEG_features_df,study_group_status):
Outer_CV.append([train_index,test_index])
Outer_CV_test.append(test_index)
Rep_Outer_CV.append(Outer_CV)
Rep_Outer_CV_test.append(Outer_CV_test)
# Check none of the repetitions are the same by using only the test sets
# The list is first flattened
Rep_Outer_CV_test_flat = [item for sublist in Rep_Outer_CV_test for item in sublist]
# All elements are converted to strings
# This makes it easier to look for uniques, and the indices are already in numerical order
Rep_Outer_CV_test_flat_str = ["".join([str(x) for x in ele])for ele in Rep_Outer_CV_test_flat]
def allUnique(x):
seen = set()
return not any(i in seen or seen.add(i) for i in x)
assert allUnique(Rep_Outer_CV_test_flat_str)
# Get current time
c_time1 = time.localtime()
c_time1 = time.strftime("%a %d %b %Y %H:%M:%S", c_time1)
print(c_time1)
for rep in range(n_repetitions):
# The outer fold CV has already been saved as lists
Outer_CV = Rep_Outer_CV[rep]
# Pre-allocate memory
model_par0 = []
final_models0 = []
final_features0 = []
final_y_preds0 = []
Outer_counter = 0
for train_index, test_index in Outer_CV:
test_df = EEG_features_df.iloc[test_index]
train_df = EEG_features_df.iloc[train_index]
# Training data will be standardized
standardizer = preprocessing.StandardScaler().fit(train_df.iloc[:,n_discrete_cols:])
train_df_standard = train_df.copy()
train_df_standard.iloc[:,n_discrete_cols:] = standardizer.transform(train_df_standard.iloc[:,n_discrete_cols:])
# Test data will also be standardized but using mean and std from training data
test_df_standard = test_df.copy()
test_df_standard.iloc[:,n_discrete_cols:] = standardizer.transform(test_df_standard.iloc[:,n_discrete_cols:])
# Get the training data
X_train = train_df_standard.copy().drop(Subject_info_cols, axis=1)
y_train = train_df_standard[Target]
# Get test data
X_test = test_df_standard.copy().drop(Subject_info_cols, axis=1)
y_test = test_df_standard[Target]
# Prepare initial filtering of feature types to alleviate imbalance in number of features
# Use a variable that is optimized using inner CV
n_feat_mRMR1 = [30,40,50,60]
# Max features from mRMR for each eye status
max_mRMR_features = n_feat_mRMR1[-1]
# mRMR on the all sparse kmeans selected features
filter_feat_selector = mRMR_feature_select(X_train.to_numpy(),y_train.to_numpy(),
num_features_to_select=max_mRMR_features,
K_MAX=1000,n_jobs=-1,verbose=False)
# Part 2 with second mRMR and classifiers in loop
k_fold = 10
skf2 = StratifiedKFold(n_splits=k_fold, shuffle=True, random_state=(rep+1)*123) # using 10% equals around 21 test subjects
# I am also using converting it to an iterable list instad of the generator to promote reuse
Inner_CV = []
for train_index2, test_index2 in skf2.split(train_df,study_group_status.iloc[train_index]):
Inner_CV.append([train_index2,test_index2])
# SVM with L2-norm (it is by default squared)
# Prepare hyper parameters
exponent = np.linspace(-3,1,9)
exponent = np.round(exponent,5) # sometimes linspace are not so exact
C_parameter_SVM = np.power(np.array([10]*len(exponent)),exponent)
kernels = ["linear"] # error when using rbf and RFECV
# rbf also overfit easier, whereas linear empirically works better in high D data
# Sequential Feature Selection for Logistic Regression
# In-built model selection CV
# The L2 is the inverse of C for LogReg
exponent = np.linspace(-3,1,9)
exponent = np.round(exponent,5) # sometimes linspace are not so exact
C_parameter_LogReg = np.power(np.array([10]*len(exponent)),exponent)
# Random forest classifier
# Prepare hyper parameters
trees = np.array([10, 100, 500, 1000])
depth = np.linspace(1,2,2) # using more depth leads to a lot of overfitting
# Prepare arrays for validation errors
val_err_svm = np.zeros((len(n_feat_mRMR1),len(C_parameter_SVM),len(kernels)))
val_feat_svm = np.zeros((len(n_feat_mRMR1),len(C_parameter_SVM),len(kernels),max_mRMR_features))
val_feat_svm.fill(np.nan)
val_err_logreg = np.zeros((len(n_feat_mRMR1),len(C_parameter_LogReg)))
val_feat_logreg = np.zeros((len(n_feat_mRMR1),len(C_parameter_LogReg),max_mRMR_features))
val_feat_logreg.fill(np.nan)
val_err_rf = np.zeros((len(n_feat_mRMR1),len(trees),len(depth)))
val_feat_rf = np.zeros((len(n_feat_mRMR1),len(trees),len(depth),max_mRMR_features))
val_feat_rf.fill(np.nan)
val_feat_import_rf = np.zeros((len(n_feat_mRMR1),len(trees),len(depth),max_mRMR_features))
val_feat_import_rf.fill(np.nan)
for n1 in range(len(n_feat_mRMR1)):
# Get selected features from mRMR
filter_features = X_train.columns[filter_feat_selector[0:n_feat_mRMR1[n1]]]
X_train_mRMR = X_train[filter_features].copy()
# SVM with recursive feature elemination
# Stratified CV to find best regularization strength and number of features
# CV is built-in RFECV
min_features_to_select = 1
# Using normal for loop as RFECV has inbuilt multiprocessing
for C in range(len(C_parameter_SVM)):
for K in range(len(kernels)):
# Define the model
svc = SVC(C=C_parameter_SVM[C], kernel=kernels[K], tol=1e-3, cache_size=4000)
# Perform recurive feature elimination with in-built CV
rfecv = RFECV(estimator=svc, n_jobs=-1, scoring="balanced_accuracy",
cv=Inner_CV,
min_features_to_select=min_features_to_select)
rfecv.fit(X_train_mRMR,y_train.to_numpy().ravel())
# Get CV score
err = rfecv.grid_scores_[rfecv.n_features_-min_features_to_select]
# Save results for hyperparameter optimization
val_err_svm[n1,C,K] = err
# Save the features chosen by RFE based on index from pre mRMR
rfe_feat_idx = np.where(X_train.columns.isin(filter_features[rfecv.support_]))[0]
val_feat_svm[n1,C,K,0:len(rfe_feat_idx)] = rfe_feat_idx
# print("Finished SVM run {} out of {}".format(C+1,np.prod([len(C_parameter_SVM),len(kernels)])))
# Logistic regression with sequential forward selection
# In-bult CV
k_features = np.arange(1,n_feat_mRMR1[n1]+1,1) # try up to all features
inner_cv_scores = np.zeros((len(C_parameter_LogReg),len(k_features)))
for C in range(len(C_parameter_LogReg)):
LogReg = LogisticRegression(penalty="l2", C=C_parameter_LogReg[C],
max_iter = 50000)
sfs = SFS(LogReg, k_features = (k_features[0],k_features[-1]),
forward = True, scoring = "balanced_accuracy",
verbose = 0, floating = False,
cv = Inner_CV, n_jobs = -1)
sfs = sfs.fit(X_train_mRMR, y_train.to_numpy().ravel())
# Save CV scores for each SFS step
for feat in range(len(k_features)):
inner_cv_scores[C,feat] = sfs.get_metric_dict()[k_features[feat]]["avg_score"]
# Find the best number of features
# I am rounding to make it more smooth and disregard small improvements
K = np.where(np.round(inner_cv_scores[C,:],2)==np.max(np.round(inner_cv_scores[C,:],2)))[0]
if len(K) > 1:
K = K[np.where(K == np.min(K))[0]]
err = inner_cv_scores[C,K]
# Save validation error and features
val_err_logreg[n1,C] = err
# Save the features chosen by RFE based on index from pre mRMR
sfs_feat_idx = np.where(X_train.columns.isin(sfs.subsets_[int(K)+1]["feature_names"]))[0] # K+1 since it starts with 1 feature
val_feat_logreg[n1,C,0:len(sfs_feat_idx)] = sfs_feat_idx
# print("Finished LogReg run {} out of {}".format(C+1,len(C_parameter_LogReg)))
# Random forest
model_test_err = np.zeros((k_fold,len(trees),len(depth)))
RF_feat_import = np.zeros((k_fold,len(trees),len(depth),X_train_mRMR.shape[1]))
counter = 0
for cv in range(k_fold):
# Retrieve CV indices
train_index_rf = Inner_CV[cv][0]
test_index_rf = Inner_CV[cv][1]
# Retrieve datasets
X_train2 = X_train_mRMR.iloc[train_index_rf]; X_test2 = X_train_mRMR.iloc[test_index_rf]
y_train2 = y_train.to_numpy().ravel()[train_index_rf]; y_test2 = y_train.to_numpy().ravel()[test_index_rf]
for t in range(len(trees)):
for d in range(len(depth)):
RF = RandomForestClassifier(n_estimators=trees[t], max_features="sqrt", max_depth=depth[d],
n_jobs=-1, random_state=None)
RF.fit(X_train2.to_numpy(),y_train2)
# Validation error
RF_y = RF.predict(X_test2.to_numpy())
err = balanced_accuracy_score(y_test2, RF_y)
# Save to array
model_test_err[cv,t,d] = err
# Save feature importance array
RF_feat_import[cv,t,d] = RF.feature_importances_
counter += 1
# print("Finished RF run {} out of {}".format(counter, k_fold))
# Average the errors over the CV folds
model_test_err_mean = np.mean(model_test_err, axis=0)
val_err_rf[n1,:,:] = model_test_err_mean
# Average the feature importances over the CV folds
RF_feat_import = np.mean(RF_feat_import, axis=0)
val_feat_import_rf[n1,:,:,0:n_feat_mRMR1[n1]] = RF_feat_import
# Save the features used by the RF based on index from pre mRMR
rf_feat_idx = np.where(X_train.columns.isin(filter_features))[0]
val_feat_rf[n1,:,:,0:len(rf_feat_idx)] = rf_feat_idx
print("Finished {}% of total run".format((n1+1)*1/len(n_feat_mRMR1)*100))
# Choose the optimal parameters
### SVM
n1, C, K = np.where(val_err_svm==np.max(val_err_svm))
if len(C) > 1:
print("There are multiple SVM runs with the same validation error")
print("Choosing run with fewest features to alleviate overfitting")
rfe_feat_len = []
for i2 in range(len(C)):
n1_temp = int(n1[i2])
C_temp = int(C[i2])
K_temp = int(K[i2])
temp_feats = val_feat_svm[n1_temp,C_temp,K_temp][~np.isnan(val_feat_svm[n1_temp,C_temp,K_temp])].astype(int)
rfe_feat_len.append(len(temp_feats))
rfe_feat_min = np.where(rfe_feat_len==np.min(rfe_feat_len))[0]
if len(rfe_feat_min) > 1:
print("Multiple SVM runs with same number of fewest features")
print("Choosing run with lowest C (highest regularization) to alleviate overfitting")
C_min = np.argmin(C[rfe_feat_min])
n1, C, K = [int(n1[rfe_feat_min][C_min]),
int(C[rfe_feat_min][C_min]), int(K[rfe_feat_min][C_min])]
else:
n1, C, K = [int(n1[int(rfe_feat_min)]),
int(C[int(rfe_feat_min)]), int(K[int(rfe_feat_min)])]
else:
n1, C, K = [int(n1), int(C), int(K)]
mRMR_chosen1 = n_feat_mRMR1[n1]
C_chosen = C_parameter_SVM[C]
K_chosen = kernels[K]
# Save the best validation erro
val_error = val_err_svm[n1, C, K]
# Get the subsets chosen
chosen_feats = val_feat_svm[n1,C,K][~np.isnan(val_feat_svm[n1,C,K])].astype(int)
rfe_features = list(X_train.columns[chosen_feats])
n_final_feat = len(rfe_features)
x_train_mRMR_rfe = X_train[rfe_features]
# Fit on all training data
model = SVC(C=C_chosen, kernel=K_chosen, tol=1e-3, cache_size=4000)
model = model.fit(x_train_mRMR_rfe,y_train.to_numpy().ravel())
# Get training output
model_train_y = model.predict(x_train_mRMR_rfe)
# Get training error
train_error = balanced_accuracy_score(y_train, model_train_y)
# Get prediction of test data
y_pred = model.predict(X_test[rfe_features])
# Use model to predict class on test data
test_error = balanced_accuracy_score(y_test, y_pred)
# Save or prepare to save
accuracy_arr[rep,Outer_counter,0,:] = [train_error,val_error,test_error]
SVM_model_par = [mRMR_chosen1,C_chosen,n_final_feat,K_chosen]
SVM_model = model
SVM_y_pred = [[model_train_y],[y_pred]]
### LogReg with SFS
n1, C = np.where(val_err_logreg==np.max(val_err_logreg))
if len(C) > 1:
print("There are multiple LogReg runs with the same validation error")
print("Choosing run with fewest features to alleviate overfitting")
sfs_feat_len = []
for i2 in range(len(C)):
n1_temp = int(n1[i2])
C_temp = int(C[i2])
temp_feats = val_feat_logreg[n1_temp,C_temp][~np.isnan(val_feat_logreg[n1_temp,C_temp])].astype(int)
sfs_feat_len.append(len(temp_feats))
sfs_feat_min = np.where(sfs_feat_len==np.min(sfs_feat_len))[0]
if len(sfs_feat_min) > 1:
print("Multiple LogReg runs with same number of fewest features")
print("Choosing run with lowest C (highest regularization) to alleviate overfitting")
C_min = np.argmin(C[sfs_feat_min])
n1, C = [int(n1[sfs_feat_min][C_min]),
int(C[sfs_feat_min][C_min])]
else:
n1, C = [int(n1[int(sfs_feat_min)]),
int(C[int(sfs_feat_min)])]
else:
n1, C = [int(n1), int(C)]
mRMR_chosen1 = n_feat_mRMR1[n1]
C_chosen = C_parameter_LogReg[C]
# Save the best validation erro
val_error = val_err_logreg[n1, C]
# Get the subsets chosen
chosen_feats = val_feat_logreg[n1,C][~np.isnan(val_feat_logreg[n1,C])].astype(int)
sfs_features = list(X_train.columns[chosen_feats])
n_final_feat = len(sfs_features)
x_train_mRMR_sfs = X_train[sfs_features]
# Fit on all training data
model = LogisticRegression(penalty="l2", C=C_chosen, max_iter = 50000)
model = model.fit(x_train_mRMR_sfs,y_train.to_numpy().ravel())
# Get training output
model_train_y = model.predict(x_train_mRMR_sfs)
# Get training error
train_error = balanced_accuracy_score(y_train, model_train_y)
# Get prediction of test data
y_pred = model.predict(X_test[sfs_features])
# Use model to predict class on test data
test_error = balanced_accuracy_score(y_test, y_pred)
# Save or prepare to save
accuracy_arr[rep,Outer_counter,1,:] = [train_error,val_error,test_error]
LogReg_model_par = [mRMR_chosen1,C_chosen,n_final_feat]
LogReg_model = model
LogReg_y_pred = [[model_train_y],[y_pred]]
### RF
n1, t, d = np.where(val_err_rf==np.max(val_err_rf))
if len(d) > 1:
print("There are multiple RF runs with the same validation error")
print("Choosing run with lowest depth to alleviate overfitting")
d_min = np.where(d==np.min(d))[0]
if len(d_min) > 1:
print("Multiple RF runs with same number number of depths")
print("Choosing run with lowest trees to alleviate overfitting")
t_min = np.argmin(t[d_min]) # argmin just takes the first
# If there are multiple with same parameters and validation error it is most likely the same
n1, t, d = [int(n1[d_min][t_min]),
int(t[d_min][t_min]), int(d[d_min][t_min])]
else:
n1, t, d = [int(n1[int(d_min)]),
int(t[int(d_min)]), int(d[int(d_min)])]
else:
n1, t, d = [int(n1), int(t), int(d)]
mRMR_chosen1 = n_feat_mRMR1[n1]
t_chosen = trees[t]
d_chosen = depth[d]
# Save the best validation error
val_error = val_err_rf[n1, t, d]
# Get the chosen features and feature importances
chosen_feats = val_feat_rf[n1,t,d][~np.isnan(val_feat_rf[n1,t,d])].astype(int)
rf_features = list(X_train.columns[chosen_feats])
n_final_feat = len(rf_features)
x_train_mRMR_rf = X_train[rf_features]
rf_feat_importances = val_feat_import_rf[n1,t,d][~np.isnan(val_feat_import_rf[n1,t,d])]
# Fit on all training data
model = RandomForestClassifier(n_estimators=t_chosen, max_features="sqrt", max_depth=d_chosen,
n_jobs=-1, random_state=None)
model = model.fit(x_train_mRMR_rf,y_train.to_numpy().ravel())
# Get training output
model_train_y = model.predict(x_train_mRMR_rf)
# Get training error
train_error = balanced_accuracy_score(y_train, model_train_y)
# Get prediction of test data
y_pred = model.predict(X_test[rf_features])
# Use model to predict class on test data
test_error = balanced_accuracy_score(y_test, y_pred)
# Save or prepare to save
accuracy_arr[rep,Outer_counter,2,:] = [train_error,val_error,test_error]
RF_model_par = [mRMR_chosen1,t_chosen,d_chosen]
RF_model = model
RF_y_pred = [[model_train_y],[y_pred]]
# Save the rest
model_par0.append([SVM_model_par,LogReg_model_par,RF_model_par])
final_models0.append([SVM_model,LogReg_model,RF_model])
final_features0.append([rfe_features, sfs_features, [rf_features,rf_feat_importances]])
final_y_preds0.append([SVM_y_pred,LogReg_y_pred,RF_y_pred])
# Move counter
Outer_counter += 1
print("Finished outer fold {} out of {} for rep: {}".format(Outer_counter,k_out,rep+1))
# Save results from all outer folds
model_par.append(model_par0)
final_models.append(final_models0)
final_features.append(final_features0)
final_y_preds.append(final_y_preds0)
# Save results to file
Rep_mRMR2_SVM_LogReg_RF = [accuracy_arr, model_par, final_models, final_features, final_y_preds]
# Run with variable feat in first and second mRMR
with open(Model_savepath+"Rep10_10x10CV_SparsKmean_mRMR_SVM_LogReg_RF_Subtype1_results_210122.pkl", "wb") as filehandle:
pickle.dump(Rep_mRMR2_SVM_LogReg_RF, filehandle)
# Get current time
c_time2 = time.localtime()
c_time2 = time.strftime("%a %d %b %Y %H:%M:%S", c_time2)
print("Started", c_time1, "\nCurrent Time",c_time2)
# Print total progress
print("Finished outer fold repetition {} out of {}".format(rep+1,n_repetitions))
# %% Each feat mRMR -> mRMR -> SVM/RF/LogReg with L2
# Using each feature (of the sparse Kmeans EC chosen ones)
# For Subtype 1
# 10 repetitions of 10 fold outer and 20 fold inner two-layer CV
# Prediction of subtype 1
EEG_features_df, EEG_features_name_list = load_features_df() # removed coh and plv and renamed Hurst in v4
with open(Feature_savepath+"All_feats_sparse_kmeans.pkl", "rb") as file:
sparcl = pickle.load(file)
# Add group status
Group_status = np.array([0]*EEG_features_df.shape[0]) # CTRL = 0
Group_status[np.array([i in cases for i in EEG_features_df["Subject_ID"]])] = 1 # PTSD = 1
EEG_features_df.insert(1, "Group_status", Group_status)
# Only keep PTSD subtype 1 by dropping subtype 2
PTSD_idx = np.array([i in cases for i in Subject_id])
CTRL_idx = np.array([not i in cases for i in Subject_id])
subtype1 = np.array(Subject_id)[PTSD_idx][sparcl["cs"]==0]
subtype2 = np.array(Subject_id)[PTSD_idx][sparcl["cs"]==1]
EEG_features_df = EEG_features_df.set_index("Subject_ID")
EEG_features_df = EEG_features_df.drop(subtype2)
EEG_features_df = EEG_features_df.reset_index()
# Check it was dropped correctly
assert all(subtype1 == EEG_features_df.loc[EEG_features_df["Group_status"] == 1,"Subject_ID"])
# Subject info columns
Subject_info_cols = list(EEG_features_df.columns[0:2])
n_subject_info_cols = len(Subject_info_cols)
n_discrete_cols = 2
# Get features from sparse kmeans
nonzero_idx = sparcl["ws"].nonzero()
sparcl_features = EEG_features_df.copy().drop(Subject_info_cols, axis=1).columns[nonzero_idx]
sum(sparcl_features.str.contains("Eyes Open"))/len(sparcl_features) # less than 3% are EO
# Only use eyes closed (2483 features)
sparcl_features = sparcl_features[sparcl_features.str.contains("Eyes Closed")]
EEG_features_df = pd.concat([EEG_features_df[Subject_info_cols],EEG_features_df[sparcl_features]],axis=1)
# To ensure proper stratification into train/test set I will stratify using group status and study status
# A variable that encodes for this is created
n_studies = 3
study_group_status = EEG_features_df["Group_status"].copy()
for i in range(n_studies):
# Get study index
study_idx = (EEG_features_df["Subject_ID"]>=(i+1)*100000)&(EEG_features_df["Subject_ID"]<(i+2)*100000)
# Assign label
study_group_status[(study_idx)&(EEG_features_df["Group_status"]==0)] = 2*i # CTRL
study_group_status[(study_idx)&(EEG_features_df["Group_status"]==1)] = 2*i+1 # PTSD
# Target variable
Target = ["Group_status"]
Target_col = EEG_features_df.iloc[:,1:].columns.isin(Target)
Target_col_idx = np.where(Target_col == True)[0]
# Make 3 models and save them to use enseemble in the end
CLF_models = ["SVM", "LogReg", "RF"]
n_models = len(CLF_models)
# Repeat the classification to see if I just got a lucky seed
n_repetitions = 10
k_out = 10
accuracy_arr = np.zeros((n_repetitions,k_out,n_models,3))
model_par = []
final_models = []
final_features = []
final_y_preds = []
np.random.seed(42)
# Prepare the splits beforehand to make sure the repetitions are not the same
Rep_Outer_CV = []
Rep_Outer_CV_test = [] # a list with only the test indices
for rep in range(n_repetitions):
skf = StratifiedKFold(n_splits=k_out, shuffle=True, random_state=(rep+1)*123) # using 10% equals around 21 test subjects
# I am also using converting it to an iterable list instad of the generator to promote reuse
Outer_CV = []
Outer_CV_test = []
for train_index, test_index in skf.split(EEG_features_df,study_group_status):
Outer_CV.append([train_index,test_index])
Outer_CV_test.append(test_index)
Rep_Outer_CV.append(Outer_CV)
Rep_Outer_CV_test.append(Outer_CV_test)
# Check none of the repetitions are the same by using only the test sets
# The list is first flattened
Rep_Outer_CV_test_flat = [item for sublist in Rep_Outer_CV_test for item in sublist]
# All elements are converted to strings
# This makes it easier to look for uniques, and the indices are already in numerical order
Rep_Outer_CV_test_flat_str = ["".join([str(x) for x in ele])for ele in Rep_Outer_CV_test_flat]
def allUnique(x):
seen = set()
return not any(i in seen or seen.add(i) for i in x)
assert allUnique(Rep_Outer_CV_test_flat_str)
# Get current time
c_time1 = time.localtime()
c_time1 = time.strftime("%a %d %b %Y %H:%M:%S", c_time1)
print(c_time1)
def Each_feat_10rep_10x10CV(rep):
# The outer fold CV has already been saved as lists
Outer_CV = Rep_Outer_CV[rep]
# Pre-allocate memory
Ind_feat_clf0 = []
accuracy_arr0 = np.zeros((k_out,len(EEG_features_name_list),n_models,3))
Outer_counter = 0
for train_index, test_index in Outer_CV:
test_df = EEG_features_df.iloc[test_index]
train_df = EEG_features_df.iloc[train_index]
# Training data will be standardized
standardizer = preprocessing.StandardScaler().fit(train_df.iloc[:,n_discrete_cols:])
train_df_standard = train_df.copy()
train_df_standard.iloc[:,n_discrete_cols:] = standardizer.transform(train_df_standard.iloc[:,n_discrete_cols:])
# Test data will also be standardized but using mean and std from training data
test_df_standard = test_df.copy()
test_df_standard.iloc[:,n_discrete_cols:] = standardizer.transform(test_df_standard.iloc[:,n_discrete_cols:])
# Get the training data
X_train = train_df_standard.copy().drop(Subject_info_cols, axis=1)
y_train = train_df_standard[Target]
# Get test data
X_test = test_df_standard.copy().drop(Subject_info_cols, axis=1)
y_test = test_df_standard[Target]
# Part 2 with second mRMR and classifiers in loop
k_fold = 10
skf2 = StratifiedKFold(n_splits=k_fold, shuffle=True, random_state=(rep+1)*123) # using 10% equals around 21 test subjects
# I am also using converting it to an iterable list instad of the generator to promote reuse
Inner_CV = []
for train_index2, test_index2 in skf2.split(train_df,study_group_status.iloc[train_index]):
Inner_CV.append([train_index2,test_index2])
# Prepare initial filtering of feature types to alleviate imbalance in number of features
# Use a variable that is optimized using inner CV
n_feat_mRMR1 = [20,30,40,50]
# Max features from mRMR for each eye status
max_mRMR_features = n_feat_mRMR1[-1]
eye_status = ["Eyes Closed", "Eyes Open"]
n_eye_status = len(eye_status)
feat_counter = []
# Perform mRMR for each feature type followed by RFE SVM/SFS LogReg/RF
Ind_feat_clf00 = []
for fset in range(len(EEG_features_name_list)):
temp_feat = EEG_features_name_list[fset]
other_feats = EEG_features_name_list[:fset]+EEG_features_name_list[fset+1:]
other_feats = [item for sublist in other_feats for item in sublist] # make the list flatten
# Retrieve the dataset for each feature
col_idx = np.zeros(len(X_train.columns), dtype=bool)
for fsub in range(len(temp_feat)):
temp_feat0 = temp_feat[fsub]
col_idx0 = X_train.columns.str.contains(temp_feat0+"_")
col_idx = np.logical_or(col_idx,col_idx0) # append all trues
temp_X_train = X_train.loc[:,col_idx]
# Check if any of the other features were wrongly chosen
for fcheck in range(len(other_feats)):
if any(temp_X_train.columns.str.contains(other_feats[fcheck]+"_")==True):
temp_X_train = temp_X_train.loc[:,np.invert(temp_X_train.columns.str.contains(other_feats[fcheck]))]
if temp_X_train.size == 0: # if no columns are left, e.g. imcoh when removing coh, then add features again
temp_X_train = X_train.loc[:,X_train.columns.str.contains(temp_feat0+"_")]
# Save number of original features fed into mRMR
feat_counter.append([temp_feat,temp_X_train.shape[1]])
# If there are no features selected, then skip the rest of the loop
if temp_X_train.shape[1] == 0:
continue
# Do not use mRMR if there are fewer than n_features
if temp_X_train.shape[1] <= max_mRMR_features:
filter_features = temp_X_train.columns
else:
# mRMR
filter_feat_selector = mRMR_feature_select(temp_X_train.to_numpy(),y_train.to_numpy(),
num_features_to_select=max_mRMR_features,
K_MAX=1000,n_jobs=1,verbose=False)
# Get selected features
filter_features = temp_X_train.columns[filter_feat_selector]
# SVM with L2-norm (it is by default squared)
# Prepare hyper parameters
exponent = np.linspace(-3,1,9)
exponent = np.round(exponent,5) # sometimes linspace are not so exact
C_parameter_SVM = np.power(np.array([10]*len(exponent)),exponent)
# rbf kernel overfit easier, whereas linear empirically works better in high D data
# Sequential Feature Selection for Logistic Regression
# In-built model selection CV
# The L2 is the inverse of C for LogReg
exponent = np.linspace(-3,1,9)
exponent = np.round(exponent,5) # sometimes linspace are not so exact
C_parameter_LogReg = np.power(np.array([10]*len(exponent)),exponent)
# Random forest classifier
# Prepare hyper parameters
trees = np.array([10, 100, 500, 1000])
depth = np.linspace(1,2,2) # using more depth leads to a lot of overfitting
# Prepare arrays for validation errors
val_err_svm = np.zeros((len(n_feat_mRMR1),len(C_parameter_SVM)))
val_feat_svm = np.zeros((len(n_feat_mRMR1),len(C_parameter_SVM),max_mRMR_features))
val_feat_svm.fill(np.nan)
val_err_logreg = np.zeros((len(n_feat_mRMR1),len(C_parameter_LogReg)))
val_feat_logreg = np.zeros((len(n_feat_mRMR1),len(C_parameter_LogReg),max_mRMR_features))
val_feat_logreg.fill(np.nan)
val_err_rf = np.zeros((len(n_feat_mRMR1),len(trees),len(depth)))
val_feat_rf = np.zeros((len(n_feat_mRMR1),len(trees),len(depth),max_mRMR_features))
val_feat_rf.fill(np.nan)
val_feat_import_rf = np.zeros((len(n_feat_mRMR1),len(trees),len(depth),max_mRMR_features))
val_feat_import_rf.fill(np.nan)
min_features_to_select = 1
for n1 in range(len(n_feat_mRMR1)):
mRMR_features1 = filter_features[0:n_feat_mRMR1[n1]]
# Use the selected filter features
temp_X_train_mRMR = temp_X_train[mRMR_features1]
# SVM with recursive feature elemination
# Stratified CV to find best regularization strength and number of features
# CV is built-in RFECV
for C in range(len(C_parameter_SVM)):
# Define the model
svc = SVC(C=C_parameter_SVM[C], kernel="linear", tol=1e-3, cache_size=4000)
# Perform recurive feature elimination with in-built CV
rfecv = RFECV(estimator=svc, n_jobs=1, scoring="balanced_accuracy",
cv=Inner_CV,
min_features_to_select=min_features_to_select)
rfecv.fit(temp_X_train_mRMR,y_train.to_numpy().ravel())
# Get CV score
err = rfecv.grid_scores_[rfecv.n_features_-min_features_to_select]
# Save results for hyperparameter optimization
val_err_svm[n1,C] = err
# Save the features chosen by RFE based on index from pre mRMR
rfe_feat_idx = np.where(X_train.columns.isin(mRMR_features1[rfecv.support_]))[0]
val_feat_svm[n1,C,0:len(rfe_feat_idx)] = rfe_feat_idx
# print("Finished SVM run {} out of {}".format(C+1,len(C_parameter_SVM)))
# Logistic regression with sequential forward selection
# In-bult CV
k_max = np.min([temp_X_train_mRMR.shape[1],n_feat_mRMR1[n1]])
k_features = np.arange(1,k_max+1,1) # try up to all features
inner_cv_scores = np.zeros((len(C_parameter_LogReg),len(k_features)))
for C in range(len(C_parameter_LogReg)):
LogReg = LogisticRegression(penalty="l2", C=C_parameter_LogReg[C],
max_iter = 50000)
sfs = SFS(LogReg, k_features = (k_features[0],k_features[-1]),
forward = True, scoring = "balanced_accuracy",
verbose = 0, floating = False,
cv = Inner_CV, n_jobs = 1)
sfs = sfs.fit(temp_X_train_mRMR, y_train.to_numpy().ravel())
# Save CV scores for each SFS step
for feat in range(len(k_features)):
inner_cv_scores[C,feat] = sfs.get_metric_dict()[k_features[feat]]["avg_score"]
# Find the best number of features
# I am rounding to make it more smooth and disregard small improvements
K = np.where(np.round(inner_cv_scores[C,:],2)==np.max(np.round(inner_cv_scores[C,:],2)))[0]
if len(K) > 1:
K = K[np.where(K == np.min(K))[0]]
err = inner_cv_scores[C,K]
# Save validation error and features
val_err_logreg[n1,C] = err
# Save the features chosen by RFE based on index from pre mRMR
sfs_feat_idx = np.where(X_train.columns.isin(sfs.subsets_[int(K)+1]["feature_names"]))[0] # K+1 since it starts with 1 feature
val_feat_logreg[n1,C,0:len(sfs_feat_idx)] = sfs_feat_idx
# print("Finished LogReg run {} out of {}".format(C+1,len(C_parameter_LogReg)))
# Random forest
model_test_err = np.zeros((k_fold,len(trees),len(depth)))
RF_feat_import = np.zeros((k_fold,len(trees),len(depth),temp_X_train_mRMR.shape[1]))
counter = 0
for cv in range(k_fold):
# Retrieve CV indices
train_index_rf = Inner_CV[cv][0]
test_index_rf = Inner_CV[cv][1]
# Retrieve datasets
X_train2 = temp_X_train_mRMR.iloc[train_index_rf]; X_test2 = temp_X_train_mRMR.iloc[test_index_rf]
y_train2 = y_train.to_numpy().ravel()[train_index_rf]; y_test2 = y_train.to_numpy().ravel()[test_index_rf]
for t in range(len(trees)):
for d in range(len(depth)):
RF = RandomForestClassifier(n_estimators=trees[t], max_features="sqrt", max_depth=depth[d],
n_jobs=1, random_state=None)
RF.fit(X_train2.to_numpy(),y_train2)
# Validation error
RF_y = RF.predict(X_test2.to_numpy())
err = balanced_accuracy_score(y_test2,RF_y)
# Save to array
model_test_err[cv,t,d] = err
# Save feature importance array
RF_feat_import[cv,t,d] = RF.feature_importances_
counter += 1
# print("Finished RF run {} out of {}".format(counter, k_fold))
# Average the errors over the CV folds
model_test_err_mean = np.mean(model_test_err, axis=0)
val_err_rf[n1,:,:] = model_test_err_mean
# Average the feature importances over the CV folds
RF_feat_import = np.mean(RF_feat_import, axis=0)
val_feat_import_rf[n1,:,:,0:RF_feat_import.shape[2]] = RF_feat_import
# Save the features used by the RF based on index from pre mRMR
rf_feat_idx = np.where(X_train.columns.isin(mRMR_features1))[0]
val_feat_rf[n1,:,:,0:len(rf_feat_idx)] = rf_feat_idx
# Print progress
current_progress = (n1+1)/(len(n_feat_mRMR1))*100
print("Finished {}% of inner fold optimization for feat: {}".format(current_progress,temp_feat[0]))
# Choose the optimal parameters
### SVM
n1, C = np.where(val_err_svm==np.max(val_err_svm))
if len(C) > 1:
print("There are multiple SVM runs with the same validation error")
print("Choosing run with fewest features to alleviate overfitting")
rfe_feat_len = []
for i2 in range(len(C)):
n1_temp = int(n1[i2])
C_temp = int(C[i2])
temp_feats = val_feat_svm[n1_temp,C_temp][~np.isnan(val_feat_svm[n1_temp,C_temp])].astype(int)
rfe_feat_len.append(len(temp_feats))
rfe_feat_min = np.where(rfe_feat_len==np.min(rfe_feat_len))[0]
if len(rfe_feat_min) > 1:
print("Multiple SVM runs with same number of fewest features")
print("Choosing run with lowest C (highest regularization) to alleviate overfitting")
C_min = np.argmin(C[rfe_feat_min])
n1, C = [int(n1[rfe_feat_min][C_min]),
int(C[rfe_feat_min][C_min])]
else:
n1, C = [int(n1[int(rfe_feat_min)]),
int(C[int(rfe_feat_min)])]
else:
n1, C = [int(n1), int(C)]
mRMR_chosen1 = n_feat_mRMR1[n1]
C_chosen = C_parameter_SVM[C]
# Save the best validation error
val_error = val_err_svm[n1, C]
# Get the subsets chosen
chosen_feats = val_feat_svm[n1,C][~np.isnan(val_feat_svm[n1,C])].astype(int)
rfe_features = list(X_train.columns[chosen_feats])
n_final_feat = len(rfe_features)
x_train_mRMR2_rfe = X_train[rfe_features]
# Fit on all training data
model = SVC(C=C_chosen, kernel="linear", tol=1e-3, cache_size=4000)
model = model.fit(x_train_mRMR2_rfe,y_train.to_numpy().ravel())
# Get training output
model_train_y = model.predict(x_train_mRMR2_rfe)
# Get training error
train_error = balanced_accuracy_score(y_train, model_train_y)
# Get prediction of test data
y_pred = model.predict(X_test[rfe_features])
# Use model to predict class on test data
test_error = balanced_accuracy_score(y_test, y_pred)
# Save or prepare to save
accuracy_arr0[Outer_counter,fset,0,:] = [train_error,val_error,test_error]
SVM_model_par = [mRMR_chosen1,C_chosen,n_final_feat]
SVM_model = model
SVM_y_pred = [[model_train_y],[y_pred]]
### LogReg with SFS
n1, C = np.where(val_err_logreg==np.max(val_err_logreg))
if len(C) > 1:
print("There are multiple LogReg runs with the same validation error")
print("Choosing run with fewest features to alleviate overfitting")
sfs_feat_len = []
for i2 in range(len(C)):
n1_temp = int(n1[i2])
C_temp = int(C[i2])
temp_feats = val_feat_logreg[n1_temp,C_temp][~np.isnan(val_feat_logreg[n1_temp,C_temp])].astype(int)
sfs_feat_len.append(len(temp_feats))
sfs_feat_min = np.where(sfs_feat_len==np.min(sfs_feat_len))[0]
if len(sfs_feat_min) > 1:
print("Multiple LogReg runs with same number of fewest features")
print("Choosing run with lowest C (highest regularization) to alleviate overfitting")
C_min = np.argmin(C[sfs_feat_min])
n1, C = [int(n1[sfs_feat_min][C_min]),
int(C[sfs_feat_min][C_min])]
else:
n1, C = [int(n1[int(sfs_feat_min)]),
int(C[int(sfs_feat_min)])]
else:
n1, C = [int(n1), int(C)]
mRMR_chosen1 = n_feat_mRMR1[n1]
C_chosen = C_parameter_LogReg[C]
# Save the best validation erro
val_error = val_err_logreg[n1, C]
# Get the subsets chosen
chosen_feats = val_feat_logreg[n1,C][~np.isnan(val_feat_logreg[n1,C])].astype(int)
sfs_features = list(X_train.columns[chosen_feats])
n_final_feat = len(sfs_features)
x_train_mRMR2_sfs = X_train[sfs_features]
# Fit on all training data
model = LogisticRegression(penalty="l2", C=C_chosen, max_iter = 50000)
model = model.fit(x_train_mRMR2_sfs,y_train.to_numpy().ravel())
# Get training output
model_train_y = model.predict(x_train_mRMR2_sfs)
# Get training error
train_error = balanced_accuracy_score(y_train, model_train_y)
# Get prediction of test data
y_pred = model.predict(X_test[sfs_features])
# Use model to predict class on test data
test_error = balanced_accuracy_score(y_test, y_pred)
# Save or prepare to save
accuracy_arr0[Outer_counter,fset,1,:] = [train_error,val_error,test_error]
LogReg_model_par = [mRMR_chosen1,C_chosen,n_final_feat]
LogReg_model = model
LogReg_y_pred = [[model_train_y],[y_pred]]
### RF
n1, t, d = np.where(val_err_rf==np.max(val_err_rf))
if len(d) > 1:
print("There are multiple RF runs with the same validation error")
print("Choosing run with lowest depth to alleviate overfitting")
d_min = np.where(d==np.min(d))[0]
if len(d_min) > 1:
print("Multiple RF runs with same number number of depths")
print("Choosing run with lowest trees to alleviate overfitting")
t_min = np.argmin(t[d_min]) # argmin just takes the first
# If there are multiple with same parameters and validation error it is most likely the same
n1, t, d = [int(n1[d_min][t_min]),
int(t[d_min][t_min]), int(d[d_min][t_min])]
else:
n1, t, d = [int(n1[int(d_min)]),
int(t[int(d_min)]), int(d[int(d_min)])]
else:
n1, t, d = [int(n1), int(t), int(d)]
mRMR_chosen1 = n_feat_mRMR1[n1]
t_chosen = trees[t]
d_chosen = depth[d]
# Save the best validation error
val_error = val_err_rf[n1, t, d]
# Get the chosen features and feature importances
chosen_feats = val_feat_rf[n1,t,d][~np.isnan(val_feat_rf[n1,t,d])].astype(int)
rf_features = list(X_train.columns[chosen_feats])
n_final_feat = len(rf_features)
x_train_mRMR2_rf = X_train[rf_features]
rf_feat_importances = val_feat_import_rf[n1,t,d][~np.isnan(val_feat_import_rf[n1,t,d])]
# Fit on all training data
model = RandomForestClassifier(n_estimators=t_chosen, max_features="sqrt", max_depth=d_chosen,
n_jobs=1, random_state=None)
model = model.fit(x_train_mRMR2_rf,y_train.to_numpy().ravel())
# Get training output
model_train_y = model.predict(x_train_mRMR2_rf)
# Get training error
train_error = balanced_accuracy_score(y_train, model_train_y)
# Get prediction of test data
y_pred = model.predict(X_test[rf_features])
# Use model to predict class on test data
test_error = balanced_accuracy_score(y_test, y_pred)
# Save or prepare to save
accuracy_arr0[Outer_counter,fset,2,:] = [train_error,val_error,test_error]
RF_model_par = [mRMR_chosen1,t_chosen,d_chosen]
RF_model = model
RF_y_pred = [[model_train_y],[y_pred]]
# Save the rest
model_par00 = [SVM_model_par,LogReg_model_par,RF_model_par]
final_models00 = [SVM_model,LogReg_model,RF_model]
final_features00 = [rfe_features, sfs_features, [rf_features,rf_feat_importances]]
final_y_preds00 = [SVM_y_pred,LogReg_y_pred,RF_y_pred]
data_splits00 = [X_train,y_train,X_test,y_test,standardizer.mean_,standardizer.scale_]
# Save the results for the specific feature
res = [temp_feat, model_par00, final_models00, final_features00, final_y_preds00, data_splits00]
Ind_feat_clf00.append(res)
print("Finished outer fold {} out of {} for rep: {} for feat {}".format(Outer_counter+1,k_out,rep+1,temp_feat))
# Check that no features were excluded when checking for wrongly chosen
total_features = 0
for feat_c in range(len(feat_counter)):
total_features += feat_counter[feat_c][1]
assert total_features == X_train.shape[1]
# Save results over outer folds
Ind_feat_clf0.append(Ind_feat_clf00)
print("Finished outer fold {} out of {} for rep: {} for all features".format(Outer_counter+1,k_out,rep+1))
# Move counter
Outer_counter += 1
out = [accuracy_arr0, Ind_feat_clf0] # [Rep][Outer CV][Feature][Variable]
# Print total progress
print("Finished outer fold repetition {} out of {}".format(rep+1,n_repetitions))
# Get current time
c_time2 = time.localtime()
c_time2 = time.strftime("%a %d %b %Y %H:%M:%S", c_time2)
print("Started", c_time1, "\nCurrent Time",c_time2)
return rep, out
Ind_feat_clf = [0]*n_repetitions
with concurrent.futures.ProcessPoolExecutor(max_workers=10) as executor:
for rep, result in executor.map(Each_feat_10rep_10x10CV, range(n_repetitions)): # Function and arguments
Ind_feat_clf[rep] = result # [acc or parameters], for parameters: [Rep][Outer CV][Feature][Variable]
# Save results to file and overwrite when new results arrive
with open(Model_savepath+"Each_feat_Rep10_10x10CV_mRMR2_SVM_LogReg_RF_Subtype1_results_250122.pkl", "wb") as filehandle:
pickle.dump(Ind_feat_clf, filehandle)
# %% For subtype 2
# The two previous sections are rerun but for subtype 2
# The only modification was indexing the subjects
# Only keep PTSD subtype 2 by dropping subtype 1
subtype1 = np.array(Subject_id)[PTSD_idx][sparcl["cs"]==0]
subtype2 = np.array(Subject_id)[PTSD_idx][sparcl["cs"]==1]
EEG_features_df = EEG_features_df.set_index("Subject_ID")
EEG_features_df = EEG_features_df.drop(subtype1)
EEG_features_df = EEG_features_df.reset_index()
# Check it was dropped correctly
assert all(subtype2 == EEG_features_df.loc[EEG_features_df["Group_status"] == 1,"Subject_ID"])