mortality_prediction.py

# -*- coding: utf-8 -*-
"""Copia de mortality_glm.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1PS_-Ls5AZlTjqV6cLAZKfcmSEAAsWCYU
"""

# Import libraries
import numpy as np
import os
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.path as path
from sklearn.metrics import accuracy_score, recall_score, precision_score
from sklearn.metrics import f1_score, confusion_matrix, roc_auc_score
from sklearn.metrics import average_precision_score, roc_curve, precision_recall_curve
from inspect import signature
import statsmodels.api as sm
from sklearn.calibration import calibration_curve, CalibratedClassifierCV

# Make pandas dataframes prettier
from IPython.display import display, HTML

import math
import seaborn as sns

from google.colab import drive
drive.mount('/content/gdrive')
base_dir = "/content/gdrive/My Drive/SCCM_Team5_Glucose_TBI"

df_w_glu = pd.read_csv(base_dir+'/df_for_model_plus_deltas.csv')

display(df_w_glu.head())

from sklearn.pipeline import Pipeline
from sklearn import impute, metrics
from sklearn import preprocessing 
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression


#predictors = ['age','apachescore','diabetes_flag','delta_glucose_mean','Glucose_mean']
predictors = ['apachescore','diabetes_flag','delta_glucose_mean','Glucose_mean']
# predictors = ['diabetes_flag','delta_glucose_mean','Glucose_mean']

outcome = 'Day1_Mortality'

df_day1 = df_w_glu[df_w_glu['DAY'] == 1].copy()

X = df_day1[predictors]
y = df_day1[outcome]

display(HTML('<H3>Dataframe size<H3>'))
print(f'Predictors: {X.shape[1]}')
print(f'Entries: {X.shape[0]}')
display(HTML('<H3>In-hospital mortality within 24 hours in the ICU<H3>'))
print(f'Class 1 (died): {sum(y)}, {(sum(y)/len(y))*100:0.2f}%')
print(f'Class 0 (alive): {len(y)- sum(y)}')

# Complete diabetes flag
X['diabetes_flag'].fillna(value=0,inplace=True)

# Find missing data
var_nan = X.columns[X.isnull().any()]
print(f'Variables with at least one missing value {list(var_nan)}')

plt.figure()
sns.heatmap(X.isnull(), cbar=False)
plt.xlabel(f'Features n={len(X.columns)}')
plt.ylabel('Observations or lines')
plt.title(f'Missingness n={X.shape[0]}')
fig = plt.gcf()
fig.set_size_inches(7, 7)

"""## Data description
For table 1
"""

display(df_day1[['age','apachescore','delta_glucose_mean','Glucose_mean']].describe())

df_day1.columns

!pip install tableone
from tableone import TableOne

df_day1['diabetes_flag'].fillna(value=0,inplace=True)

categorical = ['gender','ethnicity','diabetes_flag',
               'surgery']
groupby = 'hospitaldischargestatus'
columns = ['gender', 'age','ethnicity', 'diabetes_flag',
           'apachescore','Glucose_mean','delta_glucose_mean',
           'surgery']
glu_table = TableOne(df_day1, groupby = groupby,
                     columns = columns, categorical = categorical,
                     pval= True)
print(glu_table.tabulate(tablefmt="github"))

display(HTML('<H4>Brain problem counts<H4>'))
print(f'AIS: {sum(df_day1.brain_problem=="AIS")}')
print(f'HEM: {sum(df_day1["brain_problem"]=="HEM")}')
print(f'Other: {sum(df_day1["brain_problem"]=="Other")}')
print(f'TBI: {sum(df_day1["brain_problem"]=="TBI")}')
print(f'SZ: {sum(df_day1["brain_problem"]=="SZ")}')
print(f'''Total: {sum(df_day1.brain_problem=="AIS")+ 
sum(df_day1["brain_problem"]=="HEM")+
sum(df_day1["brain_problem"]=="Other")+
sum(df_day1["brain_problem"]=="TBI")+
sum(df_day1["brain_problem"]=="SZ")}''')

"""## Modeling"""

# Partition
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, 
                                                    random_state=10,
                                                    stratify = y)

