<jupyter_start><jupyter_text>Load prediction from CXR-based model<jupyter_code>data = pd.read_csv('/content/gdrive/MyDrive/covid19/predictions_cxr_model_4fold_models_predictions_v2.csv')
data.head()<jupyter_output><empty_output><jupyter_text>Remove unnecessary columns<jupyter_code>data_drop = data.drop(columns=['Unnamed: 0', 'SOPInstanceUID','path','APview','bbox_coordinates',])
data_drop = data_drop.set_index('patient_id')
data_drop = data_drop.sort_index()<jupyter_output><empty_output><jupyter_text>Train-test split Indicate the fold number to determine the train-test set<jupyter_code>fold_num = 0
Test_data = data_drop[(data_drop.fold == fold_num) & (data_drop.index > 0)]
# Train_data = data_drop[(data_drop.fold != fold_num) & (data_drop.index > 0)]

#Train final model
Train_data = data_drop[(data_drop.index > 0)]
# x_test = Test_data.drop(columns=['expired_30_days', 'fold', 'predicted_probability'])
x_test = Test_data.drop(columns=['expired_30_days', 'fold'])
y_test = Test_data['expired_30_days']
# x_train = Train_data.drop(columns=['expired_30_days', 'fold', 'predicted_probability'])
x_train = Train_data.drop(columns=['expired_30_days', 'fold'])
y_train = Train_data['expired_30_days']
print(x_train.shape)
print(x_test.shape)
print(y_train.shape)
print(y_test.shape)
list(x_test.columns)<jupyter_output>(1628, 46)
(414, 46)
(1628,)
(414,)<jupyter_text>Model Selection<jupyter_code>from sklearn.metrics import roc_auc_score, recall_score, make_scorer

def clipped_auc_scorer(y_true,y_pred):
    """
    Calculates the AUROC score, but returns 0 if sensitivity of the positive class is < 0.7
    """
    auc_score = roc_auc_score(y_true, y_pred)
    if recall_score(y_true,y_pred) >= 0.7:
        return auc_score
    else:
        return 0

clipped_auc = make_scorer(clipped_auc_scorer)
def model_selection(summary_dict,model_lst,param_dict,
                    x_train=x_train,y_train=y_train,x_test=x_test,y_test=y_test,
                    n_iter=None,k_fold=None,n_repeats=None):
    """
    A wrapper function for the model selection loop
    
    Parameters
    ----------
    summary_dict: dict
        An empty dictionary used to store results.
    model_lst: list
        A list of tuples containing ('model_name',model), models are sklearn estimators
    param_dict: dict
        A dictionary containing model parameter distributions - to be passed to RandomizedSearchCV
    x_train: array-like
        An array training set predictors
    y_train: array-like
        An array containing training set labels
    x_test: array-like
        An array containing test set predictors
    y_test: array-like
        An array containing test set labels
    n_iter: int
        Number of crossvalidation iterations - to be passed to RandomizedSearchCV. Defaults to n_iter parameter at top of script
    k_fold: int
        Number of crossvalidation folds - to be passed to RandomizedSearchCV. Defaults to k_fold parameter at top of script
    n_repeats: int
        Number of crossvalidation repeats - to be passed to RandomizedSearchCV. Defaults to n_repeats parameter at top of script
        
    Returns
    -------
    summary_dict: pandas.DataFrame
        A dataframe containing the best model object and associated crossvalidation results
    result_table: pandas.DataFrame
        A dataframe containing all model objects and associated crossvalidation results
    """
    iterations = n_iter
    
    #Full list of scoring metrics
    scoring = {'roc_auc':'roc_auc',
               'average_precision':'average_precision',
               'accuracy': 'accuracy',
               'f1_macro':'f1_macro',
               'clipped_auc':clipped_auc}

    #Create an empty list used to store the results
    result_list = []

    #Loop through the list of models
    for name, model in model_lst:
        #Define the cross-validation folds
        cv = RepeatedStratifiedKFold(n_splits=k_fold,n_repeats=n_repeats)
        
        #Set the optimizing metric
        #refit_score = 'clipped_auc'
        refit_score = 'f1_macro'

        
        
        #Set the estimator as the model currently being optimized
        estimator = model
        
        for sample_weight_method in ['unbalanced','balanced']:
            for method in ['sigmoid','none']:
                #Wrap the model into calibration loop
                if method != 'none':
                    clf = CalibratedClassifierCV(estimator,method=method)
                    params = param_dict
                else:
                    clf= estimator
                    new_dict = {}
                    for key in param_dict.keys():
                        subdict = {}
                        for nested_key in param_dict[key].keys():
                            new_key_name = nested_key.split('base_estimator__')[1]
                            subdict[new_key_name] = param_dict[key][nested_key]
                        new_dict[key] = subdict
                    params = new_dict

                #Create the RandomizedSearchCV object
                search = RandomizedSearchCV(clf,param_distributions=params.get(name),random_state=seed,cv=cv,n_iter=iterations,n_jobs=-1,
                                        scoring=scoring,refit=refit_score,verbose=False,return_train_score=False)

                #Begin the grid search process
                if sample_weight_method == 'balanced':
                    search.fit(x_train, y_train, sample_weight=weight_array) 
                else:
                    search.fit(x_train, y_train) 

                #Calculate some metrics on the full training dataset (purely for diagnostics)
                y_pred = search.best_estimator_.predict(x_train)

                print(f'Algorithm: {name}')
                print(f'Calibration method: {method}')
                print(f'Sample weight method: {sample_weight_method}')
                print(classification_report(y_true=y_train,y_pred=y_pred))
                print(f'CV score ({refit_score}): {search.best_score_}')
                print()
            
                #Append the results of the best model to results_list
                result_list.append((name+'_'+method+'_'+sample_weight_method,search,search.best_score_,search.cv_results_))

        ##End of loop
    
    #The following code tidies result_list in to a dataframe
    result_table = pd.DataFrame(result_list,columns=['name','model','scores','score_dict'])
    best_model_index = result_table['scores']==max(result_table['scores'])
    model_name = result_table['name'][best_model_index].values.tolist()[0]
    best_model = result_table['model'][best_model_index].values.tolist()[0]
    summary_dict['results'] = {'Model':model_name}
    metrics = ['mean_test_roc_auc','mean_test_average_precision','mean_test_accuracy']
    if hasattr(best_model,'best_score_'):
        best_score = best_model.best_score_ 
        for key in [key for key in best_model.cv_results_.keys() if key in metrics]:
            summary_dict['results'][key.split('mean_test_')[1]] = best_model.cv_results_[key][best_model.best_index_]
        summary_dict['results']['model obj'] = best_model.best_estimator_
    else:
        best_score = mean_cv_roc_auc
        summary_dict['results']['model obj'] = best_model
        for key in [key for key in best_search.cv_results_.keys() if key in metrics]:
            summary_dict['results'][key.split('mean_test_')[1]] = result_table['score_dict'][best_model_index].get(key)
    
    #Find the overall results
    print(f"Best Cross-Validation score: {best_score}")        
    return summary_dict, result_table
