Preop_Mortality_Test.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",
    "\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",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as st\n",
    "import hashlib\n",
    "\n",
    "from sklearn import metrics, tree\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.linear_model import LogisticRegression, ElasticNet, ElasticNetCV, ARDRegression, SGDClassifier\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.base import BaseEstimator, TransformerMixin\n",
    "from sklearn.externals import joblib\n",
    "\n",
    "#from fancyimpute import SoftImpute\n",
    "from ehr_utils import *\n",
    "#from xgboost import XGBClassifier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Configurations\n",
    "- main_filtered_f - this is the \"main\" data file, which we assume contains all features needed for prediction\n",
    "\n",
    "- test_or_case_id_f - this file contains all the OR_CASE_IDs that we want to predict the mortality risk for\n",
    "\n",
    "- exp_prefix - this corresponds to one of the directories in \"paper\" and this determines which set of features to use to predict mortality. \n",
    "\n",
    "- dir_to_save_files - this should usually be the full path to the directory that contains EHR_MAIN_FEATURES.csv and this is where all output results will be saved\n",
    "\n",
    "- data_dir - this directory contains a file per each categorical variable, and each file contains the allowable values for each variable\n",
    "\n",
    "- model_dir - this is the directory that contains the binary ML models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set to False to hide IDs    \n",
    "verbose = False\n",
    "\n",
    "# version of scikit-learn that was used to create model (should be in pickled model filename)\n",
    "sk_version = \"0.21.3\"\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",
    "\n",
    "MIN_ASA_STATUS=1\n",
    "MAX_ASA_STATUS=5\n",
    "MIN_AGE=18\n",
    "MAX_AGE=89\n",
    "\n",
    "# PATH SETTINGS\n",
    "#main_filtered_f = \"main_Nov21_2017_Feb_13_2018.filtered.main.txt\"\n",
    "main_filtered_f = \"/opt/data/workingdir/blhill/main_merged_w_akin_spo2.filtered.main.txt\"\n",
    "#main_filtered_f = \"/opt/genomics/workingdir/blhill/test_main_3.txt\"\n",
    "test_or_case_id_f = \"/opt/data/workingdir/blhill/or_case_id_032018.txt\"\n",
    "#main_filtered_f = \"/opt/genomics/workingdir/blhill/vali_test.txt\"\n",
    "\n",
    "# directory containing code repo\n",
    "repo_dir = \"/opt/data/workingdir/blhill/code/github/PreopMortalityPrediction\"\n",
    "\n",
    "# experiment prefix sets the set of features to use in the model\n",
    "exp_prefix = \"preop_no_lab_times\"\n",
    "dir_to_save_files = os.path.join(repo_dir, \"paper\", exp_prefix)\n",
    "\n",
    "# directory containing info about acceptable data\n",
    "data_dir = os.path.join(repo_dir, \"data\")\n",
    "\n",
    "## set path to directory containing pickled models\n",
    "#model_dir = os.path.join(\"/opt/genomics/workingdir/blhill/mortality_models\", exp_prefix)\n",
    "model_dir = dir_to_save_files\n",
    "\n",
    "if not os.path.exists(dir_to_save_files):\n",
    "    os.makedirs(dir_to_save_files)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read in data to data frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(os.path.join(home_dir, main_filtered_f), sep=\"|\", header=0)\n",
    "print(df.shape)\n",
    "print(len(df.columns))\n",
    "if verbose:\n",
    "    display(df.iloc[0:20, :])\n",
    "    \n",
    "# make sure INPT_DEATH_YN is set to boolean\n",
    "df[TARGET_VARIABLE] = df[TARGET_VARIABLE].astype(bool)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# hash OR_CASE_ID values and take only patients after March 2018"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_sha256_hash(x):\n",
    "    m = hashlib.sha256()\n",
    "    m.update(str(x).encode('utf-8'))\n",
    "    return m.hexdigest().upper()\n",
    "\n",
    "df[\"OR_CASE_ID\"] = df[\"OR_CASE_ID\"].apply(get_sha256_hash)"
   ]
  },
  {
   "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": [
    "# # get list of surgeries for testing\n",
    "test_or_case_ids = pd.read_csv(test_or_case_id_f, header=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[df[\"OR_CASE_ID\"].isin(test_or_case_ids.iloc[:,0])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# drop any rows that are exact copies of another row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df.shape)\n",
    "df.drop_duplicates(inplace=True)\n",
    "print(df.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# read in features file"
   ]
  },
  {
   "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": [
    "df[contin_vars] = df[contin_vars].astype(float)\n",
    "df[cat_vars] = df[cat_vars].astype(object)\n",
    "df[\"ASA_STATUS\"] = df[\"ASA_STATUS\"].astype(float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# case_list = [553325]\n",
    "# df = df[df.OR_CASE_ID.isin(case_list)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.options.display.max_columns = 1000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[final_features].isna().sum().sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[final_features]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df.apply(lambda x: get_sha256_hash(tuple(x[final_features])), axis = 1)"
   ]
  },
  {
   "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",
    "except ValueError:\n",
    "    pass\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": [
    "# Get acceptable values for categorical variables and filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_acceptable_vals(filename):\n",
    "    with open(os.path.join(data_dir, filename)) as f:\n",
    "        return [l.strip() for l in f.readlines()]\n",
    "\n",
    "pre_surg_location_vals = read_acceptable_vals(\"PRE_SURG_LOCATION_unique_values.txt\")\n",
    "pat_class_vals = read_acceptable_vals(\"PAT_CLASS_unique_values.txt\")\n",
    "hcup_code_vals = read_acceptable_vals(\"HCUP_CODE_unique_values.txt\")\n",
    "gender_vals = read_acceptable_vals(\"GENDER_unique_values.txt\")\n",
    "case_srv_name_vals = read_acceptable_vals(\"CASE_SRV_NAME_unique_values.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(case_srv_name_vals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[df.PRE_SURG_LOCATION.isin(pre_surg_location_vals)]\n",
    "print(df.shape)\n",
    "df = df[df.PAT_CLASS.isin(pat_class_vals)]\n",
    "print(df.shape)\n",
    "df = df[df.HCUP_CODE.astype(float).isin(hcup_code_vals)]\n",
    "print(df.shape)\n",
    "df = df[df.GENDER.isin(gender_vals)]\n",
    "print(df.shape)\n",
    "df = df[df.CASE_SRV_NAME.isin(case_srv_name_vals)]\n",
    "print(df.shape)"
   ]
  },
  {
   "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'])]\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": [
    "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": {},
   "outputs": [],
   "source": [
    "# print \"Number of Patients:\", df.shape[0]\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 \"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 \"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",
    "# if verbose:\n",
    "#     print(df[\"CASE_SRV_NAME\"].value_counts())\n",
    "#     print(df[\"CASE_SRV_NAME\"].value_counts()/df.shape[0])*100\n",
    "#     print(\"=\"*40)\n",
    "#     print(df[df[\"INPT_DEATH_YN\"] == 1][\"CASE_SRV_NAME\"].value_counts())\n",
    "#     print(df[df[\"INPT_DEATH_YN\"] == 1][\"CASE_SRV_NAME\"].value_counts()/df[\"INPT_DEATH_YN\"].sum())*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Remove outlier values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# string_cols = ['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']\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",
    "# dff = df[contin_vars]\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)"
   ]
  },
  {
   "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 "
   ]
  },
  {
   "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']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[final_features + [TARGET_VARIABLE]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# One-hot encode categorical variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for var in cat_vars:\n",
    "    try:\n",
    "        # drop_first uses k-1 dummies out of k categories\n",
    "        print(var)\n",
    "        #df = pd.get_dummies(df, columns=[var], drop_first=True)\n",
    "        df = pd.get_dummies(df, columns=[var])\n",
    "        pass\n",
    "    except ValueError:\n",
    "        pass\n",
    "    except KeyError:\n",
    "        print(var, 'already dropped')\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": [
    "# if len(feature_whitelist) == 0:\n",
    "for cat in cat_to_drop:\n",
    "    try:\n",
    "        df.drop(cat, axis=1, inplace=True)\n",
    "    except KeyError:\n",
    "        print(cat, 'already dropped')\n",
    "#print df.columns.values\n",
    "for col in sorted(df.columns.values):\n",
    "    print(col, \"\\t\\t\", df[col].dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Remove target variables from data frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df[TARGET_VARIABLE].mean())\n",
    "print(\"Column names:\", df.columns.values)\n",
    "try:\n",
    "    y = np.ravel(df[TARGET_VARIABLE])\n",
    "    #asa_status = df[\"ASA_STATUS\"]\n",
    "    df.drop(TARGET_VARIABLE, axis=1, inplace=True, errors='ignore')\n",
    "    df.drop(outcome_vars, axis=1, inplace=True, errors='ignore')\n",
    "    input_death_yn = df['INPT_DEATH_YN']\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())\n",
    "print(y.sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make sure we have all features that we used to train model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_feature_list = pd.read_csv(os.path.join(dir_to_save_files, exp_prefix + \"_final_feature_list.txt\"), header=None)\n",
    "print(final_feature_list.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in final_feature_list[0]:\n",
    "    if col not in df.columns.values:\n",
    "        print(col, \"missing from dataframe\")\n",
    "        df[col] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(final_feature_list[0].shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.columns.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df[final_feature_list[0]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Standardize training, testing data"
   ]
  },
  {
   "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 = StandardizeWithNaN()\n",
    "# scaler.fit(X_test)\n",
    "# X_test = scaler.transform(X_test)\n",
    "#scaler.fit(df)\n",
    "scaler = pickle.load(open(os.path.join(model_dir, \"StandardizeWithNaN.pkl\"), \"rb\"))\n",
    "X_test = scaler.transform(df)\n",
    "y_test = y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Impute missing values"
   ]
  },
  {
   "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": [
    "print(\"imputing X_test\")\n",
    "#print(np.isnan(X_test).any())\n",
    "#si = SoftImputeEstimator()\n",
    "si = pickle.load(open(os.path.join(model_dir, \"MeanImputer.pkl\"), \"rb\"))\n",
    "si.statistics_[np.isnan(si.statistics_)] = 0.\n",
    "print(si.statistics_)\n",
    "X_test = si.transform(X_test.replace(np.inf, np.nan))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load model(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# models = {}\n",
    "models = []\n",
    "for i in range(1):\n",
    "    model_file_name = \"Random Forest_train_sk{}.pkl\".format(sk_version)\n",
    "    model = pickle.load(open(os.path.join(model_dir, model_file_name), \"rb\"), encoding='latin1')\n",
    "    print(\"Loaded\", model_file_name)\n",
    "    #models[\"fold_{}\".format(i)] = model\n",
    "    models.append(model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "models[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Predict classes and get probability of labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# model_predictions = {k: model.predict(np.array(X_test)) for k, model in models.items()}\n",
    "# model_probs = {k: model.predict_proba(X_test) for k, model in models.items()}\n",
    "\n",
    "model_predictions = [model.predict(np.array(X_test)) for model in models]\n",
    "model_probs = [model.predict_proba(X_test) for model in models]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# model_probs_df = pd.DataFrame.from_dict({k: probs[:,1] for k, probs in model_probs.items()})\n",
    "# model_probs_df.std(axis=1).plot()\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# If we have true labels, see how well the model did "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_names = [\"Random Forest\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot_cross_val_roc_curve(model_probs, y_test, \"test_cross_val_roc.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot_accuracy_roc_auc(models, model_names, model_predictions, model_probs, y_train, y_test, \"test_accuracy_roc_auc.tif\")"
   ]
  },
  {
   "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,\"test_roc_curve.png\"))"
   ]
  },
  {
   "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, \"test_precision_recall_curve.png\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(model_probs[0])\n",
    "len(or_case_id_number)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if verbose:\n",
    "    for i in zip(model_probs[0][:,1], or_case_id_number, y):\n",
    "        print(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plt.hist(model_probs[0][y == False,1], bins=20)\n",
    "plt.hist(model_probs[0][y == True,1], bins=20)\n",
    "plt.xlabel(\"Probability\")\n",
    "plt.ylabel(\"Patient Count\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(model_probs[0][y == False,1].mean())\n",
    "print(model_probs[0][y == False,1].std())\n",
    "print(model_probs[0][y == True,1].mean())\n",
    "print(model_probs[0][y == True,1].std())"
   ]
  },
  {
   "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
}