I am writing a machine learning model to predict metastases in patients with gastric cancer based on radiomic features. I am using StratifiedKfold for splitting the training- and test set, and then I do some processing on the data such as finding the most prominent features. For the feature reduction I have used cross validation, but I have failed to implement this in the classifier part of the model. The following is my code:
from sklearn.cluster import FeatureAgglomeration
from sklearn.model_selection import GridSearchCV, train_test_split
from sklearn.linear_model import Lasso
from numpy import mean
from numpy import std
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import Pipeline
from sklearn import metrics
from sklearn.svm import SVC
from sklearn.metrics import classification_report, confusion_matrix
import matplotlib.pyplot as plt
df = pd.read_excel(r'H:\.xlsx') # Load in dataset
X = df.iloc[:,1:] # Radiomic features
y = df.iloc[:,0] # Metastases (1) or not (0)
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=40) # Define StratifiedKFold with 5 splits
for train_index, test_index in skf.split(X, y):
X_train, X_test = X.loc[train_index], X.loc[test_index] #Create train- and test set for radiomic features, stratify based on the ratio of metastases
y_train, y_test = y.loc[train_index], y.loc[test_index] #Create train- and test set for metastases, stratify based on the ratio of metastases
X_train_radiomic=pd.DataFrame(scaler.fit_transform(X_train.iloc[:,6:]))
X_test_radiomic = pd.DataFrame(scaler.fit_transform(X_test.iloc[:,6:]))
rank = X_train_radiomic.corr(method='spearman') # performs Spearman rank correlation on radiomic features, discards first column with metastases
agglo=FeatureAgglomeration(n_clusters=40).fit(rank) # clusters the features into 40 clusters and fits to Spearman rank correlation
reduced_features=[]
for i, label in enumerate(set(agglo.labels_)): # gives overview of what features belong to which cluster
features_with_label = [j for j, lab in enumerate(agglo.labels_) if lab == label]
features_with_label = features_with_label[0]
reduced_features.append(radiomic_features[features_with_label])
print('Features in agglomeration {}: {}'.format(i, features_with_label))
X_trainred=pd.DataFrame(agglo.transform(X_train_radiomic)) # Transforms radiomic features based on clustered features
X_trainred.columns=reduced_features
X_testred=pd.DataFrame(agglo.transform(X_test_radiomic)) # Transforms radiomic features based on clustered features
X_testred.columns=reduced_features
def classifying(classifier,type_model): # function for classifying the model by RF, LR and SVM for each category
total_AUC = 0 # preallocate values that will be used later
total_accuracy = 0
total_fpr = 0
total_tpr = 0
fpr_nr = 0
tpr_nr = 0
if classifier == "RF": #define what classifier to use
algorithm=RandomForestClassifier()
parameters={'classifier__n_estimators':[11]}
if classifier == "LR":
algorithm=LogisticRegression()
parameters={'classifier__solver':['liblinear'] , 'classifier__max_iter':[100]}
if classifier == "SVM":
algorithm=SVC()
parameters={'classifier__kernel':['linear'],'classifier__probability':[True]}
pipeline = Pipeline([('scaler',scaler),('model', Lasso())]) #create Pipeline for LASSO algorithm (data is already scaled so no more scaling needed)
search = GridSearchCV(pipeline, {'model__alpha': np.arange(0.001, 0.01, 0.001)}, cv=5, scoring="neg_mean_squared_error",
verbose=2, refit=True, return_train_score=True)
search.fit(X_trainred, y_train) #fit estimators on training data
coefficients = search.best_estimator_.named_steps['model'].coef_ #estimator which gave highest score
importance = np.abs(coefficients) # absolute value of coefficient, importance of the features for the model
discard=np.where(importance == 0)
X_trainred.drop(X_trainred.columns[discard],axis=1, inplace=True)
X_testred.drop(X_testred.columns[discard], axis=1, inplace=True)
X_train_full = X_trainred.reset_index(drop=True)
X_test_full = X_testred.reset_index(drop=True)
X_train_clinical = pd.DataFrame(X_train.iloc[:, :5]).reset_index(drop=True)
X_test_clinical = pd.DataFrame(X_test.iloc[:, :5]).reset_index(drop=True)
X_train_full=pd.concat([X_train_clinical, X_train_full], axis=1)
X_test_full=pd.concat([X_test_clinical, X_test_full], axis=1)
pipeline=Pipeline([('classifier', algorithm)])
clf=GridSearchCV(pipeline, parameters)
if type_model == 'clinical':
X_train_full = X_train_clinical
X_test_full = X_test_clinical
if type_model == 'radiomics':
X_train_full = X_trainred
X_test_full = X_testred
if type_model == 'combined':
X_train_full = X_train_full
X_test_full = X_test_full
y_predict_proba = clf.fit(X_train_full, y_train).predict_proba(X_test_full) #calculate probabilities for metastases or not (two column values)
y_predict=clf.fit(X_train_full, y_train).predict(X_test_full) #predict metastases based on estimators
fpr, tpr, _ = metrics.roc_curve(y_test, y_predict_proba[:,1]) # calculate false positive rate (fpr) and true positive rate (tpr) using ROC curve
total_fpr=total_fpr+np.sum(fpr) # calculate total fpr over all folds
total_tpr=total_tpr+np.sum(tpr) # calculate total tpr over all folds
fpr_nr=fpr_nr+len(fpr) # calculate the number of gives fpr's per fold, and add to total number of fpr's
tpr_nr = tpr_nr + len(tpr) # calculate the number of gives tpr's per fold, and add to total number of tpr's
total_AUC=total_AUC+metrics.roc_auc_score(y_test, y_predict_proba[:,1]) # calculate total AUC for all folds
total_accuracy=total_accuracy+metrics.accuracy_score(y_test,y_predict, normalize=True) # calculate total accuracy for all folds
mean_AUC=total_AUC/5 # mean AUC of all folds
mean_accuracy=total_accuracy/5 # mean accuracy of all folds
specificity=1-(total_fpr/fpr_nr) # specificity of all folds
sensitivity=(total_tpr/tpr_nr) # sensitivity of all folds
print(
f"The type of model is {type_model}\n"
f"The classifying model is {algorithm}\n"
f"The mean AUC is {total_AUC}\n"
f"The mean accuracy is {total_accuracy}\n"
f"The specificity is {specificity}\n"
f"The sensitivity is {sensitivity}\n"
)
clinical_RF=classifying('RF','clinical')
I have tried to implement the cross validation by using the indices of the trainingset, but since I dropped them to combined two dataframe, this does not work.
The first row of my dataset looks like this (I can add dictionary tomorrow).
The first row is the output (e.g. y_train/test), the first five are the clinical features which are categorical, and from the sixth column on there are radiomic features.
How can I implement to also use the cross validation with my training set in classifying the model (RF/LR/SVM)?