Machine learning-based identification of abnormal functional connectivity in obesity across different metabolic states.

Code ↔ Paper

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

The 5 matches

[1] § Methods › Feature selection and classification ↔ CodeScript.ipynb, lines 403–515 · score 0.74 · max_depth, n_estimators, feature selection, optimal, iterations, folds
[2] § Methods › Statistics and reproducibility ↔ CodeScript.ipynb, lines 796–825 · score 0.71 · Mann Whitney, Shapiro Wilk, accuracy
[3] § Methods › Feature selection and classification ↔ CodeScript.ipynb, lines 403–515 · score 0.64 · feature selection, consecutive, threshold, iteratively, score, RF
[4] § Methods › Feature selection and classification ↔ CodeScript.ipynb, lines 618–707 · score 0.58 · F1 score, AUC, recall, precision, metric, class
[5] § Results › Machine learning model performance ↔ CodeScript.ipynb, lines 618–707 · score 0.51 · F1 score, AUC, recall, precision, metrics, class

Paper

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

The paper is loaded when this pane is shown.

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

The authors' code

Jupyter notebook · 1,336 lines · 50 KB · CC-BY-4.0 · 5 matches

# %%
import librosa
import os as os
import pandas as pd
import re
import numpy as np
from sklearn.metrics import precision_score, recall_score, accuracy_score
from matplotlib import cm, colors, colorbar
from matplotlib import pyplot as plt
from sklearn.neural_network import MLPClassifier
from sklearn.linear_model import RidgeClassifier
rdg = RidgeClassifier(alpha=0.5)
#mlp=MLPClassifier(random_state=1,max_iter=300,activation='relu',solver='sgd',learning_rate='constant',learning_rate_init=0.0001)
mlp=MLPClassifier(random_state=1,max_iter=300,activation='relu')
from sklearn.linear_model import LogisticRegression
lgr=LogisticRegression(random_state=1,max_iter=500)
from sklearn.tree import DecisionTreeClassifier
DT = DecisionTreeClassifier(random_state=0,max_depth=10)
from sklearn.ensemble import AdaBoostClassifier
adb = AdaBoostClassifier(n_estimators=100, random_state=0)
from sklearn.ensemble import GradientBoostingClassifier
gbc= GradientBoostingClassifier(n_estimators=100, random_state=1)
from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=3)
from sklearn.linear_model import SGDClassifier
SGD=SGDClassifier(loss= 'log',random_state=1,max_iter=100,early_stopping=True,learning_rate='optimal',validation_fraction=0.2)
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import MinMaxScaler
scaler = StandardScaler()
mmscaler= MinMaxScaler()
from sklearn.ensemble import RandomForestClassifier
rf = RandomForestClassifier(n_estimators=75,max_depth=15, random_state=0)
from sklearn.svm import SVC
clf_svm=SVC(kernel='rbf')
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
from scipy.spatial import ConvexHull, convex_hull_plot_2d
from sklearn.cluster import DBSCAN
from sklearn.decomposition import PCA
from scipy.stats.mstats import mquantiles
from scipy.stats import skew
from sklearn.cluster import KMeans
from sklearn.model_selection import LeaveOneOut
pca = PCA(n_components=2, svd_solver='full')
import random
from sklearn.cluster import DBSCAN
from sklearn.model_selection import train_test_split
from sklearn.model_selection import cross_val_score
#import nilearn
#from nilearn import plotting
from matplotlib.pyplot import figure
import seaborn as sns
import statistics
import warnings
warnings.filterwarnings('ignore')
from sklearn.metrics import (recall_score, confusion_matrix, balanced_accuracy_score,
f1_score, precision_score, roc_auc_score,
matthews_corrcoef, average_precision_score)
# %% [markdown]
# ### Read Functional Connectivity Data
# %%
def extract_per_time_Lean(time,directory):
coh={}
shorten_list={}
for time,directory,in zip([time],[directory]):
os.chdir(directory)
coh[time]=os.listdir(directory)
if '.ipynb_checkpoints' in coh[time]:
coh[time].remove('.ipynb_checkpoints')
shorten_list[time]=[]
for i in range(0,len(coh[time])):
if coh[time][i][0]=='L':
shorten_list[time].append(coh[time][i].split(str(time)+'min')[0])
ns=(set(shorten_list[time])
)
ns=list(ns)
ns.sort()
#random.shuffle(ns) ### Shuffle the list
ns_l=[]
ns_ob=[]
for n in ns:
if n[0]=='L':
ns_l.append(n)
else:
ns_ob.append(n)
names_l=[]
for i in ns_l:
name=i+'0min.EDFOUT-ROIlaggedCoh-covar.txt'
names_l.append(name)
#names_l=names_l[0:27]
l_name_ls={}
for t in [time]:
sr=str(t)+'min'
print(sr)
l_name_ls[t]=[]
for n in names_l:
nn=n.replace('0min',sr)
l_name_ls[t].append(nn)
names_list={}
for t in [time]:
names_list[t]=l_name_ls[t]
names_list[t]=names_list[t][0:30]
def extract_connectivity(band,data):
Y=[]
coh_ar=np.zeros([len(data),88*88])
for i in range(0,len(data)):
m=np.loadtxt(data[i])[band*88:(band+1)*88,:] # extract only delta band
m=np.tril(m, k=-1).flatten() ## Take upper/lower Triangle of the Symetrical Coherence Matrix
coh_ar[i,:]=m
#cor_ar=cor_ar[0:60,:]
if (data[i][0])=='L':
Y.append(0)
else:
Y.append(1)
#Y=Y[0:60]
return coh_ar,Y
connectivity={}
for band in [0,1,2,3,4]: # 0-delta, 1-theta, 2-alpha, 3-beta, 4-gamma
#for band in [2]:
#print(band)
connectivity[band]=np.zeros([1,88*88])
Y=[]
for time,directory,in zip([time],[directory]):
os.chdir(directory)
data=names_list[time][0:60]
con=extract_connectivity(band,data)[0]
y=extract_connectivity(band,data)[1]
connectivity[band]=np.vstack([connectivity[band],con])
Y=Y+y
connectivity[band]=connectivity[band][1:,:]
print(connectivity[band].shape)
print(' ')
con_all_bl_lean_x=np.hstack([connectivity[0],connectivity[1],connectivity[2],connectivity[3],connectivity[4]])
#con_alpha_bl_lean_x=connectivity[band]
con_all_bl_lean_y=np.zeros([con_all_bl_lean_x.shape[0]])
return con_all_bl_lean_x,con_all_bl_lean_y
# %%
def extract_per_time_Ob(time,directory):
coh={}
shorten_list={}
for time,directory,in zip([time],[directory]):
os.chdir(directory)
coh[time]=os.listdir(directory)
if '.ipynb_checkpoints' in coh[time]:
coh[time].remove('.ipynb_checkpoints')
shorten_list[time]=[]
for i in range(0,len(coh[time])):
shorten_list[time].append(coh[time][i].split(str(time)+'min')[0])
# find those files that belong to the subjects that are not missing in any timestates
ns=(set(shorten_list[time])
)
ns=list(ns)
ns.sort()
#random.shuffle(ns) ### Shuffle the list
ns_l=[]
ns_ob=[]
for n in ns:
if n[0]=='L':
ns_l.append(n)
else:
ns_ob.append(n)
names_ob=[]
for r in ns_ob:
name=r+'0min-ROIlaggedCoh-covar.txt'
names_ob.append(name)
ob_name_ls={}
for t in [time]:
sr=str(t)+'min'
print(sr)
ob_name_ls[t]=[]
for n in names_ob:
nn=n.replace('0min',sr)
ob_name_ls[t].append(nn)
names_list={}
for t in [time]:
names_list[t]=ob_name_ls[t]
names_list[t]=names_list[t][0:30]
def extract_connectivity(band,data):
Y=[]
coh_ar=np.zeros([len(data),88*88])
for i in range(0,len(data)):
m=np.loadtxt(data[i])[band*88:(band+1)*88,:] # extract only delta band
m=np.tril(m, k=-1).flatten() ## Take upper/lower Triangle of the Symetrical Coherence Matrix
coh_ar[i,:]=m
#cor_ar=cor_ar[0:60,:]
if (data[i][0])=='L':
Y.append(0)
else:
Y.append(1)
#Y=Y[0:60]
return coh_ar,Y
connectivity={}
for band in [0,1,2,3,4]: # 0-delta, 1-theta, 2-alpha, 3-beta, 4-gamma
#for band in [2]:
#print(band)
connectivity[band]=np.zeros([1,88*88])
Y=[]
for time,directory,in zip([time], [directory]):
os.chdir(directory)
data=names_list[time][0:60]
con=extract_connectivity(band,data)[0]
y=extract_connectivity(band,data)[1]
connectivity[band]=np.vstack([connectivity[band],con])
Y=Y+y
connectivity[band]=connectivity[band][1:,:]
#print(connectivity[band].shape)
#print(' ')
con_all_pwl_ob_x=np.hstack([connectivity[0],connectivity[1],connectivity[2],connectivity[3],connectivity[4]])
con_all_pwl_ob_y=np.ones([con_all_pwl_ob_x.shape[0]])
print(con_all_pwl_ob_x.shape,con_all_pwl_ob_y.shape)
return con_all_pwl_ob_x,con_all_pwl_ob_y
# %%
Xdata={}
Ydata={}
for time in [0,15,30,45,60,90,120,180,240]:
xarray=[]
yarray=[]
for directory in [f"/datasets/Demo_MOCKED_data/3MON/T{time}" ,
f"/datasets/Demo_MOCKED_data/PostWL/T{time}"
]:
#for directory in [f"/datasets/MOCKED_data/3MON/T{time}" ,
#
# f"/datasets/MOCKED_data/PostWL/T{time}"
# ]:
xarray.append(extract_per_time_Ob(time,directory)[0])
yarray.append(extract_per_time_Ob(time,directory)[1])
xarray=np.concatenate(xarray, axis=0)
yarray=np.concatenate(yarray, axis=0)
Xdata[time]=xarray
Ydata[time]=yarray
# %%
for time in [0,15,30,45,60,90,120,180,240]:
print(Xdata[time].shape, Ydata[time].shape)
# %%
def extract_per_time_Ob_BL(time,directory):
coh={}
shorten_list={}
for time,directory,in zip([time],[directory]):
os.chdir(directory)
coh[time]=os.listdir(directory)
if '.ipynb_checkpoints' in coh[time]:
coh[time].remove('.ipynb_checkpoints')
shorten_list[time]=[]
for i in range(0,len(coh[time])):
shorten_list[time].append(coh[time][i].split(str(time)+'min')[0])
# find those files that belong to the subjects that are not missing in any timestates
ns=(set(shorten_list[time])
)
ns=list(ns)
ns.sort()
#random.shuffle(ns) ### Shuffle the list
ns_l=[]
ns_ob=[]
for n in ns:
if n[0]=='L':
ns_l.append(n)
else:
ns_ob.append(n)
names_ob=[]
for r in ns_ob:
name=r+'0min.EDFOUT-ROIlaggedCoh-covar.txt'
names_ob.append(name)
ob_name_ls={}
for t in [time]:
sr=str(t)+'min'
print(sr)
ob_name_ls[t]=[]
for n in names_ob:
nn=n.replace('0min',sr)
ob_name_ls[t].append(nn)
names_list={}
for t in [time]:
names_list[t]=ob_name_ls[t]
names_list[t]=names_list[t][0:30]
def extract_connectivity(band,data):
Y=[]
coh_ar=np.zeros([len(data),88*88])
for i in range(0,len(data)):
m=np.loadtxt(data[i])[band*88:(band+1)*88,:] # extract only delta band
m=np.tril(m, k=-1).flatten() ## Take upper/lower Triangle of the Symetrical Coherence Matrix
coh_ar[i,:]=m
#cor_ar=cor_ar[0:60,:]
if (data[i][0])=='L':
Y.append(0)
else:
Y.append(1)
#Y=Y[0:60]
return coh_ar,Y
connectivity={}
for band in [0,1,2,3,4]: # 0-delta, 1-theta, 2-alpha, 3-beta, 4-gamma
#for band in [2]:
#print(band)
connectivity[band]=np.zeros([1,88*88])
Y=[]
for time,directory,in zip([time], [directory]):
os.chdir(directory)
data=names_list[time][0:60]
con=extract_connectivity(band,data)[0]
y=extract_connectivity(band,data)[1]
connectivity[band]=np.vstack([connectivity[band],con])
Y=Y+y
connectivity[band]=connectivity[band][1:,:]
#print(connectivity[band].shape)
#print(' ')
con_all_pwl_ob_x=np.hstack([connectivity[0],connectivity[1],connectivity[2],connectivity[3],connectivity[4]])
con_all_pwl_ob_y=np.ones([con_all_pwl_ob_x.shape[0]])
print(con_all_pwl_ob_x.shape,con_all_pwl_ob_y.shape)
return con_all_pwl_ob_x,con_all_pwl_ob_y
# %%
Xdata_={}
Ydata_={}
for time in [0,15,30,45,60,90,120,180,240]:
xarray=[]
yarray=[]
for directory in [f"/datasets/Demo_MOCKED_data/BL/T{time}"]:
#for directory in [f"/datasets/MOCKED_data/BL/T{time}"]:
xarray.append(extract_per_time_Ob_BL(time,directory)[0])
yarray.append(extract_per_time_Ob_BL(time,directory)[1])
xarray=np.concatenate(xarray, axis=0)
yarray=np.concatenate(yarray, axis=0)
Xdata_[time]=xarray
Ydata_[time]=yarray
# %%
Xdata_ob={}
Ydata_ob={}
for time in [0,15,30,45,60,90,120,180,240]:
Xdata_ob[time]=np.vstack([Xdata_[time],Xdata[time]])
Ydata_ob[time]=np.hstack([Ydata_[time],Ydata[time]])
# %%
Xdata={}
Ydata={}
for time in [0,15,30,45,60,90,120,180,240]:
xarray=[]
yarray=[]
for directory in [f"/datasets/Demo_MOCKED_data/BL/T{time}"]:
xarray.append(extract_per_time_Lean(time,directory)[0])
yarray.append(extract_per_time_Lean(time,directory)[1])
xarray=np.concatenate(xarray, axis=0)
yarray=np.concatenate(yarray, axis=0)
Xdata[time]=xarray
Ydata[time]=yarray
# %%
Xdata_Lean={}
Ydata_Lean={}
for time in [0,15,30,45,60,90,120,180,240]:
Xdata_Lean[time]=Xdata[time]
Ydata_Lean[time]=Ydata[time]
# %%
Xdata_Lean_list = list(Xdata_Lean.values())
Xdata_ob_list = list(Xdata_ob.values())
Xdata_Lean_concatenated = np.concatenate(Xdata_Lean_list, axis=0)
Xdata_ob_concatenated = np.concatenate(Xdata_ob_list, axis=0)
print(Xdata_Lean_concatenated.shape)
print(Xdata_ob_concatenated.shape)
xall=np.vstack([Xdata_ob_concatenated ,Xdata_Lean_concatenated ])
yall=np.array([1]*Xdata_ob_concatenated.shape[0]+[0]*Xdata_Lean_concatenated.shape[0])
# %%
def select_top_n_feature(data,numb_feat,id_ls,y,rf_estimators,rf_depth):
dX=id_ls
Top1_feat_id=[]
rf=RandomForestClassifier(n_estimators=rf_estimators,max_depth=rf_depth, random_state=0,class_weight='balanced')
rf.fit(data[:,dX],y)
Top1_feat_id=(-rf.feature_importances_).argsort()[:numb_feat].tolist()
return (Top1_feat_id)
def kfoldCV(xdata,ydata):
cv = StratifiedKFold(n_splits=10, shuffle=True, random_state=0)
classifier = XGBClassifier(n_estimators=xgb_est, learning_rate=xgb_lr, max_depth=xgb_depth,
subsample=xgb_subsamp, colsample_bytree=xgb_colsam, reg_alpha=0.5,
random_state=42, use_label_encoder=False, eval_metric='logloss')
acc = cross_val_score(classifier, xdata, ydata, scoring='balanced_accuracy', cv=cv, n_jobs=-1)
avg_acc=sum(acc/acc.shape[0])
return (avg_acc )
def feature_selection(traindataX, traindataY,
xgb_est,xgb_lr,xgb_depth,xgb_subsamp,xgb_colsam,
rf_estimators, rf_depth, meric_Threshold,numb_runs):
existing_id_ls = []
acc_dif = {}
DF = {}
rf_filtered_id_ls = {}
Xdata_train = traindataX
Ydata_train = traindataY
all_fidls = np.arange(0, Xdata_train.shape[1])
# Select the top features based on initial selection
rf_filtered_id_ls = select_top_n_feature(Xdata_train, Xdata_train.shape[1], all_fidls, traindataY, rf_estimators, rf_depth)
# Initialize the existing_id_ls with the best feature
existing_id_ls = select_top_n_feature(Xdata_train, 1, np.arange(0, Xdata_train.shape[1]), Ydata_train, rf_estimators, rf_depth)
# Get the rest of the features that are not in existing_id_ls
rest_id_ls = [e for e in rf_filtered_id_ls if e not in existing_id_ls]
mrf_index = 0
consecutive_runs = 0
highest_acc = -np.inf # To track the highest accuracy
best_id_ls = existing_id_ls.copy() # To track the best feature set
while mrf_index < len(rest_id_ls) and mrf_index < numb_runs:
print('Run', mrf_index)
# Select the next feature to add
mrf = rest_id_ls[mrf_index]
select_id = existing_id_ls + [mrf]
# Evaluate the accuracy with the new feature added
xdata = Xdata_train[:, select_id]
ydata = np.array(Ydata_train)
acc_newf = kfoldCV(xdata, ydata)
#acc_newf = kfold_cv(xdata, ydata, valdataX[:, select_id], valdataY,xgb_est,xgb_lr,xgb_depth,xgb_subsamp,xgb_colsam)
# Evaluate the accuracy with the current set of selected features
xdata1 = Xdata_train[:, existing_id_ls]
ydata1 = np.array(Ydata_train)
acc_currentselect = kfoldCV(xdata1, ydata1)
#acc_currentselect = kfold_cv(xdata1, ydata1, valdataX[:, existing_id_ls], valdataY,xgb_est,xgb_lr,xgb_depth,xgb_subsamp,xgb_colsam)
print(acc_newf, select_id, acc_currentselect, existing_id_ls, 'p_gain:', acc_newf - acc_currentselect)
# Track the highest accuracy and the corresponding feature set
if acc_newf - acc_currentselect >meric_Threshold: # <-- Only update if the new accuracy is better than current
existing_id_ls = select_id.copy() # Update the current feature set
consecutive_runs = 0 # Reset the counter if acc_newf > acc_currentselect
# Update the best set if this is the highest accuracy seen
if acc_newf > highest_acc:
highest_acc = acc_newf
best_id_ls = select_id.copy() # Update the best feature set
else:
consecutive_runs += 1
# Feature reselection process (if new feature does not improve performance)
rf = RandomForestClassifier(n_estimators=rf_estimators, max_depth=rf_depth, class_weight='balanced', random_state=42)
rf.fit(Xdata_train[:, select_id], Ydata_train)
imp_score = rf.feature_importances_
# Remove the least important feature
toremove_id = select_id[np.argmin(imp_score)]
existing_id_ls = select_id.copy() # Always use a copy to prevent mutation
existing_id_ls.remove(toremove_id)
# Stopping criteria
if consecutive_runs == numb_runs:
break
if len(existing_id_ls) > 30:
break
mrf_index += 1
print('consecutive_runs', consecutive_runs)
print('existing_id_ls', existing_id_ls)
# Return the feature set that gave the highest accuracy across all iterations
return best_id_ls
#def get_best_n_features(feature_list, traindataX, traindataY,xval,yval,rf_estimators, rf_depth):
def get_best_n_features(feature_list, traindataX, traindataY,rf_estimators, rf_depth):
selected_features = []
accuracy_per_num_features = []
skf = StratifiedKFold(n_splits=9, shuffle=True, random_state=42)
for i, feature in enumerate(feature_list):
#print(f"Adding feature {i+1}/{len(feature_list)}: Feature {feature}")
selected_features.append(feature)
accuracies=kfoldCV(traindataX[:, selected_features],np.array(traindataY))
accuracy_per_num_features.append(accuracies)
optimal_num_features = np.argmax(accuracy_per_num_features) + 1 # Adding 1 since index starts at 0
print(f"\nOptimal number of features: {optimal_num_features} with accuracy: {accuracy_per_num_features[optimal_num_features-1]}")
best_acc=accuracy_per_num_features[optimal_num_features-1]
return optimal_num_features, best_acc,accuracy_per_num_features
# %% [markdown]
# ### Traning Set and Testing Set Data Split
# %%
from sklearn.model_selection import StratifiedKFold
from imblearn.over_sampling import RandomOverSampler
import numpy as np
# Initialize the outer StratifiedKFold
skf = StratifiedKFold(n_splits=10, shuffle=True, random_state=42)
# Initialize dictionaries for storing training, validation, and test sets
Xtv = {}
Ytv = {}
Xtest = {}
Ytest = {}
ros = RandomOverSampler(random_state=42)
for time, i in zip([0, 15, 30, 45, 60, 90, 120, 180, 240], range(9)):
X = np.vstack([Xdata_ob[time], Xdata_Lean[time]])
y = np.hstack([Ydata_ob[time], Ydata_Lean[time]])
Xtv[time] = {}
Ytv[time] = {}
Xtest[time] = {}
Ytest[time] = {}
# Initialize fold count
fold = 0
for train_index, test_index in skf.split(X, y):
print('fold', fold)
fold += 1
# Create outer training and testing sets
X_train, y_train = X[train_index], y[train_index]
print('X_train, y_train',X_train.shape,y_train.shape,sum(y_train==0),sum(y_train==1))
X_resampled, y_resampled = ros.fit_resample(X_train, y_train)
print('X_resampled, y_resampled',X_resampled.shape,y_resampled.shape,sum(y_resampled==0),sum(y_resampled==1))
X_test, y_test = X[test_index], y[test_index]
# Store the inner training, validation, and outer test sets
Xtv[time][fold] = X_resampled#X_tv
Ytv[time][fold] = y_resampled#Y_tv
Xtest[time][fold] = X_test
Ytest[time][fold] = y_test
print(' ')
# %%
#os.chdir('/home/jupy/SourceAnalysis/Journal2024_updatedMICMAC/data_Stratifysplit_on_TIME')
Xtv_alltime={}
Xtest_alltime={}
#Xval_alltime={}
Ytv_alltime={}
Ytest_alltime={}
#Yval_alltime={}
for fold in range(1,11):
print(fold)
Xtv_alltime[fold]=np.vstack([ Xtv[0][fold], Xtv[15][fold], Xtv[30][fold], Xtv[45][fold], Xtv[60][fold],
Xtv[90][fold], Xtv[120][fold],Xtv[180][fold], Xtv[240][fold] ])
Xtest_alltime[fold]=np.vstack([ Xtest[0][fold], Xtest[15][fold], Xtest[30][fold], Xtest[45][fold], Xtest[60][fold],
Xtest[90][fold], Xtest[120][fold],Xtest[180][fold], Xtest[240][fold] ])
Ytv_alltime[fold]=np.hstack([ Ytv[0][fold], Ytv[15][fold],Ytv[30][fold], Ytv[45][fold], Ytv[60][fold],
Ytv[90][fold], Ytv[120][fold],Ytv[180][fold], Ytv[240][fold] ])
Ytest_alltime[fold]=np.hstack([ Ytest[0][fold], Ytest[15][fold], Ytest[30][fold], Ytest[45][fold], Ytest[60][fold],
Ytest[90][fold], Ytest[120][fold], Ytest[180][fold], Ytest[240][fold] ])
# %%
Xtv_alltime[fold].shape
# %%
xgb_est=100
xgb_lr=0.3
xgb_depth=6
xgb_subsamp=1
xgb_colsam=1
rf_estimators=100
rf_depth=None
meric_Threshold=0.00
# %%
from xgboost import XGBClassifier
features_lst=[]
test_acc_lst=[]
best_features_lst=[]
for f in range(1,11):
features=feature_selection(Xtv_alltime[f], Ytv_alltime[f],
xgb_est,xgb_lr,xgb_depth,xgb_subsamp,xgb_colsam,
rf_estimators, rf_depth, meric_Threshold,numb_runs=100)
features_lst.append(features)
features=features[0:(get_best_n_features(features, Xtv_alltime[f], Ytv_alltime[f], rf_estimators, rf_depth)[0])]
best_features_lst.append(features)
print('best_features_lst',best_features_lst)
best_features_acc_lst=(get_best_n_features(features, Xtv_alltime[f], Ytv_alltime[f], rf_estimators, rf_depth)[2])
print('best_features_lst_ACC',best_features_acc_lst)
print('######################## ACC ######################')
# %%
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import (recall_score, confusion_matrix, balanced_accuracy_score,
f1_score, precision_score, roc_auc_score,
matthews_corrcoef, average_precision_score)
original_features =[5991, 3529, 30022, 38491, 38539, 36613]
# Initialize dictionaries to store metrics for each fold
recall_scores = {}
normal_accuracies={}
balanced_accuracies = {}
f1_scores = {}
precision_scores = {}
roc_auc_scores = {}
pr_auc_scores = {}
confusion_matrices={}
for fold in range(1, 11):
# Train the model on the specified features in the training set
rf = RandomForestClassifier(random_state=42, class_weight='balanced')
rf.fit(np.vstack([ Xtv_alltime[fold][:, original_features] ]),
np.hstack([ Ytv_alltime[fold]]) )
# Predict on the test set
y_pred = rf.predict(Xtest_alltime[fold][:, original_features])
y_prob = rf.predict_proba(Xtest_alltime[fold][:, original_features])[:, 1] # Probability for ROC-AUC and PR-AUC
recall = recall_score(Ytest_alltime[fold], y_pred, pos_label=1)
recall_scores[fold] = recall
cm = confusion_matrix(Ytest_alltime[fold], y_pred)
tn, fp, fn, tp = cm.ravel()
normal_acc = accuracy_score(Ytest_alltime[fold], y_pred)
normal_accuracies[fold] = normal_acc
balanced_acc = balanced_accuracy_score(Ytest_alltime[fold], y_pred)
balanced_accuracies[fold] = balanced_acc
f1 = f1_score(Ytest_alltime[fold], y_pred, pos_label=1)
f1_scores[fold] = f1
precision = precision_score(Ytest_alltime[fold], y_pred, pos_label=1)
precision_scores[fold] = precision
roc_auc = roc_auc_score(Ytest_alltime[fold], y_prob)
roc_auc_scores[fold] = roc_auc
pr_auc = average_precision_score(Ytest_alltime[fold], y_prob)
pr_auc_scores[fold] = pr_auc
cm = confusion_matrix(Ytest_alltime[fold], y_pred)
confusion_matrices[fold] = cm
# Display metrics for each fold
for fold in range(1, 11):
print(f"Fold {fold}:")
print(f" Recall: {recall_scores[fold]:.4f}")
print(f" Balanced Accuracy: {balanced_accuracies[fold]:.4f}")
print(f" Normal Accuracy: {normal_accuracies[fold]:.4f}")
print(f" F1 Score: {f1_scores[fold]:.4f}")
print(f" Precision: {precision_scores[fold]:.4f}")
print(f" ROC-AUC: {roc_auc_scores[fold]:.4f}")
print(f" Precision-Recall AUC: {pr_auc_scores[fold]:.4f}")
print(f" CM: {confusion_matrices[fold]}")
# Calculate averages and standard deviations for each metric
avg_recall = np.mean(list(recall_scores.values()))
std_recall = np.std(list(recall_scores.values()))
avg_normal_accuracy = np.mean(list(normal_accuracies.values()))
std_normal_accuracy = np.std(list(normal_accuracies.values()))
avg_balanced_accuracy = np.mean(list(balanced_accuracies.values()))
std_balanced_accuracy = np.std(list(balanced_accuracies.values()))
avg_f1 = np.mean(list(f1_scores.values()))
std_f1 = np.std(list(f1_scores.values()))
avg_precision = np.mean(list(precision_scores.values()))
std_precision = np.std(list(precision_scores.values()))
avg_roc_auc = np.mean(list(roc_auc_scores.values()))
std_roc_auc = np.std(list(roc_auc_scores.values()))
avg_pr_auc = np.mean(list(pr_auc_scores.values()))
std_pr_auc = np.std(list(pr_auc_scores.values()))
# Display the averages and standard deviations
print("\nAverage Metrics Across All Folds (with Standard Deviations):")
print(f" Average Recall: {avg_recall:.4f} ± {std_recall:.4f}")
print(f" Average Normal Accuracy: {avg_normal_accuracy:.4f} ± {std_normal_accuracy:.4f}")
print(f" Average Balanced Accuracy: {avg_balanced_accuracy:.4f} ± {std_balanced_accuracy:.4f}")
print(f" Average F1 Score: {avg_f1:.4f} ± {std_f1:.4f}")
print(f" Average Precision: {avg_precision:.4f} ± {std_precision:.4f}")
print(f" Average ROC-AUC: {avg_roc_auc:.4f} ± {std_roc_auc:.4f}")
print(f" Average Precision-Recall AUC: {avg_pr_auc:.4f} ± {std_pr_auc:.4f}")
# %% [markdown]
# # ------------------------------------------------------------
# %%
from sklearn.model_selection import StratifiedKFold
def cross_val_with_features(X, y):
skf = StratifiedKFold(n_splits=10, shuffle=True, random_state=42)
accuracies = []
for train_index, test_index in skf.split(X, y):
X_train, X_test = X[train_index], X[test_index]
y_train, y_test = y[train_index], y[test_index]
rf.fit(X_train, y_train)
y_pred =rf.predict(X_test)
accuracies.append(accuracy_score(y_test, y_pred))
return accuracies
X_3mon_ob=np.vstack([ Xdata_ob[0][30:60,:], Xdata_ob[15][30:60,:], Xdata_ob[30][30:60,:], Xdata_ob[45][30:60,:],
Xdata_ob[60][30:60,:],Xdata_ob[90][30:60,:], Xdata_ob[120][30:60,:],
Xdata_ob[180][30:60,:], Xdata_ob[240][30:60,:] ])
X_bl_ob=np.vstack([ Xdata_ob[0][0:30,:], Xdata_ob[15][0:30,:], Xdata_ob[30][0:30,:], Xdata_ob[45][0:30,:],
Xdata_ob[60][0:30,:],Xdata_ob[90][0:30,:], Xdata_ob[120][0:30,:],
Xdata_ob[180][0:30,:], Xdata_ob[240][0:30,:] ])
X_pwl_ob=np.vstack([ Xdata_ob[0][60:,:], Xdata_ob[15][60:,:], Xdata_ob[30][60:,:], Xdata_ob[45][60:,:],
Xdata_ob[60][60:,:],Xdata_ob[90][60:,:], Xdata_ob[120][60:,:],
Xdata_ob[180][60:,:], Xdata_ob[240][60:,:] ])
X_bl_lean=np.vstack([ Xdata_Lean[0][0:30,:], Xdata_Lean[15][0:30,:], Xdata_Lean[30][0:30,:], Xdata_Lean[45][0:30,:],
Xdata_Lean[60][0:30,:],Xdata_Lean[90][0:30,:], Xdata_Lean[120][0:30,:],
Xdata_Lean[180][0:30,:], Xdata_Lean[240][0:30,:] ])
# %% [markdown]
# ### Test minimum best model per Stage per Timepoint
# %%
best_feat=[ 7144, 6001, 38004, 6674, 37885, 7040, 38470, 37416, 7150,
38520, 7016, 37372, 1320, 7115, 7028, 7119, 36576]
best_feat_permu_min=[5991, 3529, 30022, 38491, 38539, 36613]
# %%
def acc_per_stage_per_tps(best_feat,time_point):
subj_y = np.hstack([np.ones(30), np.zeros(30)])
Xtv_id = {}
Xtest_id = {}
index_lean = {}
index_ob = {}
for s in range(10):
index_lean[s] = ([np.where(np.array(subj_y) == 1)][0][0])[s * (60 // 10 // 2):(s + 1) * (60 // 10 // 2)]
index_ob[s] = ([np.where(np.array(subj_y) == 0)][0][0])[s * (60 // 10 // 2):(s + 1) * (60 // 10 // 2)]
index_test = index_lean[s].tolist() + index_ob[s].tolist()
Xtest_id[s] = index_test
index_train = np.arange(60)[np.isin(np.arange(60), index_test, invert=True)]
Xtv_id[s] = index_train
allstage_score = []
for stage, obese_data in zip(['BL', '3MON', 'PWL'],
[Xdata_ob[time_point][:30, :], Xdata_ob[time_point][30:60, :], Xdata_ob[time_point][60:, :]]):
Xdata = np.vstack([obese_data[:, best_feat], Xdata_Lean[time_point][:30, :][:, best_feat]])
Ydata = np.hstack([np.ones(30), np.zeros(30)])
avg_score = []
for s in range(10):
Xtv = Xdata[Xtv_id[s]]
Ytv = Ydata[Xtv_id[s]]
Xtest = Xdata[Xtest_id[s]]
Ytest = Ydata[Xtest_id[s]]
rf.fit(Xtv, Ytv)
avg_score.append(accuracy_score(rf.predict(Xtest), Ytest))
allstage_score.append(avg_score)
return sum(allstage_score[0]) / 10, sum(allstage_score[1]) / 10, sum(allstage_score[2]) / 10
# %%
time_points = [0, 15, 30, 45, 60, 90, 120, 180, 240]
pwl_acc_list = []
mon3_acc_list = []
bl_acc_list = []
# Calculate the accuracy for each time point
for time_point in time_points:
bl_acc, mon3_acc, pwl_acc = acc_per_stage_per_tps(best_feat_permu_min,time_point)
pwl_acc_list.append(pwl_acc)
mon3_acc_list.append(mon3_acc)
bl_acc_list.append(bl_acc)
# %%
from scipy.stats import ttest_ind, levene, mannwhitneyu, shapiro
def compare_lists(list1, list2, name1, name2):
# Perform Shapiro-Wilk test for normality on both lists
shapiro_list1 = shapiro(list1)
shapiro_list2 = shapiro(list2)
print(f"Shapiro-Wilk test p-values: {name1} = {shapiro_list1.pvalue:.4f}, {name2} = {shapiro_list2.pvalue:.4f}")
# Check if both distributions are normally distributed (p-value > 0.05)
if shapiro_list1.pvalue > 0.05 and shapiro_list2.pvalue > 0.05:
# Check for equality of variances using Levene's Test
stat, p_value_levene = levene(list1, list2)
equal_var = p_value_levene > 0.05
# Perform an independent t-test
t_stat, p_value = ttest_ind(list1, list2, equal_var=equal_var)
test_name = "Independent t-test"
else:
# Use Mann-Whitney U test if either distribution is not normal
t_stat, p_value = mannwhitneyu(list1, list2)
test_name = "Mann-Whitney U test"
print(f"{test_name} ({name1} vs {name2}): statistic={t_stat:.4f}, p-value={p_value:.4f}\n")
# Compare the lists in pairs
compare_lists(pwl_acc_list, mon3_acc_list, "pwl_acc_list", "mon3_acc_list")
compare_lists(pwl_acc_list, bl_acc_list, "pwl_acc_list", "bl_acc_list")
compare_lists(mon3_acc_list, bl_acc_list, "mon3_acc_list", "bl_acc_list")
# %%
time_points_str = [str(tp) for tp in time_points]
plt.figure(figsize=(8, 4))
plt.plot(time_points_str, bl_acc_list, label='BL', marker='o')
plt.fill_between(time_points_str, bl_acc_list - np.std(bl_acc_list), bl_acc_list + np.std(bl_acc_list), alpha=0.15,color='blue')
plt.plot(time_points_str, pwl_acc_list, label='PWL', marker='o')
plt.fill_between(time_points_str, pwl_acc_list - np.std(pwl_acc_list), pwl_acc_list + np.std(pwl_acc_list), alpha=0.15,color='orange')
plt.plot(time_points_str, mon3_acc_list, label='3MON', marker='o')
plt.fill_between(time_points_str, mon3_acc_list - np.std(mon3_acc_list), mon3_acc_list + np.std(mon3_acc_list), alpha=0.15, color='green')
plt.xlabel('Time Point')
plt.ylabel('Accuracy')
plt.title('Accuracy per Stage per Time Point')
plt.legend()
plt.grid(True)
plt.ylim(0.2, 1)
# %% [markdown]
# #### Model Significance Test
# %%
def swap_labels(X, y, n, swap_ratio):
np.random.seed(42) # For reproducibility
X_swapped, y_swapped = X.copy(), y.copy()
for _ in range(n):
pain_indices = np.where(y_swapped == 1)[0]
healthy_indices = np.where(y_swapped == 0)[0]
num_to_swap = int(min(len(pain_indices), len(healthy_indices)) * swap_ratio)
pain_to_healthy_indices = np.random.choice(pain_indices, size=num_to_swap, replace=False)
healthy_to_pain_indices = np.random.choice(healthy_indices, size=num_to_swap, replace=False)
y_swapped[pain_to_healthy_indices], y_swapped[healthy_to_pain_indices] = \
y[healthy_to_pain_indices], y[pain_to_healthy_indices]
return X_swapped, y_swapped
# %%
num_splits=1000
swap_ratio=0.5
def permutation_test(real_accuracies, shuffled_accuracies, n_permutations=1000):
observed_diff = np.mean(real_accuracies) - np.mean(shuffled_accuracies)
all_accuracies = np.concatenate([real_accuracies, shuffled_accuracies])
extreme_count = 0
for i in range(n_permutations):
# Shuffle labels
np.random.shuffle(all_accuracies)
perm_real = all_accuracies[:len(real_accuracies)]
perm_shuffled = all_accuracies[len(real_accuracies):]
perm_diff = np.mean(perm_real) - np.mean(perm_shuffled)
# Count how often permuted difference is as extreme as observed
if abs(perm_diff) >= abs(observed_diff):
extreme_count += 1
p_value = (extreme_count + 1) / (n_permutations + 1)
return p_value, observed_diff
clf=rf
for t in range(0, 9):
shuffled_acc = []
non_shuffled_acc = []
# Adjust your xx and yy variable population based on t
xx = np.vstack([X_bl_ob[t*30:(t+1)*30,best_feat_permu_min], X_bl_lean[t*30:(t+1)*30,best_feat_permu_min]])
yy = np.hstack([np.zeros(30), np.ones(30)])
# Calculate shuffled accuracy
for i in range(num_splits):
X_train, X_test, y_train, y_test = train_test_split(xx, yy, test_size=0.2, random_state=i, stratify=yy)
X_swapped, y_swapped = swap_labels(X_train, y_train, 1, swap_ratio)
clf.fit(X_train, y_swapped)
score = accuracy_score(y_test, clf.predict(X_test))
shuffled_acc.append(score)
print('shuf acc', sum(shuffled_acc) / len(shuffled_acc), 'shuf std', statistics.stdev(shuffled_acc))
# Calculate non-shuffled accuracy
for i in range(num_splits):
X_train, X_test, y_train, y_test = train_test_split(xx, yy, test_size=0.2, random_state=i, stratify=yy)
clf.fit(X_train, y_train)
score = accuracy_score(y_test, clf.predict(X_test))
non_shuffled_acc.append(score)
print('real acc', sum(non_shuffled_acc) / len(non_shuffled_acc), 'real std', statistics.stdev(non_shuffled_acc))
p_value, diff = permutation_test(non_shuffled_acc, shuffled_acc, n_permutations=1000)
print(f"Time {t}: p-value = {p_value:.4f}, mean difference = {diff:.3f}")
# %%
#os.chdir('/home/jupy/SourceAnalysis/Journal2024_updatedMICMAC/SourceData/PWL_shuffle')
num_splits=1000
swap_ratio=0.5
def permutation_test(real_accuracies, shuffled_accuracies, n_permutations=1000):
observed_diff = np.mean(real_accuracies) - np.mean(shuffled_accuracies)
all_accuracies = np.concatenate([real_accuracies, shuffled_accuracies])
extreme_count = 0
for i in range(n_permutations):
# Shuffle labels
np.random.shuffle(all_accuracies)
perm_real = all_accuracies[:len(real_accuracies)]
perm_shuffled = all_accuracies[len(real_accuracies):]
perm_diff = np.mean(perm_real) - np.mean(perm_shuffled)
# Count how often permuted difference is as extreme as observed
if abs(perm_diff) >= abs(observed_diff):
extreme_count += 1
p_value = (extreme_count + 1) / (n_permutations + 1)
return p_value, observed_diff
clf=rf
for t in range(0, 9):
shuffled_acc = []
non_shuffled_acc = []
# Adjust your xx and yy variable population based on t
xx = np.vstack([X_pwl_ob[t*30:(t+1)*30,best_feat_permu_min], X_bl_lean[t*30:(t+1)*30,best_feat_permu_min]])
yy = np.hstack([np.zeros(30), np.ones(30)])
# Calculate shuffled accuracy
for i in range(num_splits):
X_train, X_test, y_train, y_test = train_test_split(xx, yy, test_size=0.2, random_state=i, stratify=yy)
X_swapped, y_swapped = swap_labels(X_train, y_train, 1, swap_ratio)
clf.fit(X_train, y_swapped)
score = accuracy_score(y_test, clf.predict(X_test))
shuffled_acc.append(score)
print('shuf acc', sum(shuffled_acc) / len(shuffled_acc), 'shuf std', statistics.stdev(shuffled_acc))
# Calculate non-shuffled accuracy
for i in range(num_splits):
X_train, X_test, y_train, y_test = train_test_split(xx, yy, test_size=0.2, random_state=i, stratify=yy)
clf.fit(X_train, y_train)
score = accuracy_score(y_test, clf.predict(X_test))
non_shuffled_acc.append(score)
print('real acc', sum(non_shuffled_acc) / len(non_shuffled_acc), 'real std', statistics.stdev(non_shuffled_acc))
p_value, diff = permutation_test(non_shuffled_acc, shuffled_acc, n_permutations=1000)
print(f"Time {t}: p-value = {p_value:.4f}, mean difference = {diff:.3f}")
# %%
num_splits=1000
swap_ratio=0.5
def permutation_test(real_accuracies, shuffled_accuracies, n_permutations=1000):
observed_diff = np.mean(real_accuracies) - np.mean(shuffled_accuracies)
all_accuracies = np.concatenate([real_accuracies, shuffled_accuracies])
extreme_count = 0
for i in range(n_permutations):
# Shuffle labels
np.random.shuffle(all_accuracies)
perm_real = all_accuracies[:len(real_accuracies)]
perm_shuffled = all_accuracies[len(real_accuracies):]
perm_diff = np.mean(perm_real) - np.mean(perm_shuffled)
# Count how often permuted difference is as extreme as observed
if abs(perm_diff) >= abs(observed_diff):
extreme_count += 1
p_value = (extreme_count + 1) / (n_permutations + 1)
return p_value, observed_diff
clf=rf
for t in range(0, 9):
shuffled_acc = []
non_shuffled_acc = []
# Adjust your xx and yy variable population based on t
xx = np.vstack([X_3mon_ob[t*30:(t+1)*30,best_feat_permu_min], X_bl_lean[t*30:(t+1)*30,best_feat_permu_min]])
yy = np.hstack([np.zeros(30), np.ones(30)])
# Calculate shuffled accuracy
for i in range(num_splits):
X_train, X_test, y_train, y_test = train_test_split(xx, yy, test_size=0.2, random_state=i, stratify=yy)
X_swapped, y_swapped = swap_labels(X_train, y_train, 1, swap_ratio)
clf.fit(X_train, y_swapped)
score = accuracy_score(y_test, clf.predict(X_test))
shuffled_acc.append(score)
print('shuf acc', sum(shuffled_acc) / len(shuffled_acc), 'shuf std', statistics.stdev(shuffled_acc))
# Calculate non-shuffled accuracy
for i in range(num_splits):
X_train, X_test, y_train, y_test = train_test_split(xx, yy, test_size=0.2, random_state=i, stratify=yy)
clf.fit(X_train, y_train)
score = accuracy_score(y_test, clf.predict(X_test))
non_shuffled_acc.append(score)
print('real acc', sum(non_shuffled_acc) / len(non_shuffled_acc), 'real std', statistics.stdev(non_shuffled_acc))
p_value, diff = permutation_test(non_shuffled_acc, shuffled_acc, n_permutations=1000)
print(f"Time {t}: p-value = {p_value:.4f}, mean difference = {diff:.3f}")
# %% [markdown]
# # Feature Interpretation
# %%
Original_features= [5991, 3529, 30022, 38491, 38539, 36613]
# %%
def extract_feat_info(feat_idx_ls):
os.chdir('/home/jupy/SourceAnalysis')
coordinates=np.loadtxt('88_areas-ROIcentroids.txt')[:,0:3]
sensor= np.arange(0,88)#.astype(str)
nodes1=np.zeros([2])
nodes2=np.zeros([2])
data=np.array(pd.read_csv('88ROIallBA-ROI-ROI.csv',header=None))
bands=['delta','theta','alpha','beta','gamma']
featid_info=[]
structure_info=[]
band_info=[]
for feat_idx in feat_idx_ls:
if feat_idx//7744==0:
e_color='BuGn'
if feat_idx//7744==1:
e_color='inferno'
if feat_idx//7744==2:
e_color='Oranges'
if feat_idx//7744==3:
e_color='Purples'
if feat_idx//7744==4:
e_color='Oranges'
#print(feat_idx)
i=feat_idx-feat_idx//7744*7744
chanel1=int((i)//88)
name_chan1=sensor[chanel1]
name_chan1_correctInd=name_chan1+1
chanel2=int((i)%88)
name_chan2=sensor[chanel2]
name_chan2_correctInd=name_chan2+1
#lobe
lobe1=data[np.where(data[:,-1]==name_chan1_correctInd)[0][0],3]
lobe2=data[np.where(data[:,-1]==name_chan2_correctInd)[0][0],3]
#structure
structure1=data[np.where(data[:,-1]==name_chan1_correctInd)[0][0],4]
structure2=data[np.where(data[:,-1]==name_chan2_correctInd)[0][0],4]
#Brodmann area
brod1=data[np.where(data[:,-1]==name_chan1_correctInd)[0][0],5]
brod2=data[np.where(data[:,-1]==name_chan2_correctInd)[0][0],5]
#ROI number
roinumb1=data[np.where(data[:,-1]==name_chan1_correctInd)[0][0],6]
roinumb2=data[np.where(data[:,-1]==name_chan2_correctInd)[0][0],6]
nodes1=np.vstack([nodes1,np.array([name_chan1,name_chan2])])
nodes2=np.vstack([nodes2,np.array([name_chan2,name_chan1])])
print('feature_id: ',feat_idx,';',bands[feat_idx//7744],';',lobe1,'-',lobe2,';',structure1,'-',structure2,';',
brod1,'-',brod2,';',roinumb1,'-',roinumb2)
#print(' ')
structure_info.append(f'{brod1} - {brod2}')
featid_info.append(feat_idx)
band_info.append(bands[feat_idx//7744])
return featid_info,structure_info,band_info
# %%
Original_features
featid,featStruc,feat_band=extract_feat_info(Original_features)
featStruc= ['40L - 4R','25L - 5R', '44R - 8L','32L - 22R','32L - 47R','38L - 3R',]
# %%
#### Feature Significance Test
# %%
df_simulated = pd.DataFrame({ 'Feature_ID':col_featid,'Feature_Structure':col_struc,'Feature_Value': col_featvalue,'Group': col_group})
# %%
from scipy.stats import mannwhitneyu
reject_array=[]
pvalue_array=[]
sigstar_array=[]
for tp in [0,15,30,45, 60,90,120,180,240]:
featStruc
col_struc=[]
for i in featStruc:
col_struc.append([f'{i}']*90)
col_struc=col_struc+col_struc
col_struc=np.array(col_struc).reshape(-1)
col_featvalue_ob=[]
col_featvalue_lean=[]
for feat in Original_features:
aa=Xdata_ob[tp][:, feat]
bb= Xdata_Lean[tp][:, feat]
col_featvalue_ob.append(aa)
col_featvalue_lean.append(bb)
col_featvalue_ob=np.array(col_featvalue_ob).reshape(-1)
col_featvalue_lean=np.array(col_featvalue_lean).reshape(-1)
ros = RandomOverSampler(sampling_strategy=0.33)
col_featvalue_lean_rosX, col_featvalue_lean_rosY = ros.fit_resample(col_featvalue_lean, y)
print(col_featvalue_ob.shape)
print(col_featvalue_lean.shape)
col_featvalue=np.hstack([col_featvalue_ob,col_featvalue_lean])
col_featid=[]
for feat in Original_features:
col_featid.append( np.array([feat]*90).astype(str) )
col_featid=col_featid+col_featid
col_featid=np.array(col_featid).reshape(-1)
col_featid
print(col_featid.shape)
col_group=['Obese']*int(col_featvalue.shape[0]/2)+['Lean']*int(col_featvalue.shape[0]/2)
col_group=np.array(col_group)
print(col_group.shape)
df_simulated = pd.DataFrame({ 'Feature_ID':col_featid,'Feature_Structure':col_struc,'Feature_Value': col_featvalue,'Group': col_group})
plt.rcParams.update({'font.size': 9})
import matplotlib.pyplot as plt
import seaborn as sns
from statannot import add_stat_annotation
from scipy.stats import ttest_ind
import os
from statsmodels.stats.multitest import multipletests
# Define box pairs for annotation
unique_features = df_simulated["Feature_Structure"].unique()
box_pairs = [((feature, "Obese"), (feature, "Lean")) for feature in unique_features]
box_pairs = [((feature, "Lean"), (feature, "Obese")) for feature in unique_features]
# Perform t-tests
p_values = []
for pair in box_pairs:
group1 = df_simulated[(df_simulated['Feature_Structure'] == pair[0][0]) & (df_simulated['Group'] == pair[0][1])]['Feature_Value']
group2 = df_simulated[(df_simulated['Feature_Structure'] == pair[1][0]) & (df_simulated['Group'] == pair[1][1])]['Feature_Value']
#t_stat, p_val = ttest_ind(group1, group2)
t_stat, p_val = mannwhitneyu(group1, group2, alternative='two-sided')
p_values.append(p_val)
reject, pvals_corrected, _, _ = multipletests(p_values, alpha=0.05, method='fdr_bh')
# Output the results
print("Reject null hypothesis for these comparisons:", reject)
print("Corrected p-values:", [f'{pval:.4f}' for pval in pvals_corrected])
def get_significance_level(p):
if p < 0.0001:
return "***"
elif p < 0.01:
return "**"
elif p < 0.05:
return "*"
else:
return "ns"
# Get significance level annotations for each feature
significance_levels = [get_significance_level(p) for p in pvals_corrected]
significance_levels
reject_array.append(reject)
pvalue_array.append(pvals_corrected)
sigstar_array.append(significance_levels)
pvalue_array=np.array(pvalue_array)
# %%
### Feature Distribution
# %%
# Define the feature and time indices
time_index = ["T0", "T15", "T30", "T45", "T60", "T90", "T120", "T180", "T240"]
feature_index = ['40L - 4R', '25L - 5R', '44R - 8L', '32L - 22R', '32L - 47R', '38L - 3R']
# Set custom colors for specific feature labels
custom_colors = {
'40L - 4R': '#D4A000',
'25L - 5R': '#D4A000',
'44R - 8L': 'green',
'32L - 22R': 'red',
'32L - 47R': 'red',
'38L - 3R': 'red'
}
# Create the heatmap with features on the x-axis and timepoints on the y-axis
plt.rcParams.update({'font.size': 12})
plt.figure(figsize=(12, 8))
heatmap = plt.imshow(pvalue_array, cmap='YlOrBr_r', aspect='auto', vmin=0, vmax=0.05) # Increase color resolution for low p-values
# Add a color bar to indicate the scale
cbar = plt.colorbar(heatmap)
cbar.set_label('Corrected p-values', fontsize=12) # Increase the color bar label font size
cbar.ax.tick_params(labelsize=11) # Increase the font size of the color bar ticks
# Set the tick labels with custom colors
xticks = plt.xticks(np.arange(len(feature_index)), feature_index, rotation=45, ha='right', fontsize=12)
yticks = plt.yticks(np.arange(len(time_index)), time_index, fontsize=12)
for label in plt.gca().get_xticklabels():
label.set_color(custom_colors.get(label.get_text(), 'black'))
# Annotate the heatmap with p-values or "ns" if p > 0.05, with white font and increased font size
for i in range(pvalue_array.shape[0]):
for j in range(pvalue_array.shape[1]):
value = pvalue_array[i, j]
if value > 0.05:
plt.text(j, i, 'ns', ha='center', va='center', color='black', fontsize=12)
else:
plt.text(j, i, f'{value:.4f}', ha='center', va='center', color='white', fontsize=12)
# %%
featStruc
col_struc=[]
for i in featStruc:
col_struc.append([f'{i}']*270)
col_struc=col_struc+col_struc
col_struc=np.array(col_struc).reshape(-1)
col_featvalue_ob=[]
col_featvalue_lean=[]
for feat in Original_features:
aa=X_bl_ob[:, feat]
bb=X_bl_lean[:, feat]
col_featvalue_ob.append(aa)
col_featvalue_lean.append(bb)
col_featvalue_ob=np.array(col_featvalue_ob).reshape(-1)
col_featvalue_lean=np.array(col_featvalue_lean).reshape(-1)
print(col_featvalue_ob.shape)
print(col_featvalue_lean.shape)
col_featvalue=np.hstack([col_featvalue_ob,col_featvalue_lean])
col_featid=[]
for feat in Original_features:
# col_featid.append( np.array([feat]*90).astype(str) )
col_featid.append( np.array([feat]*270).astype(str) )
col_featid=col_featid+col_featid
col_featid=np.array(col_featid).reshape(-1)
col_featid
print(col_featid.shape)
col_group=['Obese']*int(col_featvalue.shape[0]/2)+['Lean']*int(col_featvalue.shape[0]/2)
col_group=np.array(col_group)
print(col_group.shape)
df_simulated = pd.DataFrame({ 'Feature_ID':col_featid,'Feature_Structure':col_struc,'Feature_Value': col_featvalue,'Group': col_group})
plt.rcParams.update({'font.size': 9})
import matplotlib.pyplot as plt
import seaborn as sns
from statannot import add_stat_annotation
from scipy.stats import ttest_ind
import os
# Set up the plot
plt.figure(figsize=(10,7))
custom_colors = {
'40L - 4R': '#D4A000',
'25L - 5R': '#D4A000',
'44R - 8L': 'green',
'32L - 22R': 'red',
'32L - 47R': 'red',
'38L - 3R': 'red'
}
ax = sns.boxplot(x='Feature_Structure', y='Feature_Value', hue='Group', data=df_simulated)
unique_features = df_simulated["Feature_Structure"].unique()
box_pairs = [((feature, "Obese"), (feature, "Lean")) for feature in unique_features]
box_pairs = [((feature, "Lean"), (feature, "Obese")) for feature in unique_features]
p_values = []
for pair in box_pairs:
group1 = df_simulated[(df_simulated['Feature_Structure'] == pair[0][0]) & (df_simulated['Group'] == pair[0][1])]['Feature_Value']
group2 = df_simulated[(df_simulated['Feature_Structure'] == pair[1][0]) & (df_simulated['Group'] == pair[1][1])]['Feature_Value']
t_stat, p_val = mannwhitneyu(group1, group2, alternative='two-sided')
p_values.append(p_val)
reject, pvals_corrected, _, _ = multipletests(p_values, alpha=0.05, method='fdr_bh')
significance_annotations = ['***' if p < 0.0001 else
'**' if p < 0.001 else
'*' if p < 0.01 else
'ns' for p in pvals_corrected]
add_stat_annotation(ax, data=df_simulated, x='Feature_Structure', y='Feature_Value', hue='Group',
box_pairs=box_pairs, perform_stat_test=False, pvalues=p_values, test=None, text_format='star', loc='outside',
line_height=0.001, text_annot_custom=significance_annotations, fontsize=12)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.set_xlabel(None)
plt.ylabel('Feature Value', fontsize=12)
ax.legend(loc='upper left', bbox_to_anchor=(1, 1), fontsize=12)
plt.xticks(rotation=45, ha='right', fontsize=12)
for label in plt.gca().get_xticklabels():
label.set_color(custom_colors.get(label.get_text(), 'black'))
plt.yticks(fontsize=12)
plt.tight_layout()

CodeScript.ipynb, under CC-BY-4.0 · at the source

Overview

Authors: Yuan Yue¹, Patrick Manning², Dirk De Ridder³, Matthew Hall³, Divya Bharatkumar Adhia³, Samantha Ross², Daniel Alencar da Costa¹, Jeremiah D. Deng¹

ORCID iDs: Yuan Yue, Matthew Hall, Divya Bharatkumar Adhia, Jeremiah D. Deng

Institutions: University of Otago (New Zealand)

Journal: Communications medicine, volume 6, issue 1, article 241

Dates: received 7 October 2024; accepted 27 February 2026; published online 10 March 2026

Type: Research article · Language: English

License: CC BY-NC-ND

Identifiers: DOI 10.1038/s43856-026-01518-5 · PMID 41807824 · PMCID PMC13106778 · OpenAlex W7134936701

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: human (organism), other condition (population), clinical / translational (subfield)

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

Keywords: Neurology, Predictive markers, Endocrinology

Topic: Transcranial Magnetic Stimulation Studies (Neurology, Neuroscience), according to OpenAlex

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

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

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

yuan410/obesitymeta

License: Apache-2.0

State: the link answers, verified on 30 September 2026

Evidence: files inventoried

Commit: a2b735338c6bbccc8bf4a0cfdd28796904f36c02, 10 October 2025

Languages: Jupyter (1)

Size: 1,445 files, 1 script

Software Heritage: not archived

Found in: “Code availability”

Holds: README, license file, environment (Dockerfile, requirements.txt), 1 notebook

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

Tools: imbalanced-learn (1 file), Matplotlib (1 file), NumPy (1 file), pandas (1 file), scikit-learn (1 file), SciPy (1 file), seaborn (1 file), statannotations (1 file), statsmodels (1 file), XGBoost (1 file)

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

30 September 2026: the link answers

3 files

CodeScript.ipynb, Jupyter, 1,336 lines
LICENSE.md, License, 201 lines
README.md, Text, 65 lines

At the source: github.com/yuan410/ObesityMeta/tree/main

figshare 30351418

License: CC-BY-4.0

State: the link answers, verified on 30 September 2026

Evidence: files inventoried

Languages: Jupyter (1)

Size: 1 file, 1 script

Software Heritage: not checked

Found in: “Code availability”

Holds: 1 notebook

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Availability: 1 check, the latest on 30 September 2026: the link answers (HTTP 200)

30 September 2026: the link answers (HTTP 200)

1 file

CodeScript.ipynb, Jupyter, 1,336 lines, 5 matches

At the source:

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to the authors' code: figshare 30351418, yuan410/obesitymeta

Read it in the paper: doi.org/10.1038/s43856-026-01518-5.

Tracing map

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

What the map holds:

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

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

Data

Datasets cited

figshare:30351280, at figshare; found in “Data availability”

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

it points to a dataset: figshare 30351280
it says that the data are available on request

Read it in the paper: doi.org/10.1038/s43856-026-01518-5.

Versions

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

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 8 authors, 3 keywords, 86 references.

Cite

This paper

Yue, Y., Manning, P., De Ridder, D., Hall, M., Adhia, D. B., Ross, S., Alencar da Costa, D., & Deng, J. D. (2026). Machine learning-based identification of abnormal functional connectivity in obesity across different metabolic states. Communications medicine, 6(1), 241. https://doi.org/10.1038/s43856-026-01518-5

BibTeX

@article{yue2026machine, author = {Yue, Yuan and Manning, Patrick and De Ridder, Dirk and Hall, Matthew and Adhia, Divya Bharatkumar and Ross, Samantha and Alencar da Costa, Daniel and Deng, Jeremiah D.}, title = {{Machine learning-based identification of abnormal functional connectivity in obesity across different metabolic states}}, journal = {Communications medicine}, year = {2026}, month = mar, volume = {6}, number = {1}, pages = {241}, publisher = {Nature Publishing Group}, issn = {2730-664X}, doi = {10.1038/s43856-026-01518-5}, url = {https://doi.org/10.1038/s43856-026-01518-5}, pmid = {41807824}, pmcid = {PMC13106778} }

RIS

TY - JOUR AU - Yue, Yuan AU - Manning, Patrick AU - De Ridder, Dirk AU - Hall, Matthew AU - Adhia, Divya Bharatkumar AU - Ross, Samantha AU - Alencar da Costa, Daniel AU - Deng, Jeremiah D. TI - Machine learning-based identification of abnormal functional connectivity in obesity across different metabolic states T2 - Communications medicine J2 - Commun Med (Lond) PY - 2026 DA - 2026/03/10 VL - 6 IS - 1 SP - 241 SN - 2730-664X PB - Nature Publishing Group DO - 10.1038/s43856-026-01518-5 UR - https://doi.org/10.1038/s43856-026-01518-5 LA - en ER -

CSL-JSON

{ "id": "10.1038/s43856-026-01518-5", "type": "article-journal", "title": "Machine learning-based identification of abnormal functional connectivity in obesity across different metabolic states", "container-title": "Communications medicine", "author": [ { "family": "Yue", "given": "Yuan" }, { "family": "Manning", "given": "Patrick" }, { "family": "De Ridder", "given": "Dirk" }, { "family": "Hall", "given": "Matthew" }, { "family": "Adhia", "given": "Divya Bharatkumar" }, { "family": "Ross", "given": "Samantha" }, { "family": "Alencar da Costa", "given": "Daniel" }, { "family": "Deng", "given": "Jeremiah D." } ], "container-title-short": "Commun Med (Lond)", "volume": "6", "issue": "1", "page": "241", "DOI": "10.1038/s43856-026-01518-5", "PMID": "41807824", "PMCID": "PMC13106778", "ISSN": "2730-664X", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s43856-026-01518-5", "language": "en", "issued": { "date-parts": [ [ 2026, 3, 10 ] ] } }

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

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.3389/fnins.2026.1803154 [code]: Multimodal machine learning reveals neurobiological signatures of binge-type eating disorders.
Journal: Frontiers in neuroscience
In common: statannotations, XGBoost, statsmodels, 6 other tools, other condition, 2 references
[2] doi:10.1016/j.isci.2026.115329 [code]: Brain metastases converge on shared geometric architecture and transcriptomic landscape yet remain distinct from gliomas.
Journal: iScience
In common: imbalanced-learn, statannotations, XGBoost, 7 other tools, clinical / translational, other condition
[3] doi:10.3390/s26175327 [code]: Subject Identity Confounds qEEG Emotion Recognition on DEAP and DREAMER.
Journal: Sensors (Basel, Switzerland)
In common: imbalanced-learn, XGBoost, statsmodels, 6 other tools, 1 reference
[4] doi:10.1038/s41598-026-64405-y [code]: Evaluating clinical and neuroimaging predictors for cognitive-behavioral therapy outcome in obsessive-compulsive disorder.
Journal: Scientific reports
In common: imbalanced-learn, XGBoost, statsmodels, 6 other tools, clinical / translational, other condition
[5] doi:10.1016/j.nicl.2026.104012 [code]: Structural-functional multilayer brain network properties and outcome of combined repetitive transcranial magnetic stimulation and psychotherapy for obsessive-compulsive disorder.
Journal: NeuroImage. Clinical
In common: imbalanced-learn, XGBoost, statsmodels, 6 other tools, other condition
[6] doi:10.1038/s42003-026-10957-8 [code]: Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: imbalanced-learn, XGBoost, statsmodels, 6 other tools
[7] doi:10.1186/s13059-026-04125-8 [code]: MLMarker: a machine learning framework for tissue inference and biomarker discovery.
Journal: Genome biology
In common: imbalanced-learn, XGBoost, statsmodels, 6 other tools
[8] doi:10.1038/s41598-026-56688-y [code]: On the value of radiomics in addition to clinical measures in emotional conflict fMRI for predicting sertraline response in major depressive disorder.
Journal: Scientific reports
In common: imbalanced-learn, XGBoost, statsmodels, 6 other tools
[9] doi:10.1016/j.isci.2026.116825 [code]: Social hierarchy shapes behavioral and transcriptional responses to chronic stress and ketamine in male mice.
Journal: iScience
In common: imbalanced-learn, XGBoost, statsmodels, 6 other tools
[10] doi:10.7554/elife.108109 [code]: Multimodal MRI marker of cognition explains the association between cognition and mental health in the UK Biobank.
Journal: eLife
In common: seaborn, scikit-learn, pandas, 3 other tools, author Jeremiah D Deng

Contribute

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

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

Request its removal

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

Discussion, reproductions, activity

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

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

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