display(HTML('<H3>In-hospital mortality within 24 hours in the ICU<H3>'))
display(HTML('<H4>Training<H4>'))
print(f'Total entries: {len(y_train)}')
print(f'Class 1 (died): {sum(y_train)}, {(sum(y_train)/len(y_train))*100:0.2f}%')
print(f'Class 0 (alive): {len(y_train)- sum(y_train)}')
display(HTML('<H4>Test<H4>'))
print(f'Total entries: {len(y_test)}')
print(f'Class 1 (died): {sum(y_test)}, {(sum(y_test)/len(y_test))*100:0.2f}%')
print(f'Class 0 (alive): {len(y_test)- sum(y_test)}')

# upsampling mortality
display(HTML('<H4>Before upsampling<H4>'))
train_class0 = np.where(y_train == 0)[0]
print(f'Class 0 for training (original): {len(y_train.iloc[train_class0])}')
train_class1 = np.where(y_train == 1)[0]
print(f'Class 1 for training (original): {len(y_train.iloc[train_class1])}')
train_class1_upsampled = np.random.choice(train_class1, size=len(train_class0), 
                                          replace=True)
display(HTML('<H4>Ramdomly upsampling Class 1 (length)<H4>'))
print(f'Class 1 for training (up-sampled): {len(train_class1_upsampled)}')

y_train = np.concatenate((y_train.iloc[train_class0], y_train.iloc[train_class1_upsampled]))
X_train = np.concatenate((X_train.iloc[train_class0], X_train.iloc[train_class1_upsampled]))

# Review the number of cases in each set
display(HTML('<H4>Total entries<H4>'))
print("Entries for training (upsampled): {}".format(X_train.shape[0]))
print("Entries for testing (upsampled): {}".format(X_test.shape[0]))

def specificity_score(y_true, y_pred):
  # TN/TN+FP
  # Confusion matrix
  matrix = confusion_matrix(y_true, y_pred)
  spec = matrix[0,0]/(matrix[0,0]+matrix[0,1])
  return spec

model = LogisticRegression(solver='lbfgs', n_jobs= -1,random_state = 44, C = 0.5)

# Impute missing values and scale using a pipeline
estimator = Pipeline([("imputer", impute.SimpleImputer(missing_values=np.nan, strategy="mean")),
                      ("scaler", preprocessing.StandardScaler()),
                      ("logistic_regression", model)])

estimator.fit(X_train, y_train)

y_pred = estimator.predict(X_test)
prob_pos = estimator.predict_proba(X_test)

# More metrics
auc = roc_auc_score(y_test, y_pred)
aupr = average_precision_score(y_test, y_pred)
acc = accuracy_score(y_test, y_pred)
sens = recall_score(y_test, y_pred) # TP/(TP+FN)
spec = specificity_score(y_test, y_pred) # TN/TN+FP
pr = precision_score(y_test, y_pred) # TP/(TP+FP)
f1_sc = f1_score(y_test, y_pred)

display(HTML('<H4>Performance metrics of logistic regression classifier on the test set<H4>'))
print(f'AUC: {auc:0.2f}\nAUCPr: {aupr:0.2f}\nAccuracy: {acc:0.2f}')
print(f'Sensitivity (TPR): {sens:0.2f}\nSpecificity (TNR): {spec:0.2f}')
print(f'Precision: {pr:0.2f}\nF1 score: {f1_sc:0.2f}')

display(HTML('<H4>Summary report<H4>'))
print(metrics.classification_report(y_test, y_pred))

# Confusion matrix
log_conf_matrix = confusion_matrix(y_test, y_pred)

display(HTML('<H4>Confusion matrix<H4>'))
print(f'True positives {log_conf_matrix[1,1]}')
print(f'True negatives {log_conf_matrix[0,0]}')
print(f'False positives {log_conf_matrix[0,1]}')
print(f'False negatives {log_conf_matrix[1,0]}')

# Axis for AUCROC
fpr, tpr, thresholds = roc_curve(y_test, prob_pos[:,1])

# Style
sns.set_style("white",{'xtick.bottom': True,'xtick.top': False,
                'ytick.left': True, 'ytick.right': False,})
sns.set_palette(sns.color_palette("Set2"))