#This cell runs the model selection loop 
summary_dict = {}

#Create model objects
logistic = SGDClassifier(loss='log',random_state=seed)
rf = RandomForestClassifier(random_state=seed)
gbm = GradientBoostingClassifier(random_state=seed)
xgbclf = XGBClassifier(random_state=seed,objective='binary:logistic')

classifier_list = [('xgb',xgbclf),('lr',logistic)] #,('gbm',gbm),('rf',rf),

#Calculate class weights
weight_array = compute_sample_weight(class_weight="balanced",y=y_train)

#Define the hyperparameter search space
params = {'lr':{'base_estimator__alpha':uniform(1e-5,1e-3),
                'base_estimator__penalty':['l1', 'l2', 'elasticnet'],
                'base_estimator__l1_ratio':uniform(0.01,0.30)},
          'rf':{'base_estimator__bootstrap': [True, False],
                'base_estimator__max_depth':randint(3,12),
                'base_estimator__max_features': ['auto', 'sqrt'],
                'base_estimator__min_samples_split':randint(2,12),
                'base_estimator__min_samples_leaf':randint(2,12),
                'base_estimator__n_estimators':randint(200, 1000)},
          'gbm':{'base_estimator__loss':['deviance','exponential'],
                 'base_estimator__learning_rate':uniform(0.003, 0.3),
                 'base_estimator__n_estimators':randint(200, 1000),
                 'base_estimator__subsample':uniform(0.5, 0.5),
                 'base_estimator__criterion':['friedman_mse','mse','mae'],
                 'base_estimator__min_samples_split':randint(2,12),
                 'base_estimator__min_samples_leaf':randint(2,12),
                 'base_estimator__max_depth':randint(3,12),
                 'base_estimator__max_features':['sqrt', 'log2']},
          'xgb':{'base_estimator__colsample_bytree': uniform(0.5,0.5),
                 'base_estimator__eta': (0.0001,0.1),
                 'base_estimator__max_depth': randint(3,12),
                 'base_estimator__min_child_weight': randint(3,12),
                 'base_estimator__subsample': uniform(0.5,0.5)}}

#Run the model selection loop
summary_dict, results = model_selection(summary_dict=summary_dict,
                                        model_lst=classifier_list,
                                        param_dict=params,
                                        n_iter=20,k_fold=5,n_repeats=1)
best_model_index = results['scores']==max(results['scores'])
best_params = results['model'][best_model_index].values.tolist()[0]
print(best_params)<jupyter_output><empty_output><jupyter_text>Evaluation<jupyter_code>from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
%matplotlib inline
cv_results = results.sort_values('scores',ascending=False)

