Scripts to create a dataset from Redcap outputs to use for a PLS-DA classification.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

306 lines
7.9 KiB

# -*- coding: utf-8 -*-
"""
Created on Wed May 19 09:23:49 2021
@author: Dijkhofmf
"""
# Import stuff
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
#from pca import pca
from sklearn import preprocessing
from sklearn.model_selection import train_test_split, KFold, LeaveOneOut
from sklearn.cross_decomposition import PLSRegression
from sklearn.model_selection import RandomizedSearchCV, GridSearchCV, cross_val_score
from sklearn.metrics import r2_score, accuracy_score, plot_confusion_matrix, confusion_matrix, roc_curve, auc, roc_auc_score, precision_recall_fscore_support, classification_report, mean_squared_error
#%% Import data and path
Path = 'I:\Mike Dijkhof\Connecare MGP\Data\FinalFiles'
# Set path
os.chdir(Path)
#%% Create DF
FinalDF = pd.DataFrame(pd.read_csv('FinalDataset.csv'))
X = FinalDF.drop(['Pt Type'], axis=1)
X = X.drop(['Study ID'], axis=1)
# #%%
# plt.figure(dpi=720, figsize=(20,20))
# heatmap = sns.heatmap(X.corr(), vmin=-1, vmax=1, annot=True)
# plt.title('Pearson correlation heatmap')
# plt.figure(dpi=720, figsize=(20,20))
# heatmap = sns.heatmap(X.corr(method='spearman'), vmin=-1, vmax=1, annot=True)
# plt.title('Spearman correlation heatmap')
#%%
cols = pd.DataFrame(X).columns
y = FinalDF['Pt Type']
y = y.replace('Complication', int(1))
y = y.replace('Healthy', int(0))
X = X.dropna()
y = y.loc[X.index]
y = y.astype(int)
#X = X.reset_index()
#y = y.reset_index()
#y = y.drop('index', axis=1)
#X = X.drop('Study ID', axis=1)
X = preprocessing.scale(X)
#%%
#PLS = PLSRegression(n_components=2, scale=True)
# # Or reduce the data towards 2 PCs
# model = pca(n_components=8)
# # Fit transform
# results = model.fit_transform(X, y, col_labels=cols)
# # Plot explained variance
# fig, ax = model.plot()
# # Scatter first 2 PCs
# fig, ax = model.scatter()
# # Make biplot with the number of features
# fig, ax = model.biplot(n_feat=8)
#%%
def base_pls(X_train, y_train, X_test, n_components):
# Define PLS model and fit it to the train set
PLS = PLSRegression(n_components=n_components)
PLS.fit(X_train, y_train)
# Prediction on train set
y_c = PLS.predict(X_train)
# Prediction on test set
y_p = PLS.predict(X_test)
return(y_c, y_p)
#%%
cval = LeaveOneOut()
n_components = range(1,19)
X = pd.DataFrame(X)
yu = pd.DataFrame(y)
rmse_c = np.zeros(len(n_components)).astype('float32')
rmse_cv = np.zeros(len(n_components)).astype('float32')
R2_cv = np.zeros(len(n_components)).astype('float32')
Q2_cv = np.zeros(len(n_components)).astype('float32')
PRESScv = np.zeros(len(n_components)).astype('float32')
RSS_cv = np.zeros(len(n_components)).astype('float32')
WoldR = np.zeros(len(n_components)).astype('float32')
for i in n_components:
rmsec = []
rmsecv = []
PRESS = []
RSS = []
R2 = []
Q2 = []
for train, test in cval.split(X):
PLS = PLSRegression(n_components=i, scale=False)
PLS.fit_transform(X.iloc[train], y.iloc[train])
Y_pred_train = PLS.predict(X.iloc[train])
Y_pred_test = PLS.predict(X.iloc[test])
rmsec.append(np.sqrt(mean_squared_error(y.iloc[train], Y_pred_train[:,0])))
rmsecv.append(np.sqrt(mean_squared_error(y.iloc[test], Y_pred_test[:,0])))
y_value = y.iloc[test]
y_value = y_value.iloc[0]
RSS.append((Y_pred_test - y_value)**2)
PRESS.append(mean_squared_error(y.iloc[test], Y_pred_test[:,0]))
R2Y = r2_score(y.iloc[train], Y_pred_train)
rmse_c[i-1] = np.mean(np.array(rmsec))
rmse_cv[i-1] = np.mean(np.array(rmsecv))
RSS_cv[i-1] = np.mean(np.array(RSS))
PRESScv[i-1] = np.mean(np.array(PRESS))
R2_cv[i-1] = np.mean(np.array(R2Y))
if i == 1:
Q2_cv[0] = 0
else:
Q2_cv[i-1] = (1 - (PRESScv[i-1]/RSS_cv[i-2]))
if i == 1:
WoldR[i-1] = 0
elif i > 1:
WoldR[i-1] = (PRESScv[i-1]/PRESScv[i-2])
with plt.style.context(('seaborn')):
plt.plot(rmse_c, 'o-r', label = "RMSE Calibration")
plt.plot(rmse_cv, 'o-b', label = "RMSE Cross Validation")
plt.ylabel('RMSE')
plt.xlabel('Number of LV included')
plt.legend()
plt.show()
with plt.style.context(('seaborn')):
plt.plot(R2_cv, 'o-r', label = "R2_score")
plt.plot(Q2_cv, 'o-b', label = "Q2_score")
plt.ylabel('Score')
plt.xlabel('Number of LV included')
plt.legend()
plt.show()
with plt.style.context(('seaborn')):
plt.plot(WoldR, 'o-r', label = "Wold R")
plt.plot(PRESScv, 'o-b', label = "PRESS")
plt.ylabel('Score')
plt.xlabel('Number of LV included')
plt.legend()
plt.show()
#%% Select n_components based on LOOCV:
PLS = PLSRegression(n_components=6, scale=True)
PLS.fit_transform(X, y)
Y_pred = PLS.predict(X)
bin_pred = (PLS.predict(X)[:,0] > 0.5).astype('uint8')
plt.figure(dpi=720)
plt.title('Confusion Matrix PLS-DA')
cf_rf = confusion_matrix(y, bin_pred)
sns.heatmap(cf_rf, annot=True, cmap='Blues')
plt.xlabel('Predicted outcomes')
plt.ylabel('True outcomes')
plt.show()
print(classification_report(y, bin_pred))
scores_PLS = precision_recall_fscore_support(y, bin_pred)
#%%
def vip(model):
t = model.x_scores_
w = model.x_weights_
q = model.y_loadings_
p, h = w.shape
vips = np.zeros((p,))
s = np.diag(t.T @ t @ q.T @ q).reshape(h, -1)
total_s = np.sum(s)
for i in range(p):
weight = np.array([ (w[i,j] / np.linalg.norm(w[:,j]))**2 for j in range(h) ])
vips[i] = np.sqrt(p*(s.T @ weight)/total_s)
return vips
VIP = vip(PLS)
#%%
cols = ['Age (years)', 'Gender', 'Daily alcohol use', 'Medication',
'ASA-classification', 'Recurrent disease?', 'Comorb',
'Independent, with others', 'Smokes cigarettes/sigar', 'BMI', 'GFI',
'HADS_A', 'HADS Depression', 'ADL', 'iADL', 'TUG', 'Handgrip strength',
'Avg. Steps/day', 'Avg. MVPA/day']
VIP = pd.DataFrame(VIP, columns=['VIP'])
VIP['labels'] = cols
VIP = VIP.sort_values('VIP', ascending=False)
#%%
BCT = VIP['VIP'].describe()[6] + (1.5*( VIP['VIP'].describe()[6]- VIP['VIP'].describe()[4]))
#%%
VIP = VIP[VIP['VIP']>1]
#%%
import matplotlib.pylab as pylab
params = {'legend.fontsize': 'x-large',
'axes.labelsize': 'x-large',
'axes.titlesize':'x-large',
'xtick.labelsize':'x-large',
'ytick.labelsize':'x-large'}
pylab.rcParams.update(params)
plt.figure(figsize=(12,7), dpi=720)
sns.set_style('whitegrid')
plt.bar(x=VIP['labels'], height=VIP['VIP'])
plt.title('VIP scores')
plt.yticks(fontsize='x-large')
plt.xticks(rotation=-90, fontsize='x-large')
#%% Compute ROC curve and ROC area for each class
fpr = dict()
tpr = dict()
roc_auc = dict()
fpr, tpr, _ = roc_curve(y, bin_pred)
roc_auc = auc(fpr, tpr)
# Compute micro-average ROC curve and ROC area
fpr_m, tpr_m, _ = roc_curve(y.ravel(), bin_pred.ravel())
roc_auc_m = auc(fpr_m, tpr_m)
plt.figure(dpi=720)
lw = 2
plt.plot(fpr, tpr, color='darkorange',
lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)
plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC-Curve PLS-DA')
plt.legend(loc="lower right")
plt.show()
plt.figure(dpi=720)
lw = 2
plt.plot(fpr_m, tpr_m, color='darkorange',
lw=lw, label='ROC curve (area = %0.2f)' % roc_auc_m)
plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic example')
plt.legend(loc="lower right")
plt.show()