plt.figure()
plt.plot(fpr, tpr, label='Logistic Regression (AUCROC = %0.2f)' % auc)
plt.plot([0, 1], [0, 1], linestyle=':', label='Random classifier', color = 'red')
plt.xlim([-0.01, 1.0])
plt.ylim([0.0, 1.02])

plt.xlabel('False Positive Rate (1-specificity)', fontsize = 16)
plt.xticks(fontsize = 14)
plt.ylabel('True Positive Rate (recall)', fontsize = 16)
plt.yticks(fontsize = 14)
plt.title('Area under the ROC curve', fontsize = 18)
plt.legend(loc="lower right",fontsize = 14)
fig = plt.gcf()
fig.set_size_inches(10, 8)
fig.savefig(base_dir+'/AUROC.png')

precision, recall, thresholds = precision_recall_curve(y_test, 
                                                       y_pred, 
                                                       pos_label=1, 
                                                       sample_weight=None)

# Ref line
sum_wpos = sum(y_test[y_test == 1])
print(f'Positives during testing: {sum_wpos}')
sum_wneg = len(y_test)-sum_wpos
print(f'Negatives during testing: {sum_wneg}')

step_kwargs = ({'step': 'post'}
               if 'step' in signature(plt.fill_between).parameters
               else {})

plt.figure()
plt.step(recall, precision, alpha=0.8,
         where='post',color=sns.color_palette("Set2")[2])
plt.fill_between(recall, precision, alpha=0.3, 
                 color=sns.color_palette("Set2")[2],
                 **step_kwargs, 
                 label = f'AUCPRC={aupr:0.2f}')
plt.xlabel('Recall',fontsize = 16)
plt.ylabel('Precision',fontsize = 16)
plt.ylim([0.0, max(precision)])
plt.yticks(fontsize = 14)
plt.xlim([-0.02, 1.0])
plt.xticks(fontsize = 14)
plt.hlines(sum_wpos/(sum_wneg+sum_wpos),0,1,colors='r',linestyles='dashed',
           label=f'Baseline = {sum_wpos/(sum_wneg+sum_wpos):0.2f}')
plt.legend(fontsize = 14)
plt.title('Precision-Recall plot', fontsize = 18)
fig = plt.gcf()
fig.set_size_inches(10, 8)
fig.savefig(base_dir+'/AUCPr.png')

sns.set_style("darkgrid",{'xtick.bottom': True,'xtick.top': False,
                'ytick.left': True, 'ytick.right': False,})

plt.figure(figsize=(10, 10))
ax1 = plt.subplot2grid((3, 1), (0, 0), rowspan=2)
ax2 = plt.subplot2grid((3, 1), (2, 0))

ax1.plot([0, 1], [0, 1], "k:", label="Perfectly calibrated")

clf = CalibratedClassifierCV(model, cv = 2, method='isotonic')
clf.fit(prob_pos[:,1].reshape(-1, 1), y_test)
prob_pos_test = clf.predict_proba(prob_pos[:,1].reshape(-1, 1))[:,1]
fraction_of_positives, mean_predicted_value = \
        calibration_curve(y_test, prob_pos_test, n_bins=5, normalize=True)

ax1.plot(mean_predicted_value, fraction_of_positives, "--",
             label="%s" % ('LR'))

ax2.hist(prob_pos_test, range=(0, 1), bins=10, label='LR',
             histtype="step", lw=2)

ax1.set_ylabel("Fraction of positives", fontsize=16)
ax1.set_ylim([-0.02, 1.02])
ax1.set_xlim([-0.02, 1.02])
ax1.legend(loc='upper left', facecolor='white', fontsize = 12)
ax1.set_title('Calibration plots  (reliability curve)', fontsize=18)
ax1.set_xlabel("Mean predicted value", fontsize=16)

ax2.set_xlabel("Positive probabilities", fontsize=16)
ax2.set_ylabel("Count", fontsize=16)
ax2.set_xlim([-0.02, 1.02])
ax2.grid(False)
ax2.set_facecolor("white")
ax2.spines['bottom'].set_color('gray')
ax2.spines['left'].set_color('gray')
ax2.spines['top'].set_color('gray')
ax2.spines['right'].set_color('gray')
ax2.legend(loc="upper right", ncol=1, facecolor='white',fontsize = 12)

