EHR_main_TimeSplit.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import re\n",
    "from tempfile import mkdtemp\n",
    "import hashlib\n",
    "import time\n",
    "\n",
    "print(sys.version_info)\n",
    "\n",
    "home_dir = os.getenv(\"HOME\")\n",
    "print(os.getenv(\"PYTHONPATH\"))\n",
    "\n",
    "import numpy as np\n",
    "import pickle\n",
    "from numpy.linalg import svd\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import sklearn\n",
    "from sklearn import linear_model\n",
    "from sklearn.model_selection import train_test_split, GroupKFold\n",
    "from sklearn.linear_model import LogisticRegression, ElasticNet, ElasticNetCV, ARDRegression, SGDClassifier\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn import metrics, tree\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.preprocessing import StandardScaler, Imputer, LabelEncoder\n",
    "from sklearn.metrics import precision_recall_curve, average_precision_score, brier_score_loss, make_scorer\n",
    "from sklearn.model_selection import cross_validate, cross_val_predict, GridSearchCV, StratifiedKFold\n",
    "from sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier\n",
    "from sklearn.calibration import calibration_curve, CalibratedClassifierCV\n",
    "from sklearn.externals import joblib\n",
    "from sklearn.base import BaseEstimator, TransformerMixin\n",
    "from scipy.stats import mstats\n",
    "import scipy.stats as st\n",
    "from IPython.display import display, HTML\n",
    "\n",
    "from tqdm import tqdm_notebook\n",
    "from fancyimpute import SoftImpute, KNN, MICE\n",
    "from ehr_utils import *\n",
    "from imblearn.over_sampling import SMOTE, ADASYN\n",
    "from imblearn.under_sampling import NearMiss\n",
    "from imblearn.pipeline import Pipeline\n",
    "from xgboost import XGBClassifier, XGBRegressor, plot_tree\n",
    "from missingpy import MissForest\n",
    "#import shap\n",
    "#from predictive_imputer import predictive_imputer\n",
    "\n",
    "# from keras.models import model_from_json, Sequential\n",
    "# from keras.layers import Dense, Dropout\n",
    "# from keras.wrappers.scikit_learn import KerasClassifier\n",
    "# from keras.optimizers import SGD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pandas2ri.activate()\n",
    "%matplotlib inline\n",
    "pd.set_option('display.max_rows', 2000)\n",
    "pd.set_option('display.max_columns', 100)\n",
    "# colours for plotting for BJA\n",
    "COLOUR_1 = \"#006BFA\" # blue\n",
    "COLOUR_2 = \"#378E00\" # green\n",
    "COLOUR_3 = \"#E64DFF\" # pink\n",
    "COLOUR_4 = \"#FF8000\" # gold\n",
    "COLOUR_5 = \"#33FFFF\" # light blue\n",
    "COLOUR_6 = \"#FF3300\" # orange\n",
    "PLOT_COLOURS = [COLOUR_1, COLOUR_2, COLOUR_3, COLOUR_4, COLOUR_5, COLOUR_6]\n",
    "# DPI for figures\n",
    "fig_dpi = 1200\n",
    "\n",
    "data_dir=\"/data/DataNov21_2017\"\n",
    "ccs_icd_mapping_f=\"/home/blhill/external_data/ccs_dx_icd10cm_2018_1.csv\"\n",
    "\n",
    "# main_f = os.path.join(data_dir, \"Main_Data.txt\")\n",
    "# outcomes_f = os.path.join(data_dir, \"Outcomes_Data.txt\")\n",
    "# icd_codes_f = os.path.join(data_dir, \"ICD9_10Codes.txt\")\n",
    "# main_filtered_f = \"Main_Data_filtered.txt\"\n",
    "# main_filtered_working_f = \"Main_Data_filtered.working.txt\"\n",
    "asa_predictions_train_f = \"paper/asa_imp/imputed_asa_train.txt\"\n",
    "asa_predictions_test_f = \"paper/asa_imp/imputed_asa_test.txt\"\n",
    "comorbid_f = \"charlson_scores.txt\"\n",
    "pospom_f = \"/opt/data/workingdir/blhill/POSPOM_scores.txt\"\n",
    "#main_filtered_f = main_f\n",
    "main_filtered_f = \"/opt/data/workingdir/blhill/main_merged_w_akin_spo2.filtered.main.txt\"\n",
    "test_or_case_id_f = \"/opt/data/workingdir/blhill/or_case_id_032018.txt\"\n",
    "\n",
    "# changing the exp_prefix directory below decides which features to load\n",
    "# and where to save output files\n",
    "exp_prefix = \"preop_no_lab_times\"\n",
    "impute_asa = False\n",
    "dir_to_save_files = os.path.join(\"paper/\", exp_prefix)\n",
    "if not os.path.exists(dir_to_save_files):\n",
    "    os.makedirs(dir_to_save_files)\n",
    "\n",
    "# set to False to hide IDs    \n",
    "verbose = False\n",
    "# default size for plots\n",
    "default_size=(8,8)\n",
    "# model type == \"regression\" || \"classification\"\n",
    "model_type = \"classification\"\n",
    "# AKIN stage > AKIN_THRESHOLD considered an AKIN_EVENT\n",
    "AKIN_THRESHOLD=1\n",
    "# data points measured after this time will be discarded\n",
    "TIME_THRESHOLD=-60*24*0\n",
    "# if True, subsample to ensure classes are split 50/50\n",
    "USE_EQUAL_CLASS_FREQ=False\n",
    "CLASS_MULTIPLIER=1.0\n",
    "\n",
    "MIN_ASA_STATUS=1\n",
    "MAX_ASA_STATUS=5\n",
    "MIN_AGE=18\n",
    "MAX_AGE=89\n",
    "\n",
    "## this variable is the column that we will use as the target variable for the model\n",
    "TARGET_VARIABLE = 'INPT_DEATH_YN'\n",
    "#TARGET_VARIABLE = 'AKIN_STAGE1_YN'\n",
    "#TARGET_VARIABLE = 'MAX_AKIN_STAGE'\n",
    "#TARGET_VARIABLE = 'AKIN_EVENT'\n",
    "#TARGET_VARIABLE = 'FLOOR_2_ICU_YN'\n",
    "#TARGET_VARIABLE = \"ASA_STATUS\"\n",
    "\n",
    "\n",
    "df = pd.read_csv(main_filtered_f,sep=\"|\")\n",
    "print(df.shape)\n",
    "#display(df.PATIENT_CLASS.unique())\n",
    "print(len(df.columns))\n",
    "if verbose:\n",
    "    display(df.iloc[0:20, :])\n",
    "\n",
    "def get_sha256_hash(x):\n",
    "    m = hashlib.sha256()\n",
    "    m.update(str(x).encode('utf-8'))\n",
    "    return m.hexdigest().upper()\n",
    "df[\"OR_CASE_ID\"] = df[\"OR_CASE_ID\"].apply(get_sha256_hash)\n",
    "    \n",
    "# make sure binary columns are set to True/False\n",
    "try:\n",
    "    df.NITRIC_OXIDE = df.NITRIC_OXIDE.astype(bool)\n",
    "    df.ART_LINE = df.ART_LINE.astype(bool)\n",
    "    df.PA_CATHETER = df.PA_CATHETER.astype(bool)\n",
    "    df.INPT_DEATH_YN = df.INPT_DEATH_YN.astype(bool)\n",
    "except AttributeError:\n",
    "    pass\n",
    "if verbose:\n",
    "    display(df.describe(include=\"all\"))\n",
    "ids = df[\"ADMSN_ID\"]\n",
    "#display(df[ids.isin(ids[ids.duplicated()])].head(30))\n",
    "# display(df[df[\"INPT_DEATH_YN\"] == True].sort_values(by=[\"ADMSN_ID\", \"HRS_ADMSN_TO_SURGERY\"]).head(30))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.columns.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get methylation sample ID to AtLAs Biobank ID mapping\n",
    "meth_map_f = \"/opt/data/azuremethylationcontainer/2019-9110-1 Halperin Meth Epic/2019-9110 Sample sheet.xlsx\"\n",
    "meth_map = pd.read_excel(meth_map_f, dtype=str)\n",
    "meth_map[\"meth_ID\"] = meth_map[\"Chip Barcode ID\"] + \"_\" + meth_map[\"Stripe ID\"]\n",
    "meth_map[\"External Sample ID\"] = meth_map[\"External Sample ID\"].str.split(\".\", n=1, expand=True).iloc[:,0]\n",
    "meth_map.head()\n",
    "meth_map.describe(include=\"all\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get list of surgeries for testing\n",
    "test_or_case_ids = pd.read_csv(test_or_case_id_f, header=None)\n",
    "print(test_or_case_ids.shape)\n",
    "# get list of patient IDs for testing - we want to remove these from the training set\n",
    "test_pat_ids = df[df.OR_CASE_ID.isin(test_or_case_ids.iloc[:,0])][\"PAT_ID\"]\n",
    "print(test_pat_ids.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filt_df = pd.read_csv(\"/opt/data/workingdir/blhill/main_merged_w_akin_spo2.filtered.filtIDX.txt\", sep=\",\")\n",
    "print(filt_df.ANES_CASE_YN.sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "main_unfiltered_f = \"/opt/data/workingdir/blhill/main_merged_w_akin_spo2.txt\"\n",
    "df_unfiltered = pd.read_csv(main_unfiltered_f, sep=\"|\")\n",
    "df_unfiltered[\"OR_CASE_ID\"] = df_unfiltered[\"OR_CASE_ID\"].apply(get_sha256_hash)\n",
    "df_unfiltered.drop_duplicates(inplace=True)\n",
    "unfiltered_test_pat_ids = df_unfiltered[df_unfiltered.OR_CASE_ID.isin(test_or_case_ids.iloc[:,0])][\"PAT_ID\"]\n",
    "print(unfiltered_test_pat_ids.shape)\n",
    "print(df_unfiltered.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df_unfiltered.shape)\n",
    "df_unfiltered = df_unfiltered[df_unfiltered[\"ANES_CASE_YN\"] == 1]\n",
    "print(df_unfiltered.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(df_unfiltered.shape)\n",
    "df_unfiltered = df_unfiltered[df_unfiltered['LOCATION_GROUP'].isin(['RR OR', 'SM OR'])]\n",
    "#df_unfiltered = df_unfiltered[(df_unfiltered['LOCATION_GROUP'].str.startswith('RR')) | (df_unfiltered['LOCATION_GROUP'].str.startswith('SM'))]\n",
    "print(df_unfiltered.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df_unfiltered.shape)\n",
    "df_unfiltered = df_unfiltered[df_unfiltered['PAT_CLASS'].isin(['INPATIENT', 'SAME DAY ADMIT', 'EMERGENCY', 'OVERNIGHT RECOVERY'])]\n",
    "print(df_unfiltered.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df_unfiltered.shape)\n",
    "#df_unfiltered = df_unfiltered[(df_unfiltered[\"ASA_STATUS\"] <= MAX_ASA_STATUS) & (df_unfiltered[\"ASA_STATUS\"] >= MIN_ASA_STATUS)]\n",
    "df_unfiltered = df_unfiltered[(df_unfiltered[\"ASA_STATUS\"] <= MAX_ASA_STATUS)]\n",
    "print(df_unfiltered.shape)\n",
    "df_unfiltered = df_unfiltered[(df_unfiltered[\"AGE_LT_90\"] <= MAX_AGE) & (df_unfiltered[\"AGE_LT_90\"] > MIN_AGE)]\n",
    "print(df_unfiltered.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print_demographics_info(df_unfiltered[~(df_unfiltered.PAT_ID.isin(test_pat_ids))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_unfiltered.PAT_ID = df_unfiltered.PAT_ID.astype(int)\n",
    "print(df_unfiltered.PAT_ID.dtype)\n",
    "print(df_unfiltered[~(df_unfiltered.PAT_ID.isin(test_pat_ids))].shape)\n",
    "print(df_unfiltered[df_unfiltered.OR_CASE_ID.isin(test_or_case_ids.iloc[:,0])].shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # if a patient dies in an admission, set only the last surgery in the admission to 1\n",
    "# for name, group in df.groupby(\"ADMSN_ID\"):\n",
    "#     if group.shape[0] > 1 and group.iloc[-1,:][\"INPT_DEATH_YN\"] == 1.0:\n",
    "#         sorted_group = group.sort_values(\"ADMSN_SURGERY_NUMBER_W_ANES\")\n",
    "#         print(sorted_group)\n",
    "#         or_case_ids_to_change = sorted_group.iloc[:-1, :][\"OR_CASE_ID\"]\n",
    "#         print(\"=\"*30)\n",
    "#         for ocid in or_case_ids_to_change:\n",
    "#             print(ocid)\n",
    "#             df.loc[df.OR_CASE_ID == ocid, \"INPT_DEATH_YN\"] = 0.\n",
    "#         print(\"=\"*30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print df[\"_merge\"].unique()\n",
    "# print df[\"AGE_LT_90\"].isnull().sum()\n",
    "# print df[df[\"_merge\"] == 'both'].shape\n",
    "# print df.INPT_DEATH_YN.mean()\n",
    "# df2 = pd.read_csv(\"Main_Data_filtered.txt\", sep=\"\\t\")\n",
    "# print df2.shape\n",
    "\n",
    "# print df.shape\n",
    "# print df[df.OR_CASE_ID.isin(df2.OR_CASE_ID)].shape\n",
    "# print df[df.OR_CASE_ID.isin(df2.OR_CASE_ID)].INPT_DEATH_YN.mean()\n",
    "\n",
    "# df2.drop_duplicates(inplace=True)\n",
    "# print df2[\"AGE_LT_89\"].isnull().sum()\n",
    "# # df2.sort_values(by=[\"ADMSN_ID\", \"HRS_ADMSN_TO_SURGERY\"], inplace=True)\n",
    "# # df2.drop_duplicates(subset=\"ADMSN_ID\", keep=\"first\", inplace=True)\n",
    "# # df2 = df2[df2['PAT_CLASS'].isin(['INPATIENT', 'SAME DAY ADMIT', 'EMERGENCY', 'OVERNIGHT RECOVERY'])]\n",
    "# #df2 = df2[df2['LOCATION_GROUP'].isin(['RR OR', 'SM OR', 'SM SC'])]\n",
    "# #df2 = df2[(df2[\"AGE_LT_89\"] <= MAX_AGE) & (df2[\"AGE_LT_89\"] >= MIN_AGE)]\n",
    "# print df2.shape\n",
    "# print df2.INPT_DEATH_YN.mean()\n",
    "\n",
    "# outcomes_df = pd.read_csv(os.path.join(data_dir, \"Outcomes_Data.txt\"), sep='|', header=0)\n",
    "\n",
    "# old_df = pd.read_csv(\"Main_Data_filtered.working.txt\",sep=\"\\t\")\n",
    "# print \"old.df shape:\", old_df.shape\n",
    "# old_df.drop_duplicates(inplace=True)\n",
    "# print old_df[\"AGE_LT_90\"].isnull().sum()\n",
    "# #df.drop(['INPT_DEATH_YN'], axis=1, inplace=True)\n",
    "# # old_df = old_df[old_df[\"ADMSN_SURGERY_NUMBER\"] == 1]\n",
    "# # old_df = old_df[old_df[\"ANES_CASE_YN\"] == 1]\n",
    "# #old_df = old_df[(old_df[\"AGE_LT_90\"] <= MAX_AGE) & (old_df[\"AGE_LT_90\"] >= MIN_AGE)]\n",
    "# # old_df = old_df[old_df['PATIENT_CLASS'].isin(['INPATIENT', 'SAME DAY ADMIT', 'EMERGENCY', 'OVERNIGHT RECOVERY'])]\n",
    "# print old_df.shape\n",
    "# old_df = old_df.merge(outcomes_df, on=\"ADMSN_ID\", how=\"inner\")\n",
    "# print \"shape after merging outcomes:\", old_df.shape\n",
    "# print \"old df mean:\", old_df.INPT_DEATH_YN.mean()\n",
    "\n",
    "# df2 = df2.merge(old_df[[\"OR_CASE_ID\", \"GENDER\", \"HEIGHT_IN\", \"WEIGHT_KG\", \"BMI\", \"ASA_STATUS\", \"ADMSN_SURGERY_NUMBER\"]], on=\"OR_CASE_ID\", how=\"inner\")\n",
    "# print \"merged shape:\", df2.shape\n",
    "# print df2[~df2.OR_CASE_ID.isin(old_df.OR_CASE_ID)].shape\n",
    "# print df2[~df2.OR_CASE_ID.isin(old_df.OR_CASE_ID)].INPT_DEATH_YN.mean()\n",
    "\n",
    "# if verbose:\n",
    "#     display(df2.head())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FEATURES_PATH = os.path.join(dir_to_save_files, 'EHR_MAIN_FEATURES.csv')\n",
    "features_df = pd.read_csv(FEATURES_PATH)\n",
    "\n",
    "features_dict = {name:list(col.dropna()) for name,col in features_df.items()}\n",
    "print(features_dict.keys())\n",
    "\n",
    "final_features = features_dict['final_features']\n",
    "cat_to_drop = features_dict['cat_to_drop']\n",
    "outcome_vars = features_dict['outcome_vars']\n",
    "feat_to_drop = features_dict['feat_to_drop']\n",
    "cat_vars = features_dict['cat_vars']\n",
    "contin_vars = features_dict['contin_vars']\n",
    "bool_outcome_vars = features_dict['bool_outcome_vars']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove HRS_ADMSN_TO_SURGERY from feature list\n",
    "try:\n",
    "    final_features.remove('HRS_ADMSN_TO_SURGERY')\n",
    "    contin_vars.remove('HRS_ADMSN_TO_SURGERY')\n",
    "    cat_vars.remove('HCUP_CODE')\n",
    "    cat_vars.remove('PRE_SURG_LOCATION')\n",
    "except ValueError:\n",
    "    pass\n",
    "\n",
    "cat_to_drop.remove(\"CASE_SRV_NAME\")\n",
    "\n",
    "if 'HRS_ADMSN_TO_SURGERY' not in feat_to_drop:\n",
    "    feat_to_drop = feat_to_drop.append('HRS_ADMSN_TO_SURGERY')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remove duplicate rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove any people that had multiple surgeries\n",
    "print(df.shape)\n",
    "display(df[\"INPT_DEATH_YN\"].mean())\n",
    "# drop any rows that are exact copies of another row\n",
    "df.drop_duplicates(inplace=True)\n",
    "print(df.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#ids = df[\"ADMSN_ID\"]\n",
    "# display(df[ids.isin(ids[ids.duplicated()]) & (df[\"INPT_DEATH_YN\"] == True)].sort_values(by=[\"ADMSN_ID\", \"HRS_ADMSN_TO_SURGERY\"]).head(30))\n",
    "# display(df[ids.isin(ids[ids.duplicated()]) & (df[\"INPT_DEATH_YN\"] == True)].sort_values(by=[\"ADMSN_ID\", \"HRS_ADMSN_TO_SURGERY\"]).groupby(\"ADMSN_ID\").apply(lambda x: x.corr()).head(30))\n",
    "# display(df[\"INPT_DEATH_YN\"].mean())\n",
    "# display(df[df[\"INPT_DEATH_YN\"] == True & df[\"HCUP_CODE\"].isnull()].head(300))\n",
    "\n",
    "#display(df[df[\"INPT_DEATH_YN\"] == True].head(300))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# if multiple surgeries for a single ADMSN_ID, take first surgery"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if multiple surgeries for a single ADMSN_ID, take first surgery\n",
    "# print df.shape\n",
    "# df.sort_values(by=[\"ADMSN_ID\", \"HRS_ADMSN_TO_SURGERY\"], inplace=True)\n",
    "# df.drop_duplicates(subset=\"ADMSN_ID\", keep=\"first\", inplace=True)\n",
    "# display(df[\"INPT_DEATH_YN\"].mean())\n",
    "# print df.shape\n",
    "# display(df.describe(include=\"all\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(df.CASE_SRV_NAME.unique()))\n",
    "print(len(df.CASE_SRV_NAME_GROUP.unique()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create outcomes data frame "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create data frame for outcomes table\n",
    "# print outcomes_f\n",
    "# outcomes_df = pd.read_csv(outcomes_f, sep='|', header=0)\n",
    "# outcomes_df.INPT_DEATH_YN = outcomes_df.INPT_DEATH_YN.astype(bool)\n",
    "# outcomes_df.AKIN_STAGE1_YN = outcomes_df.AKIN_STAGE1_YN.apply(lambda x: True if x == 1 else False)\n",
    "# outcomes_df.AKIN_STAGE2_YN = outcomes_df.AKIN_STAGE2_YN.apply(lambda x: True if x == 1 else False)\n",
    "# outcomes_df.AKIN_STAGE3_YN = outcomes_df.AKIN_STAGE3_YN.apply(lambda x: True if x == 1 else False)\n",
    "# outcomes_df.POSTOP_REINTUBATION_YN = outcomes_df.POSTOP_REINTUBATION_YN.astype(bool)\n",
    "# outcomes_df.ADMSN_TRACH_YN = outcomes_df.ADMSN_TRACH_YN.apply(lambda x: True if x == 1 else False)\n",
    "# outcomes_df.FLOOR_2_ICU_YN = outcomes_df.FLOOR_2_ICU_YN.astype(bool)\n",
    "# create new feature AKIN_EVENT, true if MAX_AKIN_STAGE > AKIN_THRESHOLD else false\n",
    "df['INPT_DEATH_YN'] = df.INPT_DEATH_YN.astype(bool)\n",
    "df['AKIN_EVENT'] = df['MAX_AKIN_STAGE'].apply(lambda val: True if val > AKIN_THRESHOLD else False)\n",
    "# display(outcomes_df.describe(include=\"all\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check frequency of target variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# bool_outcome_vars = ['INPT_DEATH_YN', 'AKIN_STAGE1_YN', 'AKIN_STAGE2_YN', 'AKIN_STAGE3_YN', 'POSTOP_REINTUBATION_YN', \n",
    "#                      'ADMSN_TRACH_YN', 'FLOOR_2_ICU_YN', 'AKIN_EVENT']\n",
    "\n",
    "display(df[bool_outcome_vars].mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Join Main data frame with outcomes using ADMSN_ID as key\n",
    "\n",
    "To add/remove ASA_STATUS from the model, add/remove \"ASA_STATUS\" from the list of columns below\n",
    "To get the number of patients in RR OR/SM OR, add \"LOCATION_GROUP\" to list of columns below "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print df.shape\n",
    "# old_df = pd.read_csv(main_filtered_working_f,sep=\"\\t\")\n",
    "# print old_df.shape\n",
    "# #df = df.merge(old_df[[\"ADMSN_ID\", \"GENDER\"]], on=\"ADMSN_ID\", how=\"inner\")\n",
    "# # removed \n",
    "# df = df.merge(old_df[[\"OR_CASE_ID\", \"GENDER\", \"HEIGHT_IN\", \"WEIGHT_KG\", \"BMI\", \"ASA_STATUS\", \"ADMSN_SURGERY_NUMBER\"]], on=\"OR_CASE_ID\", how=\"inner\")\n",
    "# print df.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Use PRED_ASA_STATUS as a feature\n",
    "\n",
    "To use PRED_ASA_STATUS as a feature in the model, uncomment out the block of code below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if \"imp_asa\" in exp_prefix:\n",
    "    print(\"Adding training data predicted ASA\")\n",
    "    asa_train_df = pd.read_csv(asa_predictions_train_f, sep=\"\\t\")\n",
    "    print(\"asa_train_df shape:\", asa_train_df.shape)\n",
    "    \n",
    "    print(\"Adding testing data predicted ASA\")\n",
    "    asa_test_df = pd.read_csv(asa_predictions_test_f, sep=\"\\t\")\n",
    "    print(\"asa_test_df shape:\", asa_test_df.shape)\n",
    "    \n",
    "    asa_df = asa_train_df.append(asa_test_df)\n",
    "    print(\"asa_df shape:\", asa_df.shape)\n",
    "    \n",
    "    df = df.merge(asa_df[[\"OR_CASE_ID\", \"PRED_ASA_STATUS\"]], on=\"OR_CASE_ID\", how=\"left\")\n",
    "    print(df.shape)\n",
    "    df.drop_duplicates(inplace=True)\n",
    "    print(df.shape)\n",
    "\n",
    "    if verbose:\n",
    "        df.head(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Use Charlson comorbidity as an input feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if exp_prefix == \"charlson\": \n",
    "    print(df.shape)\n",
    "    comorbid_df = pd.read_csv(comorbid_f, sep=\"|\")\n",
    "    print(comorbid_df.shape)\n",
    "    comorbid_df.columns = [\"ADMSN_ID\", \"CHARLSON\"]\n",
    "    comorbid_df = comorbid_df[~(comorbid_df[\"ADMSN_ID\"].isna())]\n",
    "    comorbid_df[\"ADMSN_ID\"] = comorbid_df[\"ADMSN_ID\"].astype(int)\n",
    "    print(\"comorbid_df shape:\", comorbid_df.shape)\n",
    "    comorbid_df.drop_duplicates(inplace=True)\n",
    "    comorbid_df.drop_duplicates(subset=\"ADMSN_ID\", keep=\"first\", inplace=True)\n",
    "    print(comorbid_df.head(20))\n",
    "    print(comorbid_df.dtypes)\n",
    "    print(\"comorbid_df shape after dedup:\", comorbid_df.shape)\n",
    "    df[\"ADMSN_ID\"] = df[\"ADMSN_ID\"].astype(int)\n",
    "    df = df.merge(comorbid_df[[\"ADMSN_ID\", \"CHARLSON\"]], on=\"ADMSN_ID\", how=\"left\")\n",
    "    print(df.shape)\n",
    "    df.drop_duplicates(inplace=True)\n",
    "    print(df.shape)\n",
    "    # df.drop(asa_pred_features_to_drop, axis=1, inplace=True)\n",
    "    # print df.shape\n",
    "    #df.head(20)\n",
    "    if verbose:\n",
    "        df.describe(include=\"all\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Use POSPOM score as an input feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if exp_prefix == \"pospom\": \n",
    "    print(df.shape)\n",
    "    pospom_df = pd.read_csv(pospom_f, sep=\",\", header=0)\n",
    "    print(\"pospom_df shape:\", pospom_df.shape)\n",
    "    pospom_df.drop_duplicates(inplace=True)\n",
    "    # add surgical POSPOM score to diagnosis POSPOM score\n",
    "    #pospom_df[\"POSPOM\"] = pospom_df[\"SURGICAL_POSPOM_COMPONENT\"] + pospom_df[\"POSPOM_SCORE\"]\n",
    "    print(pospom_df[\"SURGICAL_POSPOM_COMPONENT\"].isna().mean())\n",
    "    print(pospom_df[\"POSPOM_SCORE\"].isna().mean())\n",
    "    \n",
    "    pospom_df[\"OR_CASE_ID\"] = pospom_df[\"OR_CASE_ID\"].apply(get_sha256_hash)\n",
    "    \n",
    "    pospom_df[\"POSPOM_SCORE\"].fillna(0, inplace=True)\n",
    "    print(pospom_df[\"POSPOM_SCORE\"].isna().mean())\n",
    "    #pospom_df[\"POSPOM\"] = pospom_df[\"SURGICAL_POSPOM_COMPONENT\"]\n",
    "    pospom_df[\"POSPOM\"] = pospom_df[\"SURGICAL_POSPOM_COMPONENT\"] + pospom_df[\"POSPOM_SCORE\"]\n",
    "    #pospom_df.drop([\"SURGICAL_POSPOM_COMPONENT\", \"POSPOM_SCORE\"], axis=1, inplace=True)\n",
    "    pospom_df.drop([\"POSPOM_SCORE\"], axis=1, inplace=True)\n",
    "    #pospom_df.drop_duplicates(subset=\"ADMSN_ID\", keep=\"first\", inplace=True)\n",
    "    print(\"pospom_df shape after dedup:\", pospom_df.shape)\n",
    "    df = df.merge(pospom_df, on=\"OR_CASE_ID\", how=\"left\")\n",
    "    print(df.shape)\n",
    "    df.drop_duplicates(inplace=True)\n",
    "    print(df.shape)\n",
    "    print(df[\"SURGICAL_POSPOM_COMPONENT\"].value_counts()/df.shape[0])\n",
    "    df[\"SURGICAL_POSPOM_COMPONENT\"].hist(bins=20)\n",
    "    #df[\"POSPOM\"].hist(bins=20)\n",
    "    plt.show()\n",
    "    # df.drop(asa_pred_features_to_drop, axis=1, inplace=True)\n",
    "    # print df.shape\n",
    "    #df.head(20)\n",
    "    if verbose:\n",
    "        df.describe(include=\"all\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# # join main data frame with outcomes data\n",
    "# df.drop(['INPT_DEATH_YN'], axis=1, inplace=True)\n",
    "# df = df.merge(outcomes_df, on=\"ADMSN_ID\", how=\"inner\")\n",
    "# print \"Column names:\", df.columns.values\n",
    "# print \"df shape:\", df.shape\n",
    "\n",
    "# display(df[bool_outcome_vars].mean())\n",
    "# display(df.describe(include=\"all\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remove patients with too many missing values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "missing_threshold = 0.3\n",
    "print(df.shape)\n",
    "df.dropna(thresh=int(missing_threshold*df.shape[1]), inplace=True)\n",
    "print(df.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check number of null values for each feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get count of null values in each column\n",
    "print(df.isnull().sum()/df.shape[0])\n",
    "plt.figure(figsize=(14,5))\n",
    "plt.bar(np.arange(df.shape[1]), df.isnull().sum()/df.shape[0])\n",
    "plt.xticks(range(0, len(df.columns.values)), df.columns.values)\n",
    "plt.xticks(rotation=90)\n",
    "plt.xlabel('Features')\n",
    "plt.ylabel('% NaN values')\n",
    "plt.gca().yaxis.grid(True)\n",
    "plt.title('Fraction of NaN values in Features')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(df.HCUP_CODE.unique()))\n",
    "#print df[(df[\"ASA_STATUS\"] == 2) & (df[\"INPT_DEATH_YN\"] == True)].shape\n",
    "#df.groupby([\"INPT_DEATH_YN\",\"ASA_STATUS\"]).hist(layout=(5,1),column=\"HRS_ADMSN_TO_SURGERY\", bins=range(0,100,5))\n",
    "print(df.shape)\n",
    "#df = df[df['HRS_ADMSN_TO_SURGERY'] <= 5.]\n",
    "print(df.shape)\n",
    "# plt.xlim(2000,4000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filter out surgeries that don't occur in RR or SM operating rooms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df.LOCATION_GROUP.unique())\n",
    "#df = df[df['LOCATION_GROUP'].isin(['RR OR', 'SM OR','SM SC','SM OB OR','RR OB OR'])]\n",
    "print(df.shape)\n",
    "#df = df[df['LOCATION_GROUP'].isin(['RR OR', 'SM OR', 'SM SC'])]\n",
    "df = df[df['LOCATION_GROUP'].isin(['RR OR', 'SM OR'])]\n",
    "print(df.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filter out surgeries that were not INPATIENT, SAME DAY ADMIT, EMERGENCY, or OVERNIGHT RECOVERY"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter out any surgeries that are not inpatient\n",
    "#print df['PATIENT_CLASS_INPATIENT']\n",
    "print(\"Shape before filtering out outpatient surgeries:\", df.shape)\n",
    "#df = df[df['PATIENT_CLASS'].isin(['INPATIENT', 'SAME DAY ADMIT', 'EMERGENCY', 'OVERNIGHT RECOVERY'])]\n",
    "df = df[df['PAT_CLASS'].isin(['INPATIENT', 'SAME DAY ADMIT', 'EMERGENCY', 'OVERNIGHT RECOVERY'])]\n",
    "print(\"Shape after filtering out outpatient surgeries:\", df.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filter based on ASA_STATUS, AGE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    print(\"Shape before filtering out based on ASA_STATUS:\", df.shape)\n",
    "    print(\"ASA_STATUS mean:\", df.ASA_STATUS.mean())\n",
    "    df = df[(df[\"ASA_STATUS\"] <= MAX_ASA_STATUS) & (df[\"ASA_STATUS\"] >= MIN_ASA_STATUS)]\n",
    "    print(\"Shape after filtering out based on ASA_STATUS:\", df.shape)\n",
    "    print(\"ASA_STATUS mean:\", df.ASA_STATUS.mean())\n",
    "except AttributeError:\n",
    "    pass\n",
    "print(\"===================================\")\n",
    "print(\"Mean age:\", df.AGE_LT_90.mean())\n",
    "print(\"STD age:\", df.AGE_LT_90.std())\n",
    "df = df[(df[\"AGE_LT_90\"] <= MAX_AGE) & (df[\"AGE_LT_90\"] >= MIN_AGE)]\n",
    "print(\"Mean age:\", df.AGE_LT_90.mean())\n",
    "print(\"STD age:\", df.AGE_LT_90.std())\n",
    "print(\"Shape after filtering out based on AGE_LT_90:\", df.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check demographic distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(\"Number of surgeries in admissions with in-hospital mortality: {}\".format(df[(df[\"INPT_DEATH_YN\"]==True) & (df[\"ADMSN_SURGERY_NUMBER\"]==2)].shape))\n",
    "\n",
    "# df.sort_values(by=[\"ADMSN_ID\", \"HRS_ADMSN_TO_SURGERY\"], inplace=True)\n",
    "# df.drop_duplicates(subset=\"ADMSN_ID\", keep=\"first\", inplace=True)\n",
    "# display(df[\"INPT_DEATH_YN\"].mean())\n",
    "# print df.shape\n",
    "# display(df.describe(include=\"all\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df.groupby(\"PAT_ID\").agg({'ADMSN_ID': 'count', 'OR_CASE_ID': 'count'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(\"Total number of surgeries:\", df.shape[0])\n",
    "print(\"Total number of patients:\", df[\"PAT_ID\"].unique().shape[0])\n",
    "print(\"Total mortalities: {} ({}%)\".format(df[df[\"INPT_DEATH_YN\"] == 1].shape[0], float(df[df[\"INPT_DEATH_YN\"] == 1].shape[0])/df.shape[0]*100.))\n",
    "\n",
    "# df.drop_duplicates(inplace=True)\n",
    "# df.sort_values(by=[\"ADMSN_ID\", \"HRS_ADMSN_TO_SURGERY\"], inplace=True)\n",
    "# df.drop_duplicates(subset=\"ADMSN_ID\", keep=\"first\", inplace=True)\n",
    "def print_demographics_info(df):\n",
    "    print(\"Number of unique patients:\", df[\"PAT_ID\"].unique().shape[0])\n",
    "    print(\"Number of unique admissions:\", df[\"ADMSN_ID\"].unique().shape[0])\n",
    "    print(\"Number of surgeries:\", df[\"OR_CASE_ID\"].unique().shape[0])\n",
    "    print(\"Average number of surgeries per patient:\", float(df[\"OR_CASE_ID\"].unique().shape[0])/df[\"PAT_ID\"].unique().shape[0])\n",
    "    print(\"Average number of admissions per patient:\", float(df[\"ADMSN_ID\"].unique().shape[0])/df[\"PAT_ID\"].unique().shape[0])\n",
    "    print(\"Average number of surgeries per admission:\", float(df[\"OR_CASE_ID\"].unique().shape[0])/df[\"ADMSN_ID\"].unique().shape[0])\n",
    "    print(\"Number of patients with more than one admission: {} ({}%)\".format(df.PAT_ID.unique()[df.groupby(\"PAT_ID\")[\"ADMSN_ID\"].nunique() > 1].shape[0], df.PAT_ID.unique()[df.groupby(\"PAT_ID\")[\"ADMSN_ID\"].nunique() > 1].shape[0]/float(df[\"PAT_ID\"].unique().shape[0])*100.))\n",
    "    print(\"Number of admissions with more than one surgery: {} ({}%)\".format(df.ADMSN_ID.unique()[df[\"ADMSN_ID\"].value_counts() > 1].shape[0], df.ADMSN_ID.unique()[df[\"ADMSN_ID\"].value_counts() > 1].shape[0]/float(df[\"ADMSN_ID\"].unique().shape[0])*100.))\n",
    "    print(\"Patients with in-hospital mortality: {} ({}%)\".format(df.INPT_DEATH_YN.value_counts()[1], df.INPT_DEATH_YN.value_counts(normalize=\"True\")[1]*100))\n",
    "    #print(\"Patients with kidney failure: {} ({}%)\".format(df.AKIN_EVENT.value_counts()[1], df.AKIN_EVENT.value_counts(normalize=\"True\")[1]*100))\n",
    "    print(\"Mean age:\", df.AGE_LT_90.mean(), \" std:\", df.AGE_LT_90.std())\n",
    "    print(\"Number of female patients: {} ({}%)\".format(df[df[\"GENDER\"] == \"F\"].shape[0], df[df[\"GENDER\"] == \"F\"].shape[0]/float(df.shape[0])*100))\n",
    "\n",
    "    try:\n",
    "        print(\"Number of patients in RR OR: {} ({}%)\".format(df.LOCATION_GROUP.value_counts()[\"RR OR\"], df.LOCATION_GROUP.value_counts(normalize=\"True\")[\"RR OR\"]*100))\n",
    "        print(\"Number of patients in SM OR: {} ({}%)\".format(df.LOCATION_GROUP.value_counts()[\"SM OR\"], df.LOCATION_GROUP.value_counts(normalize=\"True\")[\"SM OR\"]*100))\n",
    "        #print(\"Number of patients in SM SC: {} ({}%)\".format(df.LOCATION_GROUP.value_counts()[\"SM SC\"], df.LOCATION_GROUP.value_counts(normalize=\"True\")[\"SM SC\"]*100))\n",
    "    except AttributeError:\n",
    "        pass\n",
    "    print(\"=\"*40)\n",
    "    try:\n",
    "        print(\"ASA Status:\", df.ASA_STATUS.value_counts())\n",
    "        print(\"ASA Status (%):\", (df.ASA_STATUS.value_counts()/df.shape[0])*100)\n",
    "    except AttributeError:\n",
    "        pass\n",
    "\n",
    "    print(\"=\"*40)\n",
    "    print(\"Number of mortalities: {} ({}%)\".format(df[df[\"INPT_DEATH_YN\"] == 1].shape[0], float(df[df[\"INPT_DEATH_YN\"] == 1].shape[0])/df.shape[0]*100.))\n",
    "    print(\"Mean age of mortalities:\", df[df[\"INPT_DEATH_YN\"] == 1].AGE_LT_90.mean(), \" std:\", df[df[\"INPT_DEATH_YN\"] == 1].AGE_LT_90.std())\n",
    "    print(\"Number of female mortalities: {} ({}%)\".format(df[df[\"GENDER\"] == \"F\"][\"INPT_DEATH_YN\"].sum(), \n",
    "                                                         df[df[\"GENDER\"] == \"F\"][\"INPT_DEATH_YN\"].sum()/float(df[\"INPT_DEATH_YN\"].sum())*100))\n",
    "    print(\"Number of male mortalities: {} ({}%)\".format(df[df[\"GENDER\"] == \"M\"][\"INPT_DEATH_YN\"].sum(), \n",
    "                                                         df[df[\"GENDER\"] == \"M\"][\"INPT_DEATH_YN\"].sum()/float(df[\"INPT_DEATH_YN\"].sum())*100))\n",
    "    try:\n",
    "        print(\"=\"*40)\n",
    "        print(\"Number of mortalities stratified by location\")\n",
    "        print(df.groupby(\"LOCATION_GROUP\")[\"INPT_DEATH_YN\"].sum())                                \n",
    "        print(df.groupby(\"LOCATION_GROUP\")[\"INPT_DEATH_YN\"].sum()/float(df[\"INPT_DEATH_YN\"].sum())*100)\n",
    "    except AttributeError:\n",
    "        pass\n",
    "    print(\"=\"*40)\n",
    "    print(\"Number of mortalities stratified by ASA status\")\n",
    "    print(df.groupby(\"ASA_STATUS\")[\"INPT_DEATH_YN\"].sum())\n",
    "    print(df.groupby(\"ASA_STATUS\")[\"INPT_DEATH_YN\"].sum()/float(df[\"INPT_DEATH_YN\"].sum())*100)\n",
    "\n",
    "print(\"=\"*40)\n",
    "print(\"TRAINING DATA\")\n",
    "print(\"=\"*40)\n",
    "print_demographics_info(df[~(df.PAT_ID.isin(test_pat_ids))])\n",
    "print(\"=\"*40)\n",
    "print(\"TESTING DATA\")\n",
    "print(\"=\"*40)\n",
    "print_demographics_info(df[(df.OR_CASE_ID.isin(test_or_case_ids.iloc[:,0]))])\n",
    "    \n",
    "if verbose:\n",
    "    print(\"=\"*40)\n",
    "    print(\"TRAINING DATA\")\n",
    "    print(\"=\"*40)\n",
    "    print(df[~(df.PAT_ID.isin(test_pat_ids))][\"CASE_SRV_NAME\"].value_counts())\n",
    "    print(df[~(df.PAT_ID.isin(test_pat_ids))][\"CASE_SRV_NAME\"].value_counts()/df[~(df.PAT_ID.isin(test_pat_ids))].shape[0]*100)\n",
    "    print(\"=\"*40)\n",
    "    print(\"=\"*40)\n",
    "    print(\"TESTING DATA\")\n",
    "    print(\"=\"*40)\n",
    "    print(df[(df.OR_CASE_ID.isin(test_or_case_ids.iloc[:,0]))][\"CASE_SRV_NAME\"].value_counts())\n",
    "    print(df[(df.OR_CASE_ID.isin(test_or_case_ids.iloc[:,0]))][\"CASE_SRV_NAME\"].value_counts()/df[(df.OR_CASE_ID.isin(test_or_case_ids.iloc[:,0]))].shape[0]*100)\n",
    "# display(df.describe(include=\"all\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = df.groupby(\"CASE_SRV_NAME\")\n",
    "kv = {}\n",
    "for n, grp in g:\n",
    "    kv[n] = grp[\"INPT_DEATH_YN\"].sum()\n",
    "import operator\n",
    "sorted_x = sorted(kv.items(), key=operator.itemgetter(1), reverse=True)\n",
    "# for i in sorted_x:\n",
    "#     print(\"{}:\\t{}\".format(i[0], i[1]))\n",
    "pd.DataFrame(sorted_x, columns=[\"CASE_SRV_NAME\", \"INPT_DEATH_COUNT\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# df.boxplot(column=\"PRED_ASA_STATUS\", by=\"ASA_STATUS\")\n",
    "# plt.title(\"\")\n",
    "# plt.suptitle(\"\")\n",
    "# plt.ylabel(\"PRED_ASA_STATUS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# look at distribution of HCUP codes in dataset and in patients that die\n",
    "print(\"HCUP Codes:\") \n",
    "hcup_df = pd.DataFrame([df.HCUP_DESC.value_counts(), (df.HCUP_DESC.value_counts()/df.shape[0])*100]).transpose()\n",
    "hcup_df.columns = [\"HCUP_CODE_COUNT\", \"HCUP_CODE_PERCENT\"]\n",
    "hcup_df.HCUP_CODE_COUNT = hcup_df.HCUP_CODE_COUNT.astype('int32')\n",
    "hcup_df.HCUP_CODE_PERCENT = hcup_df.HCUP_CODE_PERCENT.round(1)\n",
    "#display(hcup_df)\n",
    "hcup_df.to_csv(os.path.join(dir_to_save_files, \"HCUP_codes_full_dataset.txt\"), sep=\"\\t\", header=True, index=True)\n",
    "\n",
    "hcup_inptdeath_df = pd.DataFrame([df[df[\"INPT_DEATH_YN\"] == True].HCUP_DESC.value_counts(), (df[df[\"INPT_DEATH_YN\"] == True].HCUP_DESC.value_counts()/df[df[\"INPT_DEATH_YN\"] == True].shape[0])*100]).transpose()\n",
    "hcup_inptdeath_df.columns = [\"HCUP_CODE_COUNT\", \"HCUP_CODE_PERCENT\"]\n",
    "hcup_inptdeath_df.HCUP_CODE_COUNT = hcup_inptdeath_df.HCUP_CODE_COUNT.astype('int32')\n",
    "hcup_inptdeath_df.HCUP_CODE_PERCENT = hcup_inptdeath_df.HCUP_CODE_PERCENT.round(1)\n",
    "#display(hcup_inptdeath_df)\n",
    "hcup_inptdeath_df.to_csv(os.path.join(dir_to_save_files, \"HCUP_codes_inpt_death.txt\"), sep=\"\\t\", header=True, index=True)\n",
    "\n",
    "df_train = df[~(df.PAT_ID.isin(test_pat_ids))]\n",
    "hcup_df_train = pd.DataFrame([df_train.HCUP_DESC.value_counts(), (df_train.HCUP_DESC.value_counts()/df_train.shape[0])*100]).transpose()\n",
    "hcup_df_train.columns = [\"HCUP_CODE_COUNT\", \"HCUP_CODE_PERCENT\"]\n",
    "hcup_df_train.HCUP_CODE_COUNT = hcup_df_train.HCUP_CODE_COUNT.astype('int32')\n",
    "hcup_df_train.HCUP_CODE_PERCENT = hcup_df_train.HCUP_CODE_PERCENT.round(1)\n",
    "hcup_df_train[\"HCUP_CODE_COUNT_TRAIN\"] = hcup_df_train.apply(lambda row: \"{} ({})\".format(row.HCUP_CODE_COUNT, row.HCUP_CODE_PERCENT), axis=1)\n",
    "#display(hcup_df_train)\n",
    "hcup_df_train.to_csv(os.path.join(dir_to_save_files, \"HCUP_codes_train_data.txt\"), sep=\"\\t\", header=True, index=True)\n",
    "\n",
    "df_test = df[df.OR_CASE_ID.isin(test_or_case_ids.iloc[:,0])]\n",
    "hcup_df_test = pd.DataFrame([df_test.HCUP_DESC.value_counts(), (df_test.HCUP_DESC.value_counts()/df_test.shape[0])*100]).transpose()\n",
    "hcup_df_test.columns = [\"HCUP_CODE_COUNT\", \"HCUP_CODE_PERCENT\"]\n",
    "hcup_df_test.HCUP_CODE_COUNT = hcup_df_test.HCUP_CODE_COUNT.astype('int32')\n",
    "hcup_df_test.HCUP_CODE_PERCENT = hcup_df_test.HCUP_CODE_PERCENT.round(1)\n",
    "hcup_df_test[\"HCUP_CODE_COUNT_TEST\"] = hcup_df_test.apply(lambda row: \"{} ({})\".format(row.HCUP_CODE_COUNT, row.HCUP_CODE_PERCENT), axis=1)\n",
    "#display(hcup_df_test)\n",
    "hcup_df_test.to_csv(os.path.join(dir_to_save_files, \"HCUP_codes_test_data.txt\"), sep=\"\\t\", header=True, index=True)\n",
    "\n",
    "hcup_df_all = hcup_df_train[[\"HCUP_CODE_COUNT_TRAIN\"]].merge(hcup_df_test[[\"HCUP_CODE_COUNT_TEST\"]], how=\"left\", left_index=True, right_index=True)\n",
    "hcup_df_all.to_csv(os.path.join(dir_to_save_files, \"HCUP_codes_combined_data.txt\"), sep=\"\\t\", header=True, index=True)\n",
    "hcup_df_all"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if verbose:\n",
    "    df[df.HCUP_CODE == 108]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remove outlier values\n",
    "Outliers are values greater than 3 standard deviations from the mean. We only remove outliers for continuous variables (not categorical). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "string_cols = ['PRE_SURG_LOCATION', 'CASE_SRV_NAME_GROUP', 'CASE_SRV_NAME', \n",
    "               'PRIMARY_CPT', 'GENDER', 'HCUP_DESC', 'CPT_DESC', 'PAT_CLASS', \n",
    "               'OR_CASE_ID', 'ADMSN_ID']\n",
    "dff = df.drop(['PRE_SURG_LOCATION', 'CASE_SRV_NAME_GROUP', 'CASE_SRV_NAME', 'PRIMARY_CPT',\n",
    "                                  'GENDER', 'HCUP_DESC', 'CPT_DESC', 'PAT_CLASS', 'OR_CASE_ID', 'ADMSN_ID'], axis=1)\n",
    "dff = df.select_dtypes(include=['float64'])\n",
    "df_string_cols = df[df.columns.difference(dff.columns.values)]\n",
    "print(df_string_cols.columns.values)\n",
    "#display(dff.describe())\n",
    "print (np.abs(st.zscore(dff, axis=1)) > 3)\n",
    "#print dff.sub(dff.mean()).div(dff.std()).abs().lt(3)\n",
    "df_no_outliers = dff[dff.sub(dff.mean()).div(dff.std()).abs().lt(4)]\n",
    "df_no_outliers[df_string_cols.columns.values] = df_string_cols\n",
    "if verbose:\n",
    "    display(df_no_outliers.describe(include=\"all\"))\n",
    "df = df_no_outliers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Remove variables related to lab times (i.e. *.HRS_2_SURGERY)\n",
    "If you want to include these variables in the model, comment this block of code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # remove variables that have to do with time\n",
    "# cols_to_keep_no_hrs2surgery = [c for c in df.columns if not c.endswith(\"HRS_2_SURGERY\")]\n",
    "# print cols_to_keep_no_hrs2surgery\n",
    "# print len(cols_to_keep_no_hrs2surgery)\n",
    "# df=df[cols_to_keep_no_hrs2surgery]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Remove unnecessary features \n",
    "These features are either not needed (ex. descriptions of surgery codes) or not wanted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save this for checking predictions over time\n",
    "admsn_surgery_number = df[\"ADMSN_SURGERY_NUMBER\"]\n",
    "print(admsn_surgery_number.shape)\n",
    "or_case_id_number = df[\"OR_CASE_ID\"]\n",
    "admsn_ids = df['ADMSN_ID']\n",
    "pat_ids = df['PAT_ID']\n",
    "asa_status = df[\"ASA_STATUS\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cat_to_drop = ['OR_CASE_ID', 'HCUP_DESC', 'CPT_DESC', 'ECHO_CASE',\n",
    "#                    'CASE_SRV_NAME', 'CASE_SRV_NAME_GROUP', 'PRIMARY_CPT']\n",
    "df = df[final_features + [TARGET_VARIABLE]]\n",
    "#print(df.drop(outcome_vars + cat_to_drop, axis=1).shape[1])\n",
    "pd.DataFrame(df[final_features].isnull().sum()/df.shape[0]).apply(lambda x: np.round(x, 3)).to_csv(os.path.join(dir_to_save_files, \"Feature_List_with_null.txt\"), sep=\"|\", header=False, index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df.columns.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[TARGET_VARIABLE].mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot feature correlation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "sns.set(style=\"white\")\n",
    "\n",
    "# Generate a mask for the upper triangle\n",
    "mask = np.zeros_like(df.corr(), dtype=np.bool)\n",
    "mask[np.triu_indices_from(mask)] = True\n",
    "\n",
    "# Set up the matplotlib figure\n",
    "f, ax = plt.subplots(figsize=(16, 13))\n",
    "\n",
    "# Generate a custom diverging colormap\n",
    "cmap = sns.diverging_palette(220, 10, as_cmap=True)\n",
    "\n",
    "# Draw the heatmap with the mask and correct aspect ratio\n",
    "sns.heatmap(df.corr(), mask=mask, cmap=cmap, vmax=1.0, center=0,\n",
    "            square=True, linewidths=.5, cbar_kws={\"shrink\": .5})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## One-hot encode categorical variables "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# one-hot encode values in categorical columns (creates 1 column per possible state)\n",
    "## removed: 'NITRIC_OXIDE', 'ART_LINE', 'CVC', 'PA_CATHETER',\n",
    "#categorical_string_cols = ['PAT_CLASS', 'HCUP_CODE', 'GENDER', 'CASE_SRV_NAME', 'PRE_SURG_LOCATION']\n",
    "#categorical_string_cols = []\n",
    "# for i in categorical_string_cols:\n",
    "#     df[i] = pd.Categorical(df[i]).codes\n",
    "# removed "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#dummy_cols = ['PAT_CLASS', 'HCUP_CODE', 'GENDER', 'CASE_SRV_NAME', 'PRE_SURG_LOCATION']\n",
    "#dummy_cols = ['GENDER']\n",
    "for var in cat_vars:\n",
    "    try:\n",
    "        # drop_first uses k-1 dummies out of k categories\n",
    "        print(var)\n",
    "        if var not in df.columns.values:\n",
    "            print(\"Missing \", var)\n",
    "        \n",
    "        df = pd.get_dummies(df, columns=[var], drop_first=True)\n",
    "        pass\n",
    "    except ValueError:\n",
    "        pass\n",
    "    except KeyError:\n",
    "        pass\n",
    "# remove categorical variables (string values)\n",
    "for var in cat_vars:\n",
    "    try:\n",
    "        df.drop(var, axis=1, inplace=True)\n",
    "        pass\n",
    "    except ValueError:\n",
    "        print(var, 'already dropped')\n",
    "    except KeyError: \n",
    "        print(var, 'already dropped')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remove features we don't want to include"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for cat in cat_to_drop:\n",
    "    try:\n",
    "        df.drop(cat, axis=1, inplace=True)\n",
    "    except ValueError:\n",
    "        print(cat, 'already dropped')\n",
    "    except KeyError:\n",
    "        print(cat, 'already dropped')\n",
    "\n",
    "#print df.columns.values\n",
    "for col in sorted(df.columns.values):\n",
    "    print(col, \"\\t\\t\", df[col].dtype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remove target variables from data frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(\"Column names:\", df.columns.values)\n",
    "\n",
    "try:\n",
    "    print(df.shape)\n",
    "    y = np.ravel(df[TARGET_VARIABLE])\n",
    "    print(y.mean())\n",
    "    input_death_yn = df['INPT_DEATH_YN']\n",
    "    df.drop(TARGET_VARIABLE, axis=1, inplace=True, errors='ignore')\n",
    "    df.drop(outcome_vars, axis=1, inplace=True, errors='ignore')\n",
    "except KeyError:\n",
    "    print(TARGET_VARIABLE, \"already dropped\")\n",
    "# http://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.Imputer.html\n",
    "# default stragegy: mean\n",
    "# if len(feature_whitelist) > 0:\n",
    "#     feature_whitelist = [c for c in feature_whitelist if not c.endswith(\".HRS_2_SURGERY\")]\n",
    "#     df = df[feature_whitelist]\n",
    "print(df.isnull().sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write final feature list to file so we have order of features\n",
    "pd.DataFrame(df.columns.values).to_csv(os.path.join(dir_to_save_files, exp_prefix + \"_final_feature_list.txt\"), header=False, index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Impute missing data\n",
    "\n",
    "Note that if running with only a single feature (i.e. only ASA status) that you will have to change \n",
    "df.iloc[:,0:] to df.iloc[:,1:] since the first column will be the ID and we don't want that"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This uses mean-value imputation\n",
    "# imputer = Imputer()\n",
    "# transformed_X = imputer.fit_transform(df.iloc[:,1:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print df.shape\n",
    "# [transformed_X, mu_hat, U, s, Vt] = emsvd(df.iloc[:,1:], k=8, tol=1E-3)\n",
    "#print mu_hat\n",
    "#df.drop([\"ASA_STATUS\"], inplace=True)\n",
    "\n",
    "# try: \n",
    "#     imputer = Imputer(strategy=\"median\")\n",
    "#     #transformed_X = imputer.fit_transform(df.iloc[:,1:])\n",
    "#     transformed_X = SoftImpute(max_iters=200).complete(df.iloc[:,0:])\n",
    "#     #transformed_X = MICE().complete(df.iloc[:,1:])\n",
    "#     #transformed_X = df.iloc[:,1:]\n",
    "# except ValueError as e:\n",
    "#     print \"ERROR:\", e\n",
    "#     transformed_X = np.array(df.iloc[:,0:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pd.DataFrame(data=transformed_X, columns=df.columns.values[0:]).describe(include=\"all\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# flatten y into 1D array\n",
    "print(len(y))\n",
    "print(\"Mean of target vector:\", y.mean())\n",
    "#tdf = pd.DataFrame(data=transformed_X)\n",
    "#tdf.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(y.sum())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(or_case_id_number.shape)\n",
    "print(y.shape)\n",
    "# Remove patients (and their surgeries) that occured after March 2018\n",
    "X_train = df[~(pat_ids.isin(test_pat_ids))]\n",
    "y_train = y[~(pat_ids.isin(test_pat_ids))]\n",
    "# X_train = df[~(or_case_id_number.isin(test_or_case_ids.iloc[:,0]))]\n",
    "# y_train = y[~(or_case_id_number.isin(test_or_case_ids.iloc[:,0]))]\n",
    "\n",
    "print(\"X_train shape:\", X_train.shape)\n",
    "print(\"y_train shape:\", y_train.shape)\n",
    "# use surgeries after March 2018 for testing\n",
    "X_test = df[or_case_id_number.isin(test_or_case_ids.iloc[:,0])]\n",
    "y_test = y[or_case_id_number.isin(test_or_case_ids.iloc[:,0])]\n",
    "print(\"X_test shape:\", X_test.shape)\n",
    "print(\"y_test shape:\", y_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save IDs for post-processing\n",
    "X_train_ids = pd.DataFrame()\n",
    "X_train_ids[\"OR_CASE_ID\"] = or_case_id_number[~(pat_ids.isin(test_pat_ids))]\n",
    "X_train_ids[\"ADMSN_ID\"] = admsn_ids[~(pat_ids.isin(test_pat_ids))]\n",
    "X_train_ids[\"PAT_ID\"] = pat_ids[~(pat_ids.isin(test_pat_ids))]\n",
    "X_train_ids[\"INPT_DEATH_YN\"] = y_train\n",
    "\n",
    "X_test_ids = pd.DataFrame()\n",
    "X_test_ids[\"OR_CASE_ID\"] = or_case_id_number[or_case_id_number.isin(test_or_case_ids.iloc[:,0])]\n",
    "X_test_ids[\"ADMSN_ID\"] = admsn_ids[or_case_id_number.isin(test_or_case_ids.iloc[:,0])]\n",
    "X_test_ids[\"PAT_ID\"] = pat_ids[or_case_id_number.isin(test_or_case_ids.iloc[:,0])]\n",
    "X_test_ids[\"INPT_DEATH_YN\"] = y_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Split data into training, testing sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "group_by_admsn = admsn_ids.astype('category').cat.codes\n",
    "print(\"number of operations:\", len(group_by_admsn))\n",
    "print(\"number of unique admission groups:\", len(np.unique(group_by_admsn)))\n",
    "group_kfold = GroupKFold(n_splits=5)\n",
    "\n",
    "# # take the first fold of the k-fold split\n",
    "# train_index, test_index = group_kfold.split(np.array(df.iloc[:,0:]), y, group_by_admsn).next()\n",
    "# #train_index, test_index = group_kfold.split(transformed_X, y, group_by_admsn).next()\n",
    "# X_train, X_test = df.iloc[train_index], df.iloc[test_index]\n",
    "# #X_train, X_test = transformed_X[train_index], transformed_X[test_index]\n",
    "# y_train, y_test = y[train_index], y[test_index]\n",
    "\n",
    "# X_train, X_test, y_train, y_test = train_test_split(transformed_X, y, test_size=0.2, random_state=43)\n",
    "#X_train, X_val, y_train, y_val = train_test_split(X_train, y_train, test_size=0.33, random_state=42)\n",
    "print(\"num X_train before SMOTE:\", len(X_train))\n",
    "print(\"num X_test before SMOTE: \", len(X_test))\n",
    "print(\"Mean of train target vector:\", y_train.mean())\n",
    "print(\"Mean of test target vector:\", y_test.mean())\n",
    "\n",
    "# print \"NUM null before imputation:\", np.count_nonzero(~np.isnan(X_train))\n",
    "# X_train_fit = predictive_imputer.PredictiveImputer(max_iter=10, initial_strategy='median',\n",
    "#                                                    f_model=\"RandomForest\").fit(X_train, y_train)\n",
    "# X_train = X_train_fit.transform(X_train)\n",
    "# print \"NUM null after imputation:\", np.count_nonzero(~np.isnan(X_train))\n",
    "# X_test = X_train_fit.transform(X_test)\n",
    "\n",
    "#X_train, y_train = SMOTE(kind='borderline1', ratio={True:20000},k_neighbors=13).fit_sample(X_train, y_train)\n",
    "# use the SMOTE call below\n",
    "#X_train, y_train = SMOTE(kind='borderline1',k_neighbors=3).fit_sample(X_train, y_train)\n",
    "#X_train, y_train = NearMiss(random_state=0, version=2).fit_sample(X_train, y_train)\n",
    "\n",
    "#X_resampled, y_resampled = SMOTE().fit_sample(transformed_X, y)\n",
    "#X_train, X_test, y_train, y_test = train_test_split(X_resampled, y_resampled, test_size=0.3, random_state=42)\n",
    "print(\"num X_train:\", len(X_train))\n",
    "print(\"num X_test: \", len(X_test))\n",
    "#print \"num X_val:\", len(X_val)\n",
    "print(\"Mean of train target vector:\", y_train.mean())\n",
    "print(\"Mean of test target vector:\", y_test.mean())\n",
    "# plt.hist(y_train)\n",
    "# plt.show()\n",
    "#print \"Mean of val target vector:\", y_val.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if verbose:\n",
    "    display(pd.DataFrame(data=X_train, columns=df.columns.values[0:]).describe(include=\"all\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if verbose:\n",
    "    pd.DataFrame(data=X_train, columns=df.columns.values[0:]).head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# try:\n",
    "#     plt.boxplot(X_train[:,0])\n",
    "#     print transformed_X[np.where(transformed_X[:,df.columns.get_loc(\"AGE_LT_89\")] < 0),0].mean()\n",
    "#     print np.max(X_train[:,df.columns.get_loc(\"AGE_LT_89\")])\n",
    "#     print np.max(X_test[:,df.columns.get_loc(\"AGE_LT_89\")])\n",
    "#     X_train\n",
    "#     print X_train.shape\n",
    "#     print X_test.shape\n",
    "\n",
    "#     print df.columns.get_loc(\"AGE_LT_89\")\n",
    "#     print np.mean(transformed_X[:,df.columns.get_loc(\"AGE_LT_89\")])\n",
    "#     print np.mean(X_train[:,df.columns.get_loc(\"AGE_LT_89\")])\n",
    "#     print np.mean(X_test[:,df.columns.get_loc(\"AGE_LT_89\")])\n",
    "#     #X_train[df.columns.get_loc(\"AGE_LT_89\")-1].std()\n",
    "# except KeyError:\n",
    "#     pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Standardize training, testing data\n",
    "\n",
    "Here we standardize the training set, and use the mean/std on the test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class StandardizeWithNaN(TransformerMixin, BaseEstimator):\n",
    "    '''This estimator is for standardizing a dataset that has missing data'''\n",
    "    def __init__(self):\n",
    "        self.X_mean = []\n",
    "        self.X_std = []\n",
    "        pass\n",
    "\n",
    "    def fit(self, X, y=None):\n",
    "        # get mean and standard deviation of columns\n",
    "        self.X_mean = np.nanmean(X, axis=0)\n",
    "        self.X_std  = np.nanstd(X, axis=0) \n",
    "        return self\n",
    "\n",
    "    def transform(self, X):\n",
    "        # subtract mean and divide by standard deviation\n",
    "        return (X - self.X_mean)/self.X_std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#scaler = StandardScaler()  \n",
    "scaler = StandardizeWithNaN()\n",
    "# Don't cheat - fit only on training data\n",
    "scaler.fit(X_train)\n",
    "# save scaler object for testing on future data\n",
    "pickle.dump(scaler, open(os.path.join(dir_to_save_files, \"StandardizeWithNaN.pkl\"), \"wb\"))\n",
    "# apply transformation to training data\n",
    "X_train = scaler.transform(X_train) \n",
    "# apply same transformation to test data\n",
    "X_test = scaler.transform(X_test) \n",
    "#X_val = scaler.transform(X_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SoftImputeEstimator(TransformerMixin, BaseEstimator):\n",
    "    '''This estimator is for wrapping the SoftImpute algorithm'''\n",
    "    def __init__(self, max_iters=200, verbose=True):\n",
    "        self.max_iters = max_iters\n",
    "        self.verbose = verbose\n",
    "        self.fit_count = 0\n",
    "        self.transform_count = 0\n",
    "        pass\n",
    "\n",
    "    def fit(self, X, y=None):\n",
    "        self.fit_count += 1\n",
    "        print(\"SoftImputeEstimator fit count: {}\".format(self.fit_count))\n",
    "        return self\n",
    "\n",
    "    def transform(self, X):\n",
    "        self.transform_count += 1\n",
    "        print(\"SoftImputeEstimator transform count: {}\".format(self.transform_count))\n",
    "        try:\n",
    "            # subtract mean and divide by standard deviation\n",
    "            return SoftImpute(max_iters=self.max_iters, verbose=self.verbose).complete(X.replace(np.inf, np.nan))\n",
    "        # ValueError raised if no values need to be imputed\n",
    "        except ValueError:\n",
    "            return np.array(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train.dropna(axis=1, how='all', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# X_train.drop([\"HCUP_CODE\", \"PRE_SURG_LOCATION\"], axis=1, inplace=True)\n",
    "# X_test.drop([\"HCUP_CODE\", \"PRE_SURG_LOCATION\"], axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train.columns.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"imputing X_train\")\n",
    "si = Imputer()\n",
    "si.fit(X_train)\n",
    "pickle.dump(si, open(os.path.join(dir_to_save_files, \"MeanImputer.pkl\"), \"wb\"))\n",
    "\n",
    "# si = SoftImputeEstimator()\n",
    "if False:\n",
    "    coded_cat_vars = []\n",
    "    for c in cat_vars:\n",
    "        coded_cat_vars = coded_cat_vars + list(X_train.columns[X_train.columns.str.startswith(c)].values)\n",
    "    print(coded_cat_vars)\n",
    "    cat_var_indices = [X_train.columns.get_loc(c) for c in coded_cat_vars]\n",
    "    print(cat_var_indices)\n",
    "\n",
    "    si = MissForest(n_estimators=500, max_depth=11, min_samples_leaf=5, verbose=0)\n",
    "    si.fit(X_train, cat_vars=cat_var_indices)\n",
    "\n",
    "X_train = si.transform(X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"imputing X_test\")\n",
    "#print(np.isnan(X_test).any())\n",
    "X_test = si.transform(X_test.replace(np.inf, np.nan))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if we want to create imputed ASA status, train on training data and \n",
    "# generate predicted ASA values for train and test set \n",
    "xgb_r = XGBRegressor(max_depth=7, n_estimators=2000, n_jobs=-1)\n",
    "\n",
    "if impute_asa:\n",
    "    assert (\"asa\" not in exp_prefix), \"STOP: cannot impute ASA status when it is included as a predictor\"\n",
    "    # get ASA status for train/test cases\n",
    "    asa_train = asa_status[~(pat_ids.isin(test_pat_ids))]\n",
    "    asa_test = asa_status[or_case_id_number.isin(test_or_case_ids.iloc[:,0])]\n",
    "    \n",
    "    #rfr_pipe = Pipeline([('scaler', scaler), ('impute', imp), ('smt', smt), ('rfr', rfr)])\n",
    "    #asa_predictions = cross_val_predict(rfr_pipe, df.iloc[:,0:], y, cv=group_kfold.split(df.iloc[:,0:], y, group_by_admsn), n_jobs=1)\n",
    "    \n",
    "    t0 = time.time()\n",
    "    xgb_r.fit(X_train, asa_train)\n",
    "    t1 = time.time()\n",
    "    print(\"Took {} seconds to train.\".format(t1-t0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if impute_asa:\n",
    "    from sklearn.metrics import mean_squared_error, r2_score\n",
    "    imp_asa_train = xgb_r.predict(X_train)\n",
    "    imp_asa_test = xgb_r.predict(X_test)\n",
    "\n",
    "    print(\"---Train ASA Imputation---\")\n",
    "    print(\"Train MSE:\", mean_squared_error(asa_train, imp_asa_train))\n",
    "    print(\"Train r^2:\", r2_score(asa_train, imp_asa_train))\n",
    "    print('---Test ASA Imputation---')\n",
    "    print(\"Test MSE:\", mean_squared_error(asa_test, imp_asa_test))\n",
    "    print(\"Train r^2\", r2_score(asa_test, imp_asa_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write predictions to file\n",
    "if impute_asa:\n",
    "    X_train_asa = X_train_ids\n",
    "    X_train_asa[\"ASA_STATUS\"] = asa_train\n",
    "    X_train_asa[\"PRED_ASA_STATUS\"] = imp_asa_train\n",
    "    X_train_asa.to_csv(os.path.join(dir_to_save_files, \"imputed_asa_train.txt\"), sep=\"\\t\", header=True, index=False)\n",
    "    \n",
    "    X_test_asa = X_test_ids\n",
    "    X_test_asa[\"ASA_STATUS\"] = asa_test\n",
    "    X_test_asa[\"PRED_ASA_STATUS\"] = imp_asa_test\n",
    "    X_test_asa.to_csv(os.path.join(dir_to_save_files, \"imputed_asa_test.txt\"), sep=\"\\t\", header=True, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, y_train = SMOTE(kind='borderline1',k_neighbors=3).fit_sample(X_train, y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train model(s)\n",
    "The code is set up to train/test/evaluate any model added to the *models* list. \n",
    "\n",
    "To add a model, simply append it to the list, and append a corresponding name \n",
    "for the model to the *model_names* list as this will be used for plotting. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neural_network import MLPClassifier\n",
    "from sklearn import svm\n",
    "\n",
    "num_threads = -1\n",
    "\n",
    "models = []\n",
    "models_cv = []\n",
    "model_names = []\n",
    "cachedir = mkdtemp()\n",
    "#weight_vec = np.multiply((float(X_train.shape[0]) / (2*np.bincount(y))), [1, 1.2])\n",
    "#print weight_vec\n",
    "\n",
    "#scaler = StandardScaler() \n",
    "scaler = StandardizeWithNaN()\n",
    "imp = SoftImputeEstimator(verbose=True)\n",
    "smt = SMOTE(kind='borderline1', k_neighbors=3)\n",
    "#smt = NearMiss(random_state=0, version=1)\n",
    "pipeline_stages = [('scaler', scaler), ('impute', imp), ('smt', smt)]\n",
    "#pipeline_stages = [('scaler', scaler), ('smt', smt)]\n",
    "\n",
    "## Logistic Regression\n",
    "#models.append(LogisticRegression(class_weight=\"balanced\"))\n",
    "if model_type == \"classification\":\n",
    "    #models.append(LogisticRegression(class_weight={0:weight_vec[0], 1:weight_vec[1]} ,n_jobs=-1))\n",
    "    lrg = LogisticRegression(class_weight=\"balanced\", solver=\"saga\", n_jobs=-1)\n",
    "    models.append(lrg)\n",
    "    models_cv.append(Pipeline(pipeline_stages + [('lrg', lrg)]))\n",
    "    model_names.append(\"Log. Regression\")\n",
    "## ElasticNet\n",
    "#models.append(ElasticNet(l1_ratio=0.0001, fit_intercept=True))\n",
    "\n",
    "elasticnet = ElasticNetCV(max_iter=3000, fit_intercept=True, n_jobs=num_threads, cv=5)\n",
    "elasticnet.fit(X_train, y_train)\n",
    "en_alpha = elasticnet.alpha_\n",
    "en_l1_ratio = elasticnet.l1_ratio_\n",
    "print(\"ElasticNet alpha:\", en_alpha)\n",
    "print(\"ElasticNet l1_ratio:\", en_l1_ratio)\n",
    "# model_names.append(\"ElasticNetCV\")\n",
    "\n",
    "if model_type == \"classification\":\n",
    "#     eln = SGDClassifier(loss=\"log\", penalty=\"elasticnet\", alpha=en_alpha, class_weight=\"balanced\",\n",
    "#                                 l1_ratio=en_l1_ratio, fit_intercept=True, max_iter=3000, n_jobs=num_threads)\n",
    "    eln = LogisticRegression(penalty=\"elasticnet\", \n",
    "                             class_weight=\"balanced\", \n",
    "                             l1_ratio=en_l1_ratio,\n",
    "                             solver=\"saga\",\n",
    "                             n_jobs=-1)\n",
    "    models.append(eln)\n",
    "    models_cv.append(Pipeline(pipeline_stages + [('eln', eln)]))\n",
    "    model_names.append(\"ElasticNet\")\n",
    "\n",
    "# if model_type == \"regression\":\n",
    "#     models.append(ElasticNet(alpha=en_alpha, l1_ratio=en_l1_ratio, fit_intercept=True))\n",
    "\n",
    "## Multilayer Perceptron (MLP)\n",
    "##models.append(MLPClassifier(solver='lbfgs', alpha=1e-5, hidden_layer_sizes=(32, 4), early_stopping=True))\n",
    "#models.append(MLPClassifier(solver='lbfgs', alpha=1e-5, hidden_layer_sizes=(4, 6), early_stopping=True))\n",
    "#model_names.append(\"MLP\")\n",
    "\n",
    "## Random Forest \n",
    "if model_type == \"classification\":\n",
    "    rfc = RandomForestClassifier(class_weight=\"balanced\", \n",
    "                                 oob_score=True, \n",
    "                                 n_estimators=2000, \n",
    "                                 n_jobs=num_threads, \n",
    "                                 max_features=\"log2\",\n",
    "                                 max_depth=7, \n",
    "                                 min_samples_leaf=3)\n",
    "#    rfc = RandomForestClassifier(class_weight={1:1.0, 2:1.0, 3:1.1, 4:1.3, 5:1.5}, oob_score=True, n_estimators=2000, n_jobs=-1, max_features=\"log2\")\n",
    "    models.append(rfc)\n",
    "    models_cv.append(Pipeline(pipeline_stages + [('rfc', rfc)]))\n",
    "    model_names.append(\"Random Forest\")\n",
    "    \n",
    "    xgb = XGBClassifier(max_depth=5, objective=\"binary:logistic\", n_estimators=2000, n_jobs=num_threads, silent=True)\n",
    "    models.append(xgb)\n",
    "    models_cv.append(Pipeline(pipeline_stages + [('xgb', xgb)]))\n",
    "    model_names.append(\"XGBClassifier\")\n",
    "    \n",
    "if model_type == \"regression\":\n",
    "    models.append(RandomForestRegressor(max_depth=5, oob_score=True, n_estimators=2000, n_jobs=num_threads, max_features=\"log2\"))\n",
    "    model_names.append(\"RanfomForestRegressor\")\n",
    "    \n",
    "    xgb = XGBRegressor(max_depth=7, n_estimators=2000, n_jobs=num_threads)\n",
    "    models.append(xgb)\n",
    "    model_names.append(\"XGBRegressor\")\n",
    "# models.append(GradientBoostingClassifier(n_estimators=1000))\n",
    "# model_names.append(\"GradientBoostedClass\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sklearn.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fit the models\n",
    "for mod in models:\n",
    "    t0 = time.time()\n",
    "    mod.fit(X_train, y_train) \n",
    "    t1 = time.time()\n",
    "    print(\"Took {} seconds to train.\".format(t1-t0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dump models to file\n",
    "for i in range(len(models)):\n",
    "    pickled_model_name = \"{}_train_sk{}.pkl\".format(model_names[i], sklearn.__version__)\n",
    "    pickle.dump(models[i], open(os.path.join(dir_to_save_files, pickled_model_name), 'wb'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examine feature relevance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import xgboost\n",
    "# shap.initjs()\n",
    "# #print df.columns.values\n",
    "# #xdf = pd.DataFrame(data=X_train, columns=df.columns.values[1:])\n",
    "# #bst = xgboost.train({\"learning_rate\": 0.01}, xgboost.DMatrix(xdf, label=y_train), 100)\n",
    "# #shap_values = bst.predict(xgboost.DMatrix(xdf), pred_contribs=True)\n",
    "# #shap.force_plot(shap_values[0], xdf[0,:])\n",
    "# print(type(models[2]))\n",
    "# shap_values = shap.TreeExplainer(models[2]).shap_values(X_train[1:500, :])\n",
    "# print shap_values.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print len(shap_values[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# shap.summary_plot(shap_values, X_train[1:500, :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_importance_df = pd.DataFrame()\n",
    "feature_importance_df[\"Features\"] = df.columns.values\n",
    "print(feature_importance_df.shape)\n",
    "\n",
    "for i in range(len(models)):\n",
    "    try:\n",
    "        \n",
    "        feature_importance_df[model_names[i]] = models[i].feature_importances_\n",
    "    except AttributeError:\n",
    "        #print(models[i].coef_[0].shape)\n",
    "        feature_importance_df[model_names[i]] = models[i].coef_[0]\n",
    "feature_importance_df = feature_importance_df[[\"Features\", \"Log. Regression\", \"ElasticNet\", \"Random Forest\", \"XGBClassifier\"]]\n",
    "feature_importance_df = feature_importance_df.sort_values(\"Random Forest\", ascending=False)\n",
    "feature_importance_df.to_csv(os.path.join(dir_to_save_files, \"model_feature_importance.txt\"), sep=\",\", header=True, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(models)):\n",
    "    if model_names[i] == 'Random Forest':\n",
    "        print(\"Random Forest Feature Importances:\")\n",
    "        display(pd.DataFrame(zip(df.columns.values[0:], np.transpose(models[i].feature_importances_))).sort_values([1], ascending=False))\n",
    "        pd.DataFrame(zip(df.columns.values[0:], np.transpose(models[i].feature_importances_))).sort_values([1], ascending=False).to_csv(os.path.join(dir_to_save_files, \"random_forest_feature_importance.txt\"), sep=\"|\", index=False,header=False)\n",
    "    elif model_names[i] == 'Log. Regression':\n",
    "        print(\"Log. Regression Coef. Weights:\")\n",
    "        display(pd.DataFrame(zip(df.columns.values[0:], np.transpose(models[i].coef_[0]))).sort_values([1], ascending=False))\n",
    "        pd.DataFrame(zip(df.columns.values[0:], np.transpose(models[i].coef_[0]))).sort_values([1], ascending=False).to_csv(os.path.join(dir_to_save_files, \"log_reg_feature_importance.txt\"), sep=\"|\", index=False,header=False)\n",
    "    elif model_names[i] == 'ElasticNet':\n",
    "        print(\"ElasticNet Coef. Weights:\")\n",
    "        display(pd.DataFrame(zip(df.columns.values[0:], np.transpose(models[i].coef_[0]))).sort_values([1], ascending=False)) \n",
    "        pd.DataFrame(zip(df.columns.values[0:], np.transpose(models[i].coef_[0]))).sort_values([1], ascending=False).to_csv(os.path.join(dir_to_save_files, \"elasticnet_feature_importance.txt\"), sep=\"|\", index=False,header=False)\n",
    "    elif model_names[i] == 'XGBClassifier':\n",
    "        print(\"XGBoost Feature Importances:\")\n",
    "        display(pd.DataFrame(zip(df.columns.values[0:], np.transpose(models[i].feature_importances_))).sort_values([1], ascending=False))\n",
    "        pd.DataFrame(zip(df.columns.values[0:], np.transpose(models[i].feature_importances_))).sort_values([1], ascending=False).to_csv(os.path.join(dir_to_save_files, \"xgboost_feature_importance.txt\"), sep=\"|\", index=False,header=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Create bar plots for feature importance\n",
    "# fig, ax = plt.subplots(figsize=(10,70))\n",
    "# eldf = zip(df.columns.values[1:], np.transpose(models[model_names.index(\"ElasticNet\")].coef_[0]))\n",
    "# eldf = sorted(eldf, key=lambda x: x[1], reverse=False)\n",
    "\n",
    "# ax.barh(range(len(eldf)), [x[1] for x in eldf])\n",
    "# ax.set_yticks(np.arange(len(eldf)))\n",
    "# ax.set_yticklabels([x[0] for x in eldf])\n",
    "# ax.set_xlabel(\"Model Weight\")\n",
    "# plt.gca().yaxis.grid(True)\n",
    "# plt.tight_layout()\n",
    "# plt.savefig('elasticnet_feature_importance_plot.png', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate class probabilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calibrate_predicted_probability(pred_prob, p_1):\n",
    "#def calibrate_predicted_probability(pred_prob, orig_freq, oversampled_freq):\n",
    "    if pred_prob == 0.0: \n",
    "        return 0.0\n",
    "    p_0 = 1.0 - p_1\n",
    "    #return 1./(1.+((1./orig_freq)-1.)/((1./oversampled_freq)-1.)*((1./pred_prob)-1.))\n",
    "    # formula below from Christine Lee's paper\n",
    "    return 1.0/(1.0+((1.0/pred_prob)-1.0)*(p_0/p_1))\n",
    "\n",
    "calibrate_predicted_probability_vec = np.vectorize(calibrate_predicted_probability)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#print calibrate_predicted_probability_vec([0.09500048], y_test.mean(), y_train.mean())\n",
    "\n",
    "def generate_class_probs(models, X_test):\n",
    "    try:\n",
    "        # generate class probabilities\n",
    "        #model_probs = [calibrate_predicted_probability_vec(mod.predict_proba(X_test), y_test.mean()) for mod in models]\n",
    "        model_probs = [mod.predict_proba(X_test) for mod in models]\n",
    "        #probs = model2.predict(X_test)\n",
    "        #print model_probs\n",
    "    except AttributeError as e:\n",
    "        print(e)\n",
    "        print(mod)\n",
    "        pass\n",
    "    return model_probs\n",
    "\n",
    "model_probs = generate_class_probs(models, np.array(X_test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# death_indices = [i for i,x in enumerate(y_test) if x]\n",
    "# print model_probs[2][death_indices][:,1]\n",
    "# plt.hist(model_probs[2][:,1], bins=100, density=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predict class labels & plot mean square error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# predict class labels for the test set\n",
    "def predict_given_threshold(probs, threshold=0.15):\n",
    "    return [True if x[1] > threshold else False for x in probs]\n",
    "\n",
    "#model_predictions = [np.rint(mod.predict(X_test)) for mod in models]\n",
    "#model_predictions = [predict_given_threshold(probs, y_test.mean()) for probs in model_probs]\n",
    "model_predictions = [mod.predict(np.array(X_test)) for mod in models]\n",
    "print(model_predictions[0][0:10])\n",
    "print(model_probs[0][0:10])\n",
    "print(y_test[0:10])\n",
    "#print np.rint(model_predictions[0][0:10])\n",
    "#mse = [metrics.mean_squared_error(y_test, predicted) for predicted in model_predictions]\n",
    "#print [x[1] for x in model_probs[0]]\n",
    "\n",
    "# def plot_mean_square_error(models, model_names, mse):\n",
    "#     print \"Mean Square Error:\"\n",
    "#     for i in range(len(models)):\n",
    "#         print model_names[i],\"\\t\\t\", mse[i]\n",
    "\n",
    "#     plt.figure(figsize=default_size)\n",
    "#     plt.gca().yaxis.grid(True)\n",
    "#     plt.bar(np.arange(len(models)), mse, align='center')\n",
    "#     plt.xticks(np.arange(len(models)), model_names)\n",
    "#     plt.xlabel('Classifier')\n",
    "#     plt.ylabel('Mean Square Error (MSE)')\n",
    "#     plt.title('Classifier MSE')\n",
    "#     plt.legend(loc=\"lower right\")\n",
    "#     plt.show()\n",
    "    \n",
    "#plot_mean_square_error(models, model_names, mse)\n",
    "#print metrics.mean_squared_error(y_test, [x[1] for x in model_probs[0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate accuracy and AUC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plot_accuracy_roc_auc(models, model_names, model_predictions, model_probs, y_train, y_test, \"accuracy_roc_auc.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot ROC curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_roc_curve(models, model_names, model_probs, y_test, os.path.join(dir_to_save_files,\"roc_curve.png\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "confusion_matrix_classification_rpt(models, model_names, model_predictions, y_test, dir_to_save_files)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot precision-recall curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_precision_recall_curve(models, model_names, model_probs, y_test, os.path.join(dir_to_save_files, \"precision_recall_curve.png\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ranking_precision_score(y_true, y_pred, k=10):\n",
    "    \"\"\"Precision at rank k\n",
    "    Parameters\n",
    "    ----------\n",
    "    y_true : array-like, shape = [n_samples]\n",
    "        Ground truth (true relevance labels).\n",
    "    y_pred : array-like, shape = [n_samples]\n",
    "        Predicted scores.\n",
    "    k : int\n",
    "        Rank.\n",
    "    Returns\n",
    "    -------\n",
    "    precision @k : float\n",
    "    \"\"\"\n",
    "    unique_y = np.unique(y_true)\n",
    "\n",
    "    if len(unique_y) > 2:\n",
    "        raise ValueError(\"Only supported for two relevance levels.\")\n",
    "\n",
    "    pos_label = unique_y[1]\n",
    "    n_pos = np.sum(y_true == pos_label)\n",
    "\n",
    "    order = np.argsort(y_pred)[::-1]\n",
    "    y_true = np.take(y_true, order[:k])\n",
    "    n_relevant = np.sum(y_true == pos_label)\n",
    "\n",
    "    # Divide by min(n_pos, k) such that the best achievable score is always 1.0.\n",
    "    return float(n_relevant) / min(n_pos, k)\n",
    "\n",
    "\n",
    "def plot_precision_at_k(models, model_names, model_probs, y_test, filename):\n",
    "    plt.figure(figsize=default_size)\n",
    "    lw = 2\n",
    "    for i in range(len(model_names)):\n",
    "        k_vals = np.arange(2, 100)\n",
    "        if isinstance(y_test, list):\n",
    "            precision_at_k = [ranking_precision_score(np.array(y_test[i]), model_probs[i][:,1], k=k) for k in k_vals]\n",
    "        else:\n",
    "            precision_at_k = [ranking_precision_score(np.array(y_test), model_probs[i][:,1], k=k) for k in k_vals]\n",
    "        #plt.plot(fpr, tpr, PLOT_COLOURS[i], lw=lw, label=model_names[i]+' (AUC = %0.4f)' % roc_auc)\n",
    "        plt.plot(k_vals, precision_at_k, PLOT_COLOURS[i], lw=lw, label=model_names[i])\n",
    "    #plt.plot(fpr, tpr, color='darkorange',\n",
    "    #         lw=lw, label='ROC curve (area = %0.4f)' % roc_auc)\n",
    "#    plt.xlim([0.0, 1.05])\n",
    "    plt.ylim([0.0, 1.05])\n",
    "    plt.yticks(np.arange(0., 1.1, 0.1))\n",
    "#    plt.xticks(np.arange(0., 1.1, 0.1))\n",
    "    ax = plt.axes()\n",
    "    # default width = 2, def length = 6\n",
    "    ax.set_yticks(np.arange(0., 1.1, 0.1), minor=True)\n",
    "#    ax.set_xticks(np.arange(0., 1.1, 0.1), minor=True)\n",
    "    ax.tick_params(direction='out', length=6, width=0.25, colors='black',labelsize=label_text_size)\n",
    "    ax.tick_params(axis = 'both', which = 'minor', width=0.25)\n",
    "    plt.xlabel('k', fontsize=label_text_size)\n",
    "    plt.ylabel('Precision', fontsize=label_text_size)\n",
    "    #plt.title('Receiver Operating Characteristic')\n",
    "    plt.legend(loc=\"lower right\",fontsize=legend_text_size)\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(filename, format=\"tif\", dpi=fig_dpi)\n",
    "    plt.show()\n",
    "\n",
    "plot_precision_at_k(models, model_names, model_probs, y_test, os.path.join(dir_to_save_files, \"precision_at_k.png\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot calibration curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot derived from: http://scikit-learn.org/stable/auto_examples/calibration/plot_compare_calibration.html#sphx-glr-auto-examples-calibration-plot-compare-calibration-py\n",
    "plt.figure(figsize=(10, 10))\n",
    "ax1 = plt.subplot2grid((3, 1), (0, 0), rowspan=2)\n",
    "ax2 = plt.subplot2grid((3, 1), (2, 0))\n",
    "\n",
    "ax1.plot([0, 1], [0, 1], \"k:\", label=\"Perfectly calibrated\")\n",
    "for i in range(len(models)):\n",
    "    #clf.fit(X_train, y_train)\n",
    "    #if hasattr(models[i], \"predict_proba\"):\n",
    "    #    prob_pos = clf.predict_proba(X_test)[:, 1]\n",
    "    #else:  # use decision function\n",
    "        #prob_pos = clf.decision_function(X_test)\n",
    "        #prob_pos = \\\n",
    "        #    (prob_pos - prob_pos.min()) / (prob_pos.max() - prob_pos.min())\n",
    "    calibrated_model_probs = [calibrate_predicted_probability(j, y_test.mean()) for j in model_probs[i][:, 1]]\n",
    "    print(model_probs[i][:, 1][0:7])\n",
    "    print(calibrated_model_probs[0:7])\n",
    "    #calibrated_model_probs = calibrate_predicted_probability_vec(model_probs[i][:, 1], y_test.mean(), y_train.mean())\n",
    "    fraction_of_positives, mean_predicted_value = \\\n",
    "        calibration_curve(y_test, calibrated_model_probs, n_bins=10)\n",
    "\n",
    "    ax1.plot(mean_predicted_value, fraction_of_positives, PLOT_COLOURS[i], \"s-\",\n",
    "             label=\"%s\" % (model_names[i], ))\n",
    "\n",
    "    ax2.hist(model_probs[i], range=(0, 1), bins=10, label=model_names[i],\n",
    "             histtype=\"step\", lw=2)\n",
    "\n",
    "ax1.set_ylabel(\"Fraction of positives\")\n",
    "ax1.set_ylim([-0.05, 1.05])\n",
    "ax1.legend(loc=\"upper left\")\n",
    "ax1.set_title('Calibration plots  (reliability curve)')\n",
    "\n",
    "ax2.set_xlabel(\"Mean predicted value\")\n",
    "ax2.set_ylabel(\"Count\")\n",
    "ax2.legend(loc=\"upper center\", ncol=2)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(os.path.join(dir_to_save_files, \"calibration_plots.png\"), format=\"png\", dpi=fig_dpi)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate Brier Score Loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Brier Score Loss:\")\n",
    "for i in range(len(models)):\n",
    "    \n",
    "    brier_loss = brier_score_loss(y_test, model_probs[i][:, 1])\n",
    "    print(model_names[i],\"\\t\\t\", brier_loss)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Write predicted probabilities to file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(models)):\n",
    "    print(model_names[i])\n",
    "    X_test_ids[model_names[i] + \"_probs\"] = model_probs[i][:,1]\n",
    "    X_test_ids[model_names[i] + \"_preds\"] = model_predictions[i]\n",
    "X_test_ids.to_csv(os.path.join(dir_to_save_files, \"model_pred_probs_v3.txt\"), sep=\"\\t\", header=True, index=False)\n",
    "if verbose:\n",
    "    X_test_ids.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bootstrap the predicted probabilities, grouping by patient to preserve correlation structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Until we have _n_ predictions, sample patients (with replacement) and keep count of how many total predictions we've sampled\n",
    "2. Repeat (1) for 1000 trials, so that we have 1000 _n_ predictions\n",
    "3. For each of the 1000 trials, compute the metrics of interest\n",
    "4. For each metric of interest, sort the values and take the 25th and 975th values as the lower and upper confidence intervals. Take the mean of all values to get the bootstrapped mean. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(X_test_ids.columns.values)\n",
    "X_test_ids.set_index(keys=\"PAT_ID\", inplace=True, drop=False)\n",
    "X_test_ids = X_test_ids.sort_index()\n",
    "print(X_test_ids.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from tqdm import tqdm_notebook\n",
    "import time\n",
    "n = X_test_ids.shape[0]\n",
    "#n = 100\n",
    "print(\"n:\", n)\n",
    "\n",
    "bootstrap_samples = pd.DataFrame()\n",
    "trials = 100\n",
    "\n",
    "# cache for performance\n",
    "unique_pat_ids = X_test_ids[\"PAT_ID\"].unique()\n",
    "\n",
    "for trial in range(trials):\n",
    "    bootstrap_trial = []\n",
    "    # continue until we have n predictions for this trial\n",
    "    i = 0\n",
    "    num_surg = 0\n",
    "    #while bootstrap_trial.shape[0] < n:\n",
    "    while num_surg < n:\n",
    "        t0 = time.time()\n",
    "        # get patient from list of unique patients\n",
    "        sample_patient_id = np.random.choice(unique_pat_ids)\n",
    "        t1 = time.time()\n",
    "        sample_preds = X_test_ids[X_test_ids.PAT_ID.isin([sample_patient_id])]\n",
    "#        sample_preds = X_test_ids.loc[[sample_patient_id]]\n",
    "#        print(sample_preds.values)\n",
    "#        print(sample_preds.values.shape[0])\n",
    "        t2 = time.time()\n",
    "#        bootstrap_trial.append(sample_preds.values)\n",
    "        for p in sample_preds.values:\n",
    "#             print(p)\n",
    "            bootstrap_trial.append(p)\n",
    "#         print(\"append:\", time.time()-t2)\n",
    "#         print(\"random choice:\", t1-t0)\n",
    "#         print(\"data access:\", t2-t1)\n",
    "        i += 1\n",
    "        num_surg += sample_preds.shape[0]\n",
    "    bootstrap_trial = pd.DataFrame.from_records(bootstrap_trial, columns=X_test_ids.columns)\n",
    "    bootstrap_trial[\"trial_num\"] = trial\n",
    "    print(bootstrap_trial.shape)\n",
    "    bootstrap_samples = bootstrap_samples.append(bootstrap_trial)\n",
    "bootstrap_samples.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# calculate model performance metrics and save to file\n",
    "\n",
    "def tp(y_true, y_pred): return metrics.confusion_matrix(y_true, y_pred)[1, 1]\n",
    "def tn(y_true, y_pred): return metrics.confusion_matrix(y_true, y_pred)[0, 0]\n",
    "def fp(y_true, y_pred): return metrics.confusion_matrix(y_true, y_pred)[0, 1]\n",
    "def fn(y_true, y_pred): return metrics.confusion_matrix(y_true, y_pred)[1, 0]\n",
    "\n",
    "def ci(results, alpha=0.05):\n",
    "    results = sorted(results)\n",
    "    max_index = int(np.floor((1.-alpha/2.)*len(results)))\n",
    "    min_index = int(np.floor((alpha/2.)*len(results)))\n",
    "    #print(len(results))\n",
    "    #print(\"min index:\", min_index)\n",
    "    #print(\"max_index:\", max_index)\n",
    "    return results[min_index], results[max_index]\n",
    "\n",
    "metrics_dict = {}\n",
    "metric_string = \"{:0.3f} ({:0.3f}-{:0.3f})\"\n",
    "\n",
    "for i in range(len(model_names)):\n",
    "    metrics_dict[model_names[i]] = []\n",
    "    \n",
    "    model_probs = bootstrap_samples.groupby(\"trial_num\")[model_names[i] + \"_probs\"].apply(list)\n",
    "    model_preds = bootstrap_samples.groupby(\"trial_num\")[model_names[i] + \"_preds\"].apply(list)\n",
    "    y_test = bootstrap_samples.groupby(\"trial_num\")[\"INPT_DEATH_YN\"].apply(list)\n",
    "    \n",
    "#     model_probs = model_probs_cv[i]\n",
    "#     model_preds = model_preds_cv[i]\n",
    "#     y_test = testing_labels_cv[i]\n",
    "    \n",
    "    tp_vec = [tp(y_true, y_pred) for (y_true, y_pred) in zip(y_test, model_preds)]\n",
    "    tn_vec = [tn(y_true, y_pred) for (y_true, y_pred) in zip(y_test, model_preds)]\n",
    "    fp_vec = [fp(y_true, y_pred) for (y_true, y_pred) in zip(y_test, model_preds)]\n",
    "    fn_vec = [fn(y_true, y_pred) for (y_true, y_pred) in zip(y_test, model_preds)]\n",
    "    \n",
    "    fpr_vec = np.array([metrics.roc_curve(y_true, probs)[0] for (y_true, probs) in zip(y_test, model_probs)])\n",
    "    tpr_vec = np.array([metrics.roc_curve(y_true, probs)[1] for (y_true, probs) in zip(y_test, model_probs)])\n",
    "    roc_auc_vec = np.array([metrics.auc(fpr, tpr) for (fpr, tpr) in zip(fpr_vec, tpr_vec)])\n",
    "    precision_vec = np.array([metrics.precision_score(y_true, y_pred) for (y_true, y_pred) in zip (y_test, model_preds)])\n",
    "    recall_vec = np.array([metrics.recall_score(y_true, y_pred) for (y_true, y_pred) in zip (y_test, model_preds)])\n",
    "    f1_vec = np.array([metrics.f1_score(y_true, y_pred) for (y_true, y_pred) in zip (y_test, model_preds)])\n",
    "    accuracy_vec = np.array([metrics.accuracy_score(y_true, y_pred) for (y_true, y_pred) in zip (y_test, model_preds)])\n",
    "    specificity_vec = np.array([(float(t_n)/(t_n + f_p)) for (t_n, f_p) in zip(tn_vec, fp_vec)])\n",
    "    brier_score_vec = np.array([brier_score_loss(y_true, probs) for (y_true, probs) in zip(y_test, model_probs)])\n",
    "    print(\"==================================================\")\n",
    "    print(\"MODEL: \", model_names[i])\n",
    "    print(\"ROC AUC:\\t\", metric_string.format(roc_auc_vec.mean(), ci(roc_auc_vec)[0], ci(roc_auc_vec)[1]))\n",
    "    print(\"specificity:\\t\", metric_string.format(specificity_vec.mean(), ci(specificity_vec)[0], ci(specificity_vec)[1]))\n",
    "    print(\"precision:\\t\", metric_string.format(precision_vec.mean(), ci(precision_vec)[0], ci(precision_vec)[1]))\n",
    "    print(\"recall:\\t\\t\", metric_string.format(recall_vec.mean(), ci(recall_vec)[0], ci(recall_vec)[1]))\n",
    "    print(\"f1 score:\\t\", metric_string.format(f1_vec.mean(), ci(f1_vec)[0], ci(f1_vec)[1]))\n",
    "    print(\"accuracy:\\t\", metric_string.format(accuracy_vec.mean(), ci(accuracy_vec)[0], ci(accuracy_vec)[1]))\n",
    "    print(\"brier score:\\t\", metric_string.format(brier_score_vec.mean(), ci(brier_score_vec)[0], ci(brier_score_vec)[1]))\n",
    "    \n",
    "    metrics_dict[model_names[i]].append(metric_string.format(roc_auc_vec.mean(), ci(roc_auc_vec)[0], ci(roc_auc_vec)[1]))\n",
    "    metrics_dict[model_names[i]].append(metric_string.format(specificity_vec.mean(), ci(specificity_vec)[0], ci(specificity_vec)[1]))\n",
    "    metrics_dict[model_names[i]].append(metric_string.format(precision_vec.mean(), ci(precision_vec)[0], ci(precision_vec)[1]))\n",
    "    metrics_dict[model_names[i]].append(metric_string.format(recall_vec.mean(), ci(recall_vec)[0], ci(recall_vec)[1]))\n",
    "    metrics_dict[model_names[i]].append(metric_string.format(f1_vec.mean(), ci(f1_vec)[0], ci(f1_vec)[1]))\n",
    "    metrics_dict[model_names[i]].append(metric_string.format(accuracy_vec.mean(), ci(accuracy_vec)[0], ci(accuracy_vec)[1]))\n",
    "    metrics_dict[model_names[i]].append(metric_string.format(brier_score_vec.mean(), ci(brier_score_vec)[0], ci(brier_score_vec)[1]))\n",
    "    \n",
    "metrics_df = pd.DataFrame.from_dict(metrics_dict, orient=\"index\")\n",
    "metrics_df.columns =[\"ROC AUC\", \"Specificity\", \"Precision\", \"Recall\", \"F1 Score\", \"Accuracy\", \"Brier Score\"]\n",
    "metrics_df = metrics_df[[\"Accuracy\", \"F1 Score\", \"Precision\", \"Recall\", \"Specificity\", \"ROC AUC\", \"Brier Score\"]]\n",
    "metrics_df = metrics_df.transpose()\n",
    "metrics_df = metrics_df[[\"Log. Regression\", \"ElasticNet\", \"Random Forest\", \"XGBClassifier\"]]\n",
    "metrics_df.to_csv(os.path.join(dir_to_save_files, \"model_scores_manual_cv.txt\"), sep=\",\", header=True, index=True)\n",
    "metrics_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "foo bar baz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross Validation of Model Metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_folds = 5\n",
    "\n",
    "group_by_admsn = admsn_ids.astype('category').cat.codes\n",
    "print \"number of operations:\", len(group_by_admsn)\n",
    "print \"number of unique admission groups:\", len(np.unique(group_by_admsn))\n",
    "group_kfold = GroupKFold(n_splits=num_folds)\n",
    "\n",
    "# save GroupKFold object as pickle object\n",
    "pickle.dump(group_kfold, open( os.path.join(dir_to_save_files, \"group_kfold.p\"), \"wb\" ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# keep track of these during cross-validation\n",
    "model_probs_cv = []\n",
    "model_preds_cv = []\n",
    "testing_labels_cv = []\n",
    "admsn_ids_cv = []\n",
    "or_case_ids_cv = []\n",
    "test_index_cv = []\n",
    "admsn_surgery_num_cv = []\n",
    "for model in model_names:\n",
    "    model_probs_cv.append([])\n",
    "    model_preds_cv.append([])\n",
    "    testing_labels_cv.append([])\n",
    "\n",
    "cv_num = 1\n",
    "for train_index, test_index in group_kfold.split(np.array(df.iloc[:,0:]), y, group_by_admsn):\n",
    "    X_train, X_test = df.iloc[train_index], df.iloc[test_index]\n",
    "    y_train, y_test = y[train_index], y[test_index]\n",
    "    \n",
    "    # keep track of these for writing to file\n",
    "    admsn_ids_cv = admsn_ids_cv + list(admsn_ids.iloc[test_index])\n",
    "    or_case_ids_cv = or_case_ids_cv + list(or_case_id_number.iloc[test_index])\n",
    "    admsn_surgery_num_cv = admsn_surgery_num_cv + list(admsn_surgery_number.iloc[test_index])\n",
    "    test_index_cv = test_index_cv + [cv_num] * len(list(test_index))\n",
    "    \n",
    "    print \"admsn_ids_cv len:\", len(admsn_ids_cv)\n",
    "    print \"or_case_ids_cv len:\", len(or_case_ids_cv)\n",
    "    print \"admsn_surgery_num_cv len:\", len(admsn_surgery_num_cv)\n",
    "    print \"test_index_cv len:\", len(test_index_cv)\n",
    "    \n",
    "    # scale the training and testing data\n",
    "    scaler = StandardizeWithNaN()\n",
    "    # fit on training data\n",
    "    scaler.fit(X_train)\n",
    "    X_train = scaler.transform(X_train) \n",
    "    # apply same transformation to test data\n",
    "    X_test = scaler.transform(X_test) \n",
    "    \n",
    "    print(\"imputing X_train\")\n",
    "    si = SoftImputeEstimator(verbose=False)\n",
    "    X_train = si.transform(X_train)\n",
    "    print(\"imputing X_test\")\n",
    "    X_test = si.transform(X_test.replace(np.inf, np.nan))\n",
    "    \n",
    "    # oversample the training dataset\n",
    "    X_train, y_train = SMOTE(kind='borderline1',k_neighbors=3).fit_sample(X_train, y_train)\n",
    "    \n",
    "    i = 0\n",
    "    for mod in models:\n",
    "        print(\"fitting {} model\".format(model_names[i]))\n",
    "        mod.fit(X_train, y_train)\n",
    "        print(\"finished fitting...\")\n",
    "        model_probs = mod.predict_proba(np.array(X_test))\n",
    "        model_predictions = mod.predict(np.array(X_test))\n",
    "        model_probs_cv[i].append(model_probs)\n",
    "        model_preds_cv[i].append(model_predictions)\n",
    "        testing_labels_cv[i].append(y_test)\n",
    "        i += 1\n",
    "    cv_num += 1\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate model performance metrics and save to file\n",
    "\n",
    "def tp(y_true, y_pred): return metrics.confusion_matrix(y_true, y_pred)[1, 1]\n",
    "def tn(y_true, y_pred): return metrics.confusion_matrix(y_true, y_pred)[0, 0]\n",
    "def fp(y_true, y_pred): return metrics.confusion_matrix(y_true, y_pred)[0, 1]\n",
    "def fn(y_true, y_pred): return metrics.confusion_matrix(y_true, y_pred)[1, 0]\n",
    "\n",
    "def ci(results):\n",
    "    return st.t.interval(0.95, len(results)-1, loc=np.mean(results), scale=st.sem(results))\n",
    "\n",
    "metrics_dict = {}\n",
    "metric_string = \"{:0.3f} ({:0.3f}-{:0.3f})\"\n",
    "\n",
    "for i in range(len(model_names)):\n",
    "    metrics_dict[model_names[i]] = []\n",
    "    \n",
    "    model_probs = model_probs_cv[i]\n",
    "    model_preds = model_preds_cv[i]\n",
    "    y_test = testing_labels_cv[i]\n",
    "    \n",
    "    tp_vec = [tp(y_true, y_pred) for (y_true, y_pred) in zip(y_test, model_preds)]\n",
    "    tn_vec = [tn(y_true, y_pred) for (y_true, y_pred) in zip(y_test, model_preds)]\n",
    "    fp_vec = [fp(y_true, y_pred) for (y_true, y_pred) in zip(y_test, model_preds)]\n",
    "    fn_vec = [fn(y_true, y_pred) for (y_true, y_pred) in zip(y_test, model_preds)]\n",
    "    \n",
    "    fpr_vec = np.array([metrics.roc_curve(y_true, probs[:, 1])[0] for (y_true, probs) in zip(y_test, model_probs)])\n",
    "    tpr_vec = np.array([metrics.roc_curve(y_true, probs[:, 1])[1] for (y_true, probs) in zip(y_test, model_probs)])\n",
    "    roc_auc_vec = np.array([metrics.auc(fpr, tpr) for (fpr, tpr) in zip(fpr_vec, tpr_vec)])\n",
    "    precision_vec = np.array([metrics.precision_score(y_true, y_pred) for (y_true, y_pred) in zip (y_test, model_preds)])\n",
    "    recall_vec = np.array([metrics.recall_score(y_true, y_pred) for (y_true, y_pred) in zip (y_test, model_preds)])\n",
    "    f1_vec = np.array([metrics.f1_score(y_true, y_pred) for (y_true, y_pred) in zip (y_test, model_preds)])\n",
    "    accuracy_vec = np.array([metrics.accuracy_score(y_true, y_pred) for (y_true, y_pred) in zip (y_test, model_preds)])\n",
    "    specificity_vec = np.array([(float(t_n)/(t_n + f_p)) for (t_n, f_p) in zip(tn_vec, fp_vec)])\n",
    "    print \"==================================================\"\n",
    "    print \"MODEL: \", model_names[i]\n",
    "    print \"ROC AUC:\\t\", metric_string.format(roc_auc_vec.mean(), ci(roc_auc_vec)[0], ci(roc_auc_vec)[1])\n",
    "    print \"specificity:\\t\", metric_string.format(specificity_vec.mean(), ci(specificity_vec)[0], ci(specificity_vec)[1])\n",
    "    print \"precision:\\t\", metric_string.format(precision_vec.mean(), ci(precision_vec)[0], ci(precision_vec)[1])\n",
    "    print \"recall:\\t\\t\", metric_string.format(recall_vec.mean(), ci(recall_vec)[0], ci(recall_vec)[1])\n",
    "    print \"f1 score:\\t\", metric_string.format(f1_vec.mean(), ci(f1_vec)[0], ci(f1_vec)[1])\n",
    "    print \"accuracy:\\t\", metric_string.format(accuracy_vec.mean(), ci(accuracy_vec)[0], ci(accuracy_vec)[1])\n",
    "    \n",
    "    metrics_dict[model_names[i]].append(metric_string.format(roc_auc_vec.mean(), ci(roc_auc_vec)[0], ci(roc_auc_vec)[1]))\n",
    "    metrics_dict[model_names[i]].append(metric_string.format(specificity_vec.mean(), ci(specificity_vec)[0], ci(specificity_vec)[1]))\n",
    "    metrics_dict[model_names[i]].append(metric_string.format(precision_vec.mean(), ci(precision_vec)[0], ci(precision_vec)[1]))\n",
    "    metrics_dict[model_names[i]].append(metric_string.format(recall_vec.mean(), ci(recall_vec)[0], ci(recall_vec)[1]))\n",
    "    metrics_dict[model_names[i]].append(metric_string.format(f1_vec.mean(), ci(f1_vec)[0], ci(f1_vec)[1]))\n",
    "    metrics_dict[model_names[i]].append(metric_string.format(accuracy_vec.mean(), ci(accuracy_vec)[0], ci(accuracy_vec)[1]))\n",
    "    \n",
    "metrics_df = pd.DataFrame.from_dict(metrics_dict, orient=\"index\")\n",
    "metrics_df.columns =[\"ROC AUC\", \"Specificity\", \"Precision\", \"Recall\", \"F1 Score\", \"Accuracy\"]\n",
    "metrics_df.head()\n",
    "metrics_df.to_csv(os.path.join(dir_to_save_files, \"model_scores_manual_cv.txt\"), sep=\",\", header=True, index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_probs_dict = {\"OR_CASE_ID\": or_case_ids_cv,\n",
    "                    \"ADMSN_ID\": admsn_ids_cv, \n",
    "                    \"ADMSN_SURGERY_NUMBER\": admsn_surgery_num_cv, \n",
    "                    \"test_index\": test_index_cv,\n",
    "                    \"INPT_DEATH_YN\": []}\n",
    "\n",
    "for model_name in model_names:\n",
    "    model_probs_dict[model_name] = []\n",
    "\n",
    "for j in range(group_kfold.get_n_splits()):\n",
    "    for i in range(len(model_names)):\n",
    "        print(len(model_probs_cv[i][j]))\n",
    "        model_probs_dict[model_names[i]] = model_probs_dict[model_names[i]] + list(model_probs_cv[i][j][:,1])\n",
    "        if i == 0:\n",
    "            model_probs_dict[\"INPT_DEATH_YN\"] = model_probs_dict[\"INPT_DEATH_YN\"] + list(testing_labels_cv[i][j])\n",
    "            \n",
    "for k,v in model_probs_dict.iteritems():\n",
    "    print k, len(v)\n",
    "            \n",
    "model_probs_df = pd.DataFrame(data=model_probs_dict)\n",
    "model_probs_df.to_csv(os.path.join(dir_to_save_files, \"model_pred_probs_v2.txt\"), sep=\"\\t\", header=True, index=False)\n",
    "if verbose:\n",
    "    model_probs_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "foobar baz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# take the first fold of the k-fold split\n",
    "# train_index, test_index = group_kfold.split(transformed_X, y, group_by_admsn).next()\n",
    "# X_train, X_test = transformed_X[train_index], transformed_X[test_index]\n",
    "# y_train, y_test = y[train_index], y[test_index]\n",
    "\n",
    "# ddd = {\"OR_CASE_ID\": or_case_id_number,\n",
    "#        \"ADMSN_ID\": admsn_ids, \n",
    "#        \"ADMSN_SURGERY_NUMBER\": admsn_surgery_number,\n",
    "#        \"INPT_DEATH_YN\": y}\n",
    "\n",
    "cv = StratifiedKFold(n_splits=num_folds, shuffle=True, random_state=42)\n",
    "model_probs_cv = []\n",
    "model_count = 0\n",
    "for model in models_cv:\n",
    "    print(\"starting model\", model_count)\n",
    "    probs = cross_val_predict(model, df.iloc[:,0:], y, \n",
    "                      cv=group_kfold.split(df.iloc[:,0:], y, group_by_admsn), \n",
    "                      n_jobs=1, method='predict_proba', verbose=3) \n",
    "    print(\"finished model\", model_count)\n",
    "#     ddd[model_names[model_count]] = model_probs_cv[model_count][:,1]\n",
    "#     new_df = pd.DataFrame(data=ddd)\n",
    "#     new_df.to_csv(os.path.join(dir_to_save_files, \"model_pred_probs.txt\"), sep=\"\\t\", header=True, index=False)\n",
    "    model_count += 1\n",
    "    model_probs_cv.append(probs)\n",
    "\n",
    "#model_probs_cv = [cross_val_predict(model, df.iloc[:,0:], y, cv=group_kfold.split(df.iloc[:,0:], y, group_by_admsn), n_jobs=5, method='predict_proba') for model in models_cv]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ddd = {\"ADMSN_ID\": admsn_ids, \n",
    "#        \"PRED_DEATH_LOGREG\": model_probs_cv[0][:,1], \n",
    "#        \"PRED_DEATH_ELASTICNET\": model_probs_cv[1][:,1], \n",
    "#        \"PRED_DEATH_RF\": model_probs_cv[2][:,1],\n",
    "#        \"PRED_DEATH_XGB\": model_probs_cv[3][:,1],\n",
    "#        \"INPT_DEATH_YN\": y}\n",
    "ddd = {\"OR_CASE_ID\": or_case_id_number,\n",
    "       \"ADMSN_ID\": admsn_ids, \n",
    "       \"ADMSN_SURGERY_NUMBER\": admsn_surgery_number,\n",
    "       \"INPT_DEATH_YN\": y}\n",
    "for i in range(len(model_names)):\n",
    "    ddd[model_names[i]] = model_probs_cv[i][:,1]\n",
    "new_df = pd.DataFrame(data=ddd)\n",
    "new_df.to_csv(os.path.join(dir_to_save_files, \"model_pred_probs.txt\"), sep=\"\\t\", header=True, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# new_df.sort_values(by=[\"ADMSN_ID\", \"ADMSN_SURGERY_NUMBER\"], inplace=True)\n",
    "# new_df.head(5)\n",
    "# print new_df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# print admsn_ids.shape\n",
    "# print transformed_X.shape\n",
    "# lldf = pd.DataFrame(data=np.concatenate((np.expand_dims(admsn_ids, axis=1), transformed_X), axis=1), columns=[\"ADMSN_ID\"] + list(df.columns.values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# total_admissions = new_df[\"ADMSN_ID\"].nunique()\n",
    "# num_mult_admission = 0\n",
    "# num_live_increase_risk_2nd_surgery = 0\n",
    "# num_live_increase_risk_last_surgery = 0\n",
    "\n",
    "# num_death = 0\n",
    "# num_death_increase_risk_2nd_surgery = 0\n",
    "# num_death_increase_risk_last_surgery = 0\n",
    "\n",
    "# for name,group in new_df.groupby(\"ADMSN_ID\"):\n",
    "# #     print name[0:10]\n",
    "# #     print group.shape\n",
    "#     # if we have patients with multiple surgeries\n",
    "#     if group.shape[0] > 1:\n",
    "#         num_mult_admission += 1\n",
    "#         #print group\n",
    "#         if group[\"INPT_DEATH_YN\"].iloc[0] == False:\n",
    "#             # check if risk increased from first to second surgery\n",
    "#             if group[\"PRED_DEATH_RF\"].iloc[1] > group[\"PRED_DEATH_RF\"].iloc[0]:\n",
    "#                 num_live_increase_risk_2nd_surgery += 1\n",
    "#             # check if risk increased from first to last surgery\n",
    "#             if group[\"PRED_DEATH_RF\"].iloc[-1] > group[\"PRED_DEATH_RF\"].iloc[0]:\n",
    "#                 num_live_increase_risk_last_surgery += 1\n",
    "#         else:\n",
    "#             num_death += 1\n",
    "#             # check if risk increased from first to second surgery\n",
    "#             if group[\"PRED_DEATH_RF\"].iloc[1] > group[\"PRED_DEATH_RF\"].iloc[0]:\n",
    "#                 num_death_increase_risk_2nd_surgery += 1\n",
    "#             # check if risk increased from first to last surgery\n",
    "#             if group[\"PRED_DEATH_RF\"].iloc[-1] > group[\"PRED_DEATH_RF\"].iloc[0]:\n",
    "#                 num_death_increase_risk_last_surgery += 1\n",
    "            \n",
    "# print \"total admissions:\", total_admissions\n",
    "# print \"number of deaths: {} ({}%)\".format(num_death, num_death * 100. / total_admissions)\n",
    "# print \"number of survivors that increased in risk from 1st to 2nd surgery: {} ({}%)\".format(\n",
    "#     num_live_increase_risk_2nd_surgery, num_live_increase_risk_2nd_surgery * 100. / (num_mult_admission-num_death))\n",
    "# print \"number of survivors that increased in risk from 1st to LAST surgery: {} ({}%)\".format(\n",
    "#     num_live_increase_risk_last_surgery, num_live_increase_risk_last_surgery * 100. / (num_mult_admission-num_death))\n",
    "# print \"number of deaths that increased in risk from 1st to 2nd surgery: {} ({}%)\".format(\n",
    "#     num_death_increase_risk_2nd_surgery, num_death_increase_risk_2nd_surgery * 100. / num_death)\n",
    "# print \"number of deaths that increased in risk from 1st to LAST surgery: {} ({}%)\".format(\n",
    "#     num_death_increase_risk_last_surgery, num_death_increase_risk_last_surgery * 100. / num_death)\n",
    "\n",
    "# i = 0\n",
    "# plt.hold(True)\n",
    "# for name,group in new_df.groupby(\"ADMSN_ID\"):\n",
    "# #     print name[0:10]\n",
    "# #     print group.shape\n",
    "#     # if we have patients with multiple surgeries\n",
    "#     if group.shape[0] > 1:\n",
    "#         #print group\n",
    "#         plot_color = 'b' if group[\"INPT_DEATH_YN\"].iloc[0] == False else 'r'\n",
    "#         if group[\"INPT_DEATH_YN\"].iloc[0] == False:\n",
    "#             plt.plot(range(group.shape[0]), group[\"PRED_DEATH_RF\"], plot_color)\n",
    "#             i = i+1\n",
    "#         if i > 100: \n",
    "#             break\n",
    "# i = 0\n",
    "# for name,group in new_df.groupby(\"ADMSN_ID\"):\n",
    "# #     print name[0:10]\n",
    "# #     print group.shape\n",
    "#     # if we have patients with multiple surgeries\n",
    "#     if group.shape[0] > 1:\n",
    "#         #print group\n",
    "#         plot_color = 'b' if group[\"INPT_DEATH_YN\"].iloc[0] == False else 'r'\n",
    "#         if group[\"INPT_DEATH_YN\"].iloc[0] == True:\n",
    "#             plt.plot(range(group.shape[0]), group[\"PRED_DEATH_RF\"], plot_color)\n",
    "#             i = i+1\n",
    "#         if i > 100: \n",
    "#             break\n",
    "# plt.title(\"Change in mortality prediction over multiple surgeries\")\n",
    "# plt.xlabel(\"Surgery Number\")\n",
    "# plt.ylabel(\"P(death)\")\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_roc_curve(models, model_names, model_probs_cv, y, os.path.join(dir_to_save_files, \"roc_curve_cv.tif\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_precision_recall_curve(models, model_names, model_probs_cv, y, os.path.join(dir_to_save_files, \"precision_recall_cv.tif\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print \"Brier Score Loss:\"\n",
    "with open(os.path.join(dir_to_save_files, \"brier_score_cv.txt\"), \"w\") as text_file:\n",
    "    for i in range(len(models)):\n",
    "        brier_loss = brier_score_loss(y, model_probs_cv[i][:, 1])\n",
    "        print model_names[i],\"\\t\\t\", brier_loss\n",
    "        text_file.write(\"{},{}\\n\".format(model_names[i], brier_loss))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_folds = 5\n",
    "\n",
    "# group_by_admsn = admsn_ids.astype('category').cat.codes\n",
    "# print \"number of operations:\", len(group_by_admsn)\n",
    "# print \"number of unique admission groups:\", len(np.unique(group_by_admsn))\n",
    "# group_kfold = GroupKFold(n_splits=num_folds)\n",
    "\n",
    "#cv = StratifiedKFold(n_splits=num_folds, shuffle=True, random_state=42)\n",
    "\n",
    "def tp(y_true, y_pred): return metrics.confusion_matrix(y_true, y_pred)[1, 1]\n",
    "def tn(y_true, y_pred): return metrics.confusion_matrix(y_true, y_pred)[0, 0]\n",
    "def fp(y_true, y_pred): return metrics.confusion_matrix(y_true, y_pred)[0, 1]\n",
    "def fn(y_true, y_pred): return metrics.confusion_matrix(y_true, y_pred)[1, 0]\n",
    "\n",
    "scoring = {'tp' : make_scorer(tp), 'tn' : make_scorer(tn),\n",
    "           'fp' : make_scorer(fp), 'fn' : make_scorer(fn),\n",
    "            'accuracy': 'accuracy', 'roc_auc': 'roc_auc',\n",
    "             'average_precision': 'average_precision',\n",
    "          'precision': 'precision', 'recall': 'recall',\n",
    "          'f1_score': 'f1'}\n",
    "\n",
    "#scoring = ['accuracy', 'roc_auc', 'average_precision']\n",
    "# model_scores_cv = [cross_validate(model, df.iloc[:,0:], y, scoring=scoring, \n",
    "#                                   cv=group_kfold.split(df.iloc[:,0:], y, group_by_admsn), \n",
    "#                                   n_jobs=1, verbose=2) for model in models_cv]\n",
    "\n",
    "model_scores_cv = []\n",
    "model_count = 0\n",
    "for model in models_cv:\n",
    "    print(\"starting model\", model_count)\n",
    "#     print(model.named_steps)\n",
    "    mod_scores = cross_validate(model, df.iloc[:,0:], y, scoring=scoring, \n",
    "                                  cv=group_kfold.split(df.iloc[:,0:], y, group_by_admsn), \n",
    "                                  n_jobs=1, verbose=10)\n",
    "    print(\"finished model\", model_count)\n",
    "    model_count += 1\n",
    "    model_scores_cv.append(mod_scores)\n",
    "\n",
    "# model_accuracy_cv = [cross_val_score(model, X_train, y_train, \n",
    "#                                    scoring='accuracy', cv=cv, n_jobs=-1) for model in models]\n",
    "# model_roc_auc_cv = [cross_val_score(model, X_train, y_train, \n",
    "#                                    scoring='roc_auc', cv=cv, n_jobs=-1) for model in models]\n",
    "# model_avg_precision_cv = [cross_val_score(model, X_train, y_train, \n",
    "#                                    scoring='average_precision', cv=cv, n_jobs=-1) for model in models]\n",
    "print model_scores_cv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=default_size)\n",
    "index = np.arange(len(models))\n",
    "bar_width = 0.35\n",
    "opacity = 0.8\n",
    " \n",
    "rects1 = plt.bar(index, [m['test_accuracy'].mean() for m in model_scores_cv], bar_width,\n",
    "                 alpha=opacity,\n",
    "                 color=PLOT_COLOURS[0],\n",
    "                 yerr=[np.std(m['test_accuracy']) for m in model_scores_cv],\n",
    "                 label='Mean Accuracy')\n",
    " \n",
    "rects2 = plt.bar(index + bar_width, [m['test_roc_auc'].mean() for m in model_scores_cv], bar_width,\n",
    "                 alpha=opacity,\n",
    "                 color=PLOT_COLOURS[1],\n",
    "                 yerr=[np.std(m['test_roc_auc']) for m in model_scores_cv],\n",
    "                 label='Mean ROC AUC')\n",
    "\n",
    "rects3 = plt.bar(index + 2*bar_width, [m['test_average_precision'].mean() for m in model_scores_cv], bar_width,\n",
    "                 alpha=opacity,\n",
    "                 color=PLOT_COLOURS[2],\n",
    "                 yerr=[np.std(m['test_average_precision']) for m in model_scores_cv],\n",
    "                 label='Mean Avg. Precision')\n",
    "\n",
    " \n",
    "plt.xlabel('Model')\n",
    "plt.ylabel('Score')\n",
    "plt.title('Model Accuracy/ROC AUC/Avg. Precision with ASA_STATUS (CV=%d)' % num_folds)\n",
    "plt.xticks(index + bar_width, model_names)\n",
    "plt.ylim([0.0, 1.05])\n",
    "plt.yticks(np.arange(0., 1.1, 0.05))\n",
    "plt.legend(loc=\"lower right\")\n",
    "plt.gca().yaxis.grid(True)\n",
    "plt.plot([1.0 - y_test.mean()]*(len(models)), color='red', lw=2, linestyle='--')\n",
    "plt.tight_layout()\n",
    "plt.savefig(os.path.join(dir_to_save_files, \"accuracy_roc_auc_cv.tif\"), format=\"tif\", dpi=fig_dpi)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ci(results):\n",
    "    return st.t.interval(0.95, len(results)-1, loc=np.mean(results), scale=st.sem(results))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "with open(os.path.join(dir_to_save_files, \"model_scores_cv.txt\"), \"w\") as text_file:\n",
    "    text_file.write(\"Model,\" + \",\".join(scoring.keys()) + \",specificity\" + \"\\n\")\n",
    "    for mod in model_names:\n",
    "        print \"Model:\", mod\n",
    "        text_file.write(mod + \",\")\n",
    "        for score in scoring.keys():\n",
    "            score_vals = model_scores_cv[model_names.index(mod)]['test_' + score]\n",
    "            print \"\\t\", score, \"\\t{:0.3f} ({:0.3f}-{:0.3f})\".format(score_vals.mean(), ci(score_vals)[0], ci(score_vals)[1])\n",
    "            text_file.write(\"{:0.3f} ({:0.3f}-{:0.3f}),\".format(score_vals.mean(), ci(score_vals)[0], ci(score_vals)[1]))\n",
    "            #print \"mod:\", mod, \"score:\", score\n",
    "\n",
    "        tn = model_scores_cv[model_names.index(mod)]['test_tn']\n",
    "        fp = model_scores_cv[model_names.index(mod)]['test_fp']\n",
    "        specificity_vals = np.array(tn, dtype=\"float\") / (np.array(tn) + np.array(fp))\n",
    "        print \"\\tspecificity:\\t\", \"\\t{:0.3f} ({:0.3f}-{:0.3f})\".format(specificity_vals.mean(), ci(specificity_vals)[0], ci(specificity_vals)[1])\n",
    "        text_file.write(\"{:0.3f} ({:0.3f}-{:0.3f})\\n\".format(specificity_vals.mean(), ci(specificity_vals)[0], ci(specificity_vals)[1]))\n",
    "        precision_vals = score_vals = model_scores_cv[model_names.index(mod)]['test_precision']\n",
    "        recall_vals = score_vals = model_scores_cv[model_names.index(mod)]['test_recall']\n",
    "        f1_score_vals = np.nan_to_num(2.0*(precision_vals*recall_vals)/(precision_vals + recall_vals))\n",
    "        print \"\\tf1 score:\\t\", \"\\t{:0.3f} ({:0.3f}-{:0.3f})\".format(f1_score_vals.mean(), ci(f1_score_vals)[0], ci(f1_score_vals)[1])\n",
    "        #text_file.write(\"\\t{:0.3f} ({:0.3f}-{:0.3f})\\n\".format(f1_score_vals.mean(), ci(f1_score_vals)[0], ci(f1_score_vals)[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rf_cv_auc = model_scores_cv[model_names.index('Random Forest')]['test_roc_auc']\n",
    "print \"Random Forest CV AUC\", rf_cv_auc\n",
    "print rf_cv_auc.mean()\n",
    "print ci(rf_cv_auc)\n",
    "rf_cv_precision = model_scores_cv[model_names.index('Random Forest')]['test_precision']\n",
    "print \"Random Forest CV Precision\", rf_cv_precision\n",
    "print rf_cv_precision.mean()\n",
    "print ci(rf_cv_precision)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "en_cv_auc = model_scores_cv[model_names.index('ElasticNet')]['test_roc_auc']\n",
    "print \"ElasticNet CV AUC\", en_cv_auc\n",
    "print en_cv_auc.mean()\n",
    "print ci(en_cv_auc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_cv_auc = model_scores_cv[model_names.index('Log. Regression')]['test_roc_auc']\n",
    "print \"Log. Regression CV AUC\", lr_cv_auc\n",
    "print lr_cv_auc.mean()\n",
    "print ci(lr_cv_auc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lr_cv_auc = model_scores_cv[model_names.index('XGBClassifier')]['test_roc_auc']\n",
    "print \"XGBClassifier CV AUC\", lr_cv_auc\n",
    "print lr_cv_auc.mean()\n",
    "print ci(lr_cv_auc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "######################################################\n",
    "# Code for imputing ASA status\n",
    "######################################################\n",
    "\n",
    "#scaler = StandardScaler() \n",
    "scaler = StandardizeWithNaN()\n",
    "imp = SoftImputeEstimator(verbose=True)\n",
    "smt = SMOTE(kind='borderline1', k_neighbors=3)\n",
    "\n",
    "cv = StratifiedKFold(n_splits=10, shuffle=True, random_state=42)\n",
    "#rfr = RandomForestRegressor(oob_score=True, n_estimators=2000, n_jobs=1, max_features=\"log2\")\n",
    "rfr = XGBRegressor(max_depth=7, n_estimators=2000, n_jobs=-1)\n",
    "rfr_pipe = Pipeline([('scaler', scaler), ('impute', imp), ('smt', smt), ('rfr', rfr)])\n",
    "asa_predictions = cross_val_predict(rfr_pipe, df.iloc[:,0:], y, cv=group_kfold.split(df.iloc[:,0:], y, group_by_admsn), n_jobs=1)\n",
    "#xgb = XGBClassifier(max_depth=5, objective=\"binary:logistic\", n_estimators=2000, n_jobs=-1)\n",
    "#xgb_pipe = Pipeline([('smt', smt), ('scaler', scaler), ('rfr', xgb)])\n",
    "#asa_predictions = cross_val_predict(xgb_pipe, transformed_X, y, cv=cv, n_jobs=5, method='predict_proba')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print \"Number of predictions:\", len(asa_predictions)\n",
    "print df.shape\n",
    "print asa_predictions[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print y[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_death_yn[0:10]\n",
    "death_indices = [i for i,x in enumerate(input_death_yn) if x]\n",
    "live_indices = [i for i,x in enumerate(input_death_yn) if not x]\n",
    "print death_indices[0:10]\n",
    "print len(death_indices)\n",
    "print len(live_indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print asa_predictions[death_indices]\n",
    "print y[death_indices]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(y[death_indices], asa_predictions[death_indices])\n",
    "plt.xlabel(\"True ASA_STATUS\")\n",
    "plt.ylabel(\"Predicted ASA_STATUS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "print mean_squared_error(y[death_indices], asa_predictions[death_indices])\n",
    "print r2_score(y[death_indices], asa_predictions[death_indices])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(y[live_indices], asa_predictions[live_indices])\n",
    "plt.xlabel(\"True ASA_STATUS\")\n",
    "plt.ylabel(\"Predicted ASA_STATUS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "print mean_squared_error(y[live_indices], asa_predictions[live_indices])\n",
    "print r2_score(y[live_indices], asa_predictions[live_indices])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ddd = {\"ADMSN_ID\": admsn_ids, \"OR_CASE_ID\": or_case_id_number, \"PRED_ASA_STATUS\": asa_predictions, \"ASA_STATUS\": y}\n",
    "#ddd = {\"ADMSN_ID\": admsn_ids, \"PRED_DEATH\": asa_predictions[:,1], \"INPT_DEATH_YN\": y}\n",
    "#ddd = {\"ADMSN_ID\": admsn_ids, \"PRED_AKIN_EVENT\": model_probs_cv[2][:,1], \"AKIN_EVENT\": y}\n",
    "#ddd = {\"ADMSN_ID\": admsn_ids, \"PRED_DEATH\": model_probs_cv[2][:,1], \"INPT_DEATH_YN\": y}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print len(y)\n",
    "print len(admsn_ids)\n",
    "#print len(model_probs_cv[2][:,1])\n",
    "new_df = pd.DataFrame(data=ddd)\n",
    "print df.shape\n",
    "\n",
    "#print len(asa_predictions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_df.boxplot(column=\"PRED_ASA_STATUS\", by=\"ASA_STATUS\")\n",
    "plt.title(\"\")\n",
    "plt.suptitle(\"\")\n",
    "plt.ylabel(\"PRED_ASA_STATUS\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_df.head(30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_df.to_csv(os.path.join(dir_to_save_files, \"ASA_predictions_noTime.txt\"), sep=\"\\t\", header=True, index=False)\n",
    "if False:\n",
    "    #pass\n",
    "    new_df.to_csv(\"AKIN_EVENT_predictions.txt\", sep=\"\\t\", header=True, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_df.describe(include=\"all\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:blhill-preop3] *",
   "language": "python",
   "name": "conda-env-blhill-preop3-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}