clf_list = list(zip(cv_results['name'].iloc[0:4],cv_results['model'].iloc[0:4]))
plt.figure(figsize=(6,5),dpi=100)
plt.title('Receiver Operating Characteristic')

color_list = ['#581845','#900C3F','#C70039','#FF5733']
col = 0

from sklearn import metrics

for clf_name,clf in clf_list:
    clf = clf.best_estimator_
    y_prob = clf.predict_proba(x_test)
    prob = y_prob[:,1]
    fpr, tpr, threshold = roc_curve(y_test, prob)
    roc_auc = metrics.auc(fpr, tpr)
    plt.plot(fpr, tpr, color_list[col], label = f'{clf_name} = %0.2f' % roc_auc)
    plt.legend(loc = 'lower right')
    plt.plot([0, 1], [0, 1],'k--')
    plt.xlim([0, 1])
    plt.ylim([0, 1])
    plt.ylabel('True Positive Rate')
    plt.xlabel('False Positive Rate')
    col += 1
plt.show();
import os
# Save the best 4 estimaters
model_path = './gdrive/My Drive/covid19/models/'
if not os.path.isdir(model_path):
    os.mkdir(model_path)

count = 1       
for name,best_model in clf_list:
    best_model = best_model.best_estimator_
    y_score = best_model.predict_proba(x_test)[:,1]
    average_precision = average_precision_score(y_test, y_score)
    disp = plot_precision_recall_curve(best_model, x_test, y_test)
    disp.ax_.set_title(f'{name} 2-class Precision-Recall curve: '
                    'AP={0:0.2f}'.format(average_precision))
    # joblib.dump(best_model, os.path.join(model_path,name+'_bestf1_'+str(count)+'_fold' + str(fold_num) + '.pkl'), compress = 1)
    # joblib.dump(best_model, os.path.join(model_path,name+'_str_only_bestf1'+str(count)+'_fold' + str(fold_num) + '.pkl'), compress = 1)
    joblib.dump(best_model, os.path.join(model_path,name+'_str_only_bestf1'+str(count)+'.pkl'), compress = 1)
    count = count + 1<jupyter_output><empty_output><jupyter_text>Compute final result<jupyter_code>import glob

#reading back in

for fold_num in range(4):
  fpath = glob.glob(model_path+'*_str_only_bestf11_fold'+str(fold_num)+'.pkl') 
  # fpath = glob.glob(model_path+'*_bestf1_1_fold'+str(fold_num)+'.pkl') 

  clf = joblib.load(fpath[0]) 
  ypred = clf.predict(x_test)
  print(classification_report(y_true=y_test,y_pred=ypred))
  print(clf.score)

  if (fold_num == 0):
    y_train_0 = y_test.copy()
    ypred_0 = ypred.copy()
    y_prob_0 = y_prob.copy()
  if (fold_num == 1):
    y_train_1 = y_test.copy()
    ypred_1 = ypred.copy()
    y_prob_1 = y_prob.copy()
  if (fold_num == 2):
    y_train_2 = y_test.copy()
    ypred_2 = ypred.copy()
    y_prob_2 = y_prob.copy()
  if (fold_num == 3):
    y_train_3 = y_test.copy()
    ypred_3 = ypred.copy()
    y_prob_3 = y_prob.copy() 
y_test_c = np.concatenate((y_train_0, y_train_1, y_train_2, y_train_3), axis=0)
ypred_c = np.concatenate((ypred_0, ypred_1, ypred_2, ypred_3), axis=0)
y_prob_c = np.concatenate((y_prob_0, y_prob_1, y_prob_2, y_prob_3), axis=0)<jupyter_output><empty_output><jupyter_text>Save the prediction results<jupyter_code># np.savetxt("./gdrive/My Drive/covid19/f1_y_test_c_fuse.csv", y_test_c, delimiter=",")
# np.savetxt("./gdrive/My Drive/covid19/f1_ypred_c_fuse.csv", ypred_c, delimiter=",")
# np.savetxt("./gdrive/My Drive/covid19/f1_y_prob_c_fuse.csv", y_prob_c, delimiter=",")

np.savetxt("./gdrive/My Drive/covid19/f1_y_test_c_str_only.csv", y_test_c, delimiter=",")
np.savetxt("./gdrive/My Drive/covid19/f1_ypred_c_str_only.csv", ypred_c, delimiter=",")
np.savetxt("./gdrive/My Drive/covid19/f1_y_prob_c_str_only.csv", y_prob_c, delimiter=",")
from numpy import genfromtxt