plt.tight_layout()

imp_mean = impute.SimpleImputer(missing_values=np.nan, strategy="mean")
imp_mean.fit(X_train)
X_imp = imp_mean.transform(X_train)
scaler = preprocessing.StandardScaler()
scaler.fit(X_imp)
X_new = scaler.transform(X_imp)

glm_binom = sm.GLM(y_train, X_new, family=sm.families.Binomial())
res = glm_binom.fit(method='lbfgs')
print(res.summary())

display(HTML('<H4>Adjusted Odds Ratio (AOR)<H4>'))
for x in range(res.params.size):
  print(f'{X_test.columns[x]}: {math.exp(res.params[x]):0.2f}')

display(HTML('<H4>t-values<H4>'))
for x in range(res.params.size):
  print(f'{X_test.columns[x]}: {res.tvalues[x]:0.2f}')

display(HTML('<H4>p-values<H4>'))
for x in range(res.params.size):
  if ((res.pvalues[x] < 0.001) & (res.pvalues[x] != 0)):
    print(f'{X_test.columns[x]}: < 0.001')
  else:
    print(f'{X_test.columns[x]}: {res.pvalues[x]:0.2f}')

display(HTML('<H4>95% confidence interval for AOR<H4>'))
for x in range(res.params.size):
  print(f'{X_test.columns[x]}: {math.exp(res.conf_int()[x,0]):0.2f} to {math.exp(res.conf_int()[x,1]):0.2f}')

# separate by diagnosis. may not have enough datapoints in mortality class if number of days is small
from sklearn.pipeline import Pipeline
from sklearn import impute
from sklearn import preprocessing

predictors = ['age','apachescore','diabetes_flag','delta_glucose_mean','Glucose_mean']
outcome = 'Day1_Mortality'
df_day1 = df_w_glu[df_w_glu['DAY'] == 1]

for prob in set(df_day1['brain_problem']):
  print(prob)
  df_subset = df_day1[df_day1['brain_problem'] == prob]
  # print(df_subset)
  X = df_subset[predictors]
  y = df_subset[outcome]
  X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=10)
  # upsampling mortality
  train_class0 = np.where(y_train == 0)[0]
  train_class1 = np.where(y_train == 1)[0]
  train_class1_upsampled = np.random.choice(train_class1, size=len(train_class0), replace=True)

  y_train = np.concatenate((y_train.iloc[train_class0], y_train.iloc[train_class1_upsampled]))
  X_train = np.concatenate((X_train.iloc[train_class0], X_train.iloc[train_class1_upsampled]))

  # Review the number of cases in each set
  print("Train data: {}".format(len(X_train)))
  print("Test data: {}".format(len(X_test)))
  model = LogisticRegression(solver='lbfgs', n_jobs= -1 )

  # Impute missing values and scale using a pipeline
  estimator = Pipeline([("imputer", impute.SimpleImputer(missing_values=np.nan, strategy="mean")),
                        ("scaler", preprocessing.StandardScaler()),
                        ("logistic_regression", model)])

  estimator.fit(X_train, y_train)

  y_pred = estimator.predict(X_test)
  print('Accuracy of logistic regression classifier on the test set: {:.2f}'.format(estimator.score(X_test.iloc[y_test], y_test)))
  print(metrics.classification_report(y_test, y_pred))

  logit_roc_auc = metrics.roc_auc_score(y_test, estimator.predict(X_test))
  fpr, tpr, thresholds = metrics.roc_curve(y_test, estimator.predict_proba(X_test)[:,1])

  plt.figure()
  plt.plot(fpr, tpr, label='Logistic Regression (area = %0.2f)' % logit_roc_auc)
  plt.plot([0, 1], [0, 1],'r--')
  plt.xlim([0.0, 1.0])
  plt.ylim([0.0, 1.05])

  plt.xlabel('False Positives')
  plt.ylabel('True Positives')
  plt.title('Receiver operating characteristic')
  plt.legend(loc="lower right")
  plt.savefig('Log_ROC')
  plt.show()