y_test_c = genfromtxt("./gdrive/My Drive/covid19/y_test_c.csv", delimiter=',')
ypred_c = genfromtxt("./gdrive/My Drive/covid19/ypred_c.csv", delimiter=',')
y_prob_c = genfromtxt("./gdrive/My Drive/covid19/y_prob_c.csv", delimiter=',')

# y_test_c = genfromtxt("./gdrive/My Drive/covid19/f1_y_test_c_fuse.csv", delimiter=',')
# ypred_c = genfromtxt("./gdrive/My Drive/covid19/f1_ypred_c_fuse.csv", delimiter=',')
# y_prob_c = genfromtxt("./gdrive/My Drive/covid19/f1_y_prob_c_fuse.csv", delimiter=',')

# y_test_c = genfromtxt("./gdrive/My Drive/covid19/f1_y_test_c_str_only.csv", delimiter=',')
# ypred_c = genfromtxt("./gdrive/My Drive/covid19/f1_ypred_c_str_only.csv", delimiter=',')
# y_prob_c = genfromtxt("./gdrive/My Drive/covid19/f1_y_prob_c_str_only.csv", delimiter=',')

prob = y_prob_c[:,1]<jupyter_output><empty_output><jupyter_text>Draw ROC curve for CXR-based model<jupyter_code># Test_data_all = data_drop[(data_drop.index > 0)]
# prob=Test_data_all['predicted_probability']
# y_test_c = Test_data_all['expired_30_days']<jupyter_output><empty_output><jupyter_text>Draw ROC curve for EHR-based and fusion model<jupyter_code>plt.figure(figsize=(6,5),dpi=100)
plt.title('Receiver Operating Characteristic')

color_list = ['#581845','#900C3F','#C70039','#FF5733']
col = 0

from sklearn import metrics

# for clf_name,clf in clf_list:
#     clf = clf.best_estimator_
# y_prob = clf.predict_proba(x_test)

fpr, tpr, threshold = roc_curve(y_test_c, prob)
roc_auc = metrics.auc(fpr, tpr)
plt.plot(fpr, tpr, color_list[col], label = f'AUC = %0.2f' % roc_auc)
plt.legend(loc = 'lower right')
plt.plot([0, 1], [0, 1],'k--')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
col += 1
plt.savefig("output.svg")
plt.show();

print(classification_report(y_true=y_test_c,y_pred=ypred_c))
import sklearn.metrics

fpr, tpr, thresholds = roc_curve(y_test_c, prob)
optimal_idx = np.argmin((1-tpr)**2+ (fpr)**2 )
optimal_threshold = thresholds[optimal_idx]

## use optimal cut-point 
predicted_class = np.zeros(len(prob), dtype=int)
predicted_class[prob < optimal_threshold] = 0
predicted_class[prob >= optimal_threshold] = 1

## use 0.5 as cut-point 
# predicted_class[predictions < 0.5] = 0
# predicted_class[predictions >= 0.5] = 1
# predicted_class_name = np.array(CLASSES)[predicted_class].tolist()

# use softmax
# predicted_class = np.around(predictions).astype(int)
# predicted_class_name = np.array(CLASSES)[predicted_class].tolist()


tn, fp, fn, tp = sklearn.metrics.confusion_matrix(y_test_c, predicted_class).ravel()
roc_auc = sklearn.metrics.roc_auc_score(y_test_c, prob)
sensitivity = tp / (tp + fn) # recall
specificity = tn / (tn + fp)
PPV = tp / (tp + fp) #precision
NPV = tn / (tn + fn) #NPV
f1 = 2 * ((tp / (tp + fp)) * (tp / (tp + fn)) )/ ((tp / (tp + fp)) + (tp / (tp + fn)))
acc = (tp + tn) / (tn + fp + fn + tp)

print(f'AUROC: {roc_auc:.2f}')
print(f'Sensitivity: {sensitivity:.2f}')
print(f'Specificity: {specificity:.2f}')
print(f'PPV: {PPV:.2f}')
print(f'NPV: {NPV:.2f}')    
print(f'F1-score: {f1:.2f}')
print(f'Accuracy: {acc:.2f}')
print(f'tn, fp, fn, tp: {tn, fp, fn, tp}')<jupyter_output>AUROC: 0.82
Sensitivity: 0.83
Specificity: 0.73
PPV: 0.30
NPV: 0.97
F1-score: 0.44
Accuracy: 0.75
tn, fp, fn, tp: (1048, 380, 32, 160)
{'fraction_non_alphanumeric': 0.13158819395643007, 'fraction_numerical': 0.02166783790114781, 'mean_word_length': 3.733093126385809, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '568127', 'license': 'MIT', 'n_tokens_mistral': 6355, 'n_tokens_neox': 6031, 'n_words': 1067, 'path': '3.training/madrid_structured_and_fusion_4fold.ipynb', 'repo_name': 'theonesp/fusion_model_mortality_covid', 'size': '17076'}