mimic-iv-demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "055e570d",
   "metadata": {},
   "source": [
    "# Introduction to the MIMIC Database and Machine Learning with scikit-learn\n",
    "\n",
    "**Author: Meghan R. Hutch**\n",
    "\n",
    "**Date: Februrary 8th, 2024**\n",
    "\n",
    "Inspired and adapted from Dr. Garrett Eickelberg's workshop on MIMIC-III and scikit learn: https://github.com/geickelb/MIMIC-III_to_Model/tree/master\n",
    "\n",
    "---\n",
    "\n",
    "### **Today's lecture will be divided into two parts:**\n",
    "\n",
    "* Introduction to working with the MIMIC-IV dataset\n",
    "<br></br>\n",
    "* Introduction to machine learning with scikit learn\n",
    "\n",
    "We will discuss both of these in the context of a framework for conducting clinical research"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "084c22f1",
   "metadata": {},
   "source": [
    "# Framework for Conducting Clinical Research \n",
    "\n",
    "---\n",
    "\n",
    "1. Research Question Specification (hypothesis)\n",
    "<br></br>\n",
    "2. Cohort Specification\n",
    "<br></br>\n",
    "3. Data Extraction\n",
    "<br></br>\n",
    "4. Data Pre-processing (cleaning!)\n",
    "<br></br>\n",
    "5. Model Training\n",
    "<br></br>\n",
    "5. Model Evaluation\n",
    "<br></br>\n",
    "\n",
    "**Today, we will largely focus on steps 3-6**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1529e06d",
   "metadata": {},
   "source": [
    "# Introduction to MIMIC-IV\n",
    "\n",
    "MIMIC-IV dataset is the most recently updated version of the MIMIC dataset. MIMIC curates the data of patients admitted to BIDMC emergency department or an ICU between 2008-2019. The dataset has been de-identified and contains the following (as of the MIMIC-IV v2.2 release):\n",
    "\n",
    "* 299,712 patients\n",
    "* 431,231 admissions\n",
    "* 73,181 icustays\n",
    "\n",
    "Useful [documentation](https://mimic.mit.edu/docs/iv/) for using the MIMIC-IV dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e9e4be5",
   "metadata": {},
   "source": [
    "## Demo-Dataset\n",
    "\n",
    "Today we will use the previously curated [demo-dataset](https://physionet.org/content/mimic-iv-demo/2.2/)\n",
    "\n",
    "The demo-dataset contains the data from a random subset of 100 hospitalized patients. This dataset was curated for the purposes of workshops and tutorials. Thus, having appropriate CITI training is not needed. The full dataset is available on the class quest server, but a physionet account and [CITI training](https://physionet.org/about/citi-course/) should be completed before accessing it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23991120",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import our packages\n",
    "import os\n",
    "\n",
    "import pandas as pd \n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abc9914d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set our working directory - this notebook assume you are in the folder where the MIMIC-IV demo dataset is located\n",
    "os.chdir('/projects/e30766/mimic-iv-demo/2.2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edb5b820",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make sure you are in the working directory: '/projects/e30766/mimic-iv-demo/2.2'\n",
    "os.getcwd()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ef1578d",
   "metadata": {},
   "source": [
    "# Data Extraction\n",
    "\n",
    "- The demo dataset has been downloaded to the classroom's quest allocation.\n",
    "- We will further load this data into our notebook\n",
    "\n",
    "### Load in MIMIC data\n",
    "\n",
    "First, we will load the list of subject identifiers (`subject_id`) for the demo cohort"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5064c9af",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data\n",
    "subject_id = pd.read_csv('demo_subject_id.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "254b15fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# view the first 5 rows of data\n",
    "subject_id.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c0131df",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check how many unique subjects\n",
    "subject_id.nunique()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0941b1f",
   "metadata": {},
   "source": [
    "## MIMIC-IV Tables\n",
    "\n",
    "For this demo, we will explore data found in the following two MIMIC-IV modules:\n",
    "\n",
    "\n",
    "### [Hosp](https://mimic.mit.edu/docs/iv/modules/hosp/)\n",
    "\n",
    "> \"The Hosp module provides all data acquired from the hospital wide electronic health record. Information covered includes patient and admission information, laboratory measurements, microbiology, medication administration, and billed diagnoses.\"\n",
    "\n",
    "### [ICU](https://mimic.mit.edu/docs/iv/modules/icu/)\n",
    "\n",
    "> \"The ICU module contains information collected from the clinical information system used within the ICU. Documented data includes intravenous administrations, ventilator settings, and other charted items.\"\n",
    "\n",
    "As part of the 100 patient demo dataset, we have also downloaded all of the accompanying tables within the `hosp` and `icu` modules. Data/csv files can be found: `/projects/e30766/mimic-iv-demo/2.2`\n",
    "\n",
    "MIMIC-IV also contains additional modules for data from emergency department (ED), chest X-rays (CXR), and clinical notes (Note)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22091455",
   "metadata": {},
   "source": [
    "### Load in the `patient` table to aquire demographic data and admission date information\n",
    "\n",
    "Important note about MIMIC's pre-processing of [dates](https://mimic.mit.edu/docs/iv/modules/hosp/patients/#anchor_age-anchor_year-anchor_year_group):\n",
    "\n",
    "* Dates are randomly shifted consistently for all patients\n",
    "\n",
    "* `anchor_year` - shifted year for the patient\n",
    "\n",
    "* `anchor_year_group` -  is a range of years - the patient’s anchor_year occurred during this range\n",
    "\n",
    "* `anchor_age` - the patient’s age in the `anchor_year` (year of admission). If a patient’s `anchor_age` is over 89 in the `anchor_year` then their `anchor_age` is set to 91, regardless of how old they actually were.\n",
    "\n",
    "> Example: a patient has an `anchor_year` of 2153, `anchor_year_group` of 2008 - 2010, and an `anchor_age` of 60.\n",
    "The year 2153 for the patient corresponds to 2008, 2009, or 2010.\n",
    "The patient was 60 in the shifted year of 2153, i.e. they were 60 in 2008, 2009, or 2010.\n",
    "A patient admission in 2154 will occur in 2009-2011, an admission in 2155 will occur in 2010-2012, and so on.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e68d33d",
   "metadata": {},
   "outputs": [],
   "source": [
    "patients = pd.read_csv('hosp/patients.csv.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab710a62",
   "metadata": {},
   "outputs": [],
   "source": [
    "# review the first 5 rows\n",
    "patients.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92cb7891",
   "metadata": {},
   "outputs": [],
   "source": [
    "# review the number of unique values in each dataset\n",
    "patients.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc78ad6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(patients)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eeaf294e",
   "metadata": {},
   "source": [
    "### Load data files into python\n",
    "\n",
    "Here, we will load a few of the tables we will use to demonstate the workflow of aquiring, pre-processing, and analyzing data from the MIMIC database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f4a6a73",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load several tables from the hosp module\n",
    "omr = pd.read_csv('hosp/omr.csv.gz')\n",
    "labevents = pd.read_csv('hosp/labevents.csv.gz')\n",
    "d_labitems = pd.read_csv('hosp/d_labitems.csv.gz')\n",
    "\n",
    "# load several tables from the icu module\n",
    "chartevents = pd.read_csv('icu/chartevents.csv.gz')\n",
    "d_items = pd.read_csv('icu/d_items.csv.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10c13328",
   "metadata": {},
   "source": [
    "### omr table\n",
    "\n",
    "> \"The Online Medical Record (OMR) table contains miscellaneous information from the EHR.\"\n",
    "\n",
    "[omr documentation](https://mimic.mit.edu/docs/iv/modules/hosp/omr) of columns:\n",
    "\n",
    "* `chartdate` - the date on which the observation was recorded\n",
    "* `seq_num`- a monotonically increasing integer which uniquely distinguishes results of the same type recorded on the same day. For example, if two blood pressure measurements occur on the same day, seq_num orders them chronologically.\n",
    "* `result_name` - human interpretable description of the observation\n",
    "* `result_value` - the avalue associated with the observation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f9b6c6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# review the first 5 rows\n",
    "omr.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6e5dc8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# review unique result_names (observations)\n",
    "# here we can see the unique variables/measurements contained in this table\n",
    "omr[['result_name']].drop_duplicates()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bab7d02",
   "metadata": {},
   "outputs": [],
   "source": [
    "omr.isnull().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c73f4330",
   "metadata": {},
   "source": [
    "# Data Pre-processing (Data Cleaning)\n",
    "\n",
    "---\n",
    "\n",
    "### What does it mean to clean data?\n",
    "\n",
    "Simply, data cleaning is the act of preparing data for analysis. In a clinical context, we also want to make sure that the data we are using is clinically relevant and accurate, especially in regard to the specific clinical question or problem we are trying to solve.\n",
    "\n",
    "### How do we \"clean\" data:\n",
    "\n",
    "* Ensure data are formatted into the right data types (numeric, character, date, etc)\n",
    "\n",
    "* Review summary statistics (how many observations, how many patients have each observation, mean/median/sd/variance of values, proportion of missing data)\n",
    "\n",
    "* Review the distribution of data\n",
    "\n",
    "* Are there possible data entry problems or unit conversions needed?\n",
    "\n",
    "* How should we handle missing data?\n",
    "\n",
    "* How many observations does each patient have? How should we handle repeated observations?\n",
    "\n",
    "---\n",
    "\n",
    "“Data science, he says, involves multiple “very small decisions” — data cleaning and filtering steps, for instance, which are crucially important, but difficult to document. And journal page limits preclude exposition. But by blending code, data and text in a single document, researchers can show just how their results were generated.” - Ben Marwick; [Perkel, JM](https://www.nature.com/articles/d41586-022-00563-z)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "766a63a1",
   "metadata": {},
   "source": [
    "### Let's focus on the variables `Height` and `Weight`\n",
    "\n",
    "We will filter the `omr` table to only contains results for height and weight. During data cleaning, each variable will be saved in its own dataframe for ease of pre-processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cfb1122",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create separate dataframes for each variable\n",
    "height = omr[omr['result_name']==\"Height (Inches)\"]\n",
    "\n",
    "weight = omr[omr['result_name']==\"Weight (Lbs)\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6db2b8ac",
   "metadata": {},
   "source": [
    "### Height (Inches)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d578015b",
   "metadata": {},
   "outputs": [],
   "source": [
    "height.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7f7bae6",
   "metadata": {},
   "source": [
    "### Review variable types \n",
    "\n",
    "It looks like `result_value` is a continuous variable. Let's make sure that Python also has recognized it as such."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98a6482b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check variable data types\n",
    "height.dtypes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d9b1b48",
   "metadata": {},
   "source": [
    "It looks like Python read `result_value` as a object (or character) variable. Let's convert to numeric. This is an important step to make sure any downstream functions (calculating the mean or median, plotting the distribution, etc) work correctly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80ac6731",
   "metadata": {},
   "outputs": [],
   "source": [
    "height['result_value'] = pd.to_numeric(height['result_value'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f73725a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check variable data types again to make sure our transformation and code worked\n",
    "# float64 indicates a continuous variable, so this looks good!\n",
    "height.dtypes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20d9236e",
   "metadata": {},
   "source": [
    "#### Summary statistics\n",
    "\n",
    "Next, let's calculate summary statistics to get a better sense of our variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8c506f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# how many unique patients have a height recording\n",
    "# 61 patients - not all of our 100 patients had a record for height\n",
    "height.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21a35b72",
   "metadata": {},
   "outputs": [],
   "source": [
    "# how many observations - 378 observations, thus some patients must have height recorded more than once\n",
    "len(height)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9a8974f",
   "metadata": {},
   "source": [
    "### Visualize distribution of values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dee820d",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.distplot(height['result_value'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c74ab3a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# use describe() to compute summary statistics of the distribution\n",
    "height['result_value'].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bccaf75",
   "metadata": {},
   "source": [
    "**Potential red flag:** Is there really a patient who might be 5 inches tall? This is a clear case I'd like to investigate more. Perhaps this is a data entry issue or a newborn?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fba4840",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sort values from lowest to highest - I can then retrieve the subject_id of the patient with a height of 5 inches\n",
    "# I can also see if there are any other patients with unusually low heights\n",
    "height.sort_values(['result_value'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7fd8904",
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's see if this patient has other recorded height measurements\n",
    "height[height['subject_id']==10012853]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "059cd240",
   "metadata": {},
   "source": [
    "### Investigate patients with multiple observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a258a3aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for each unique patient, we will calculate the standard deviation of the patient's height measurements\n",
    "# ddof = 0 indicates that std for patients with one measurement will be displayed as 0 rather than NaN\n",
    "height['std'] = height.groupby('subject_id')['result_value'].transform(np.std, ddof=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53e4b0d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "height.sort_values('std', ascending = False)[['subject_id', 'std']].drop_duplicates().head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fe650c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "height[height['subject_id']==10005909]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f48057fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "height[height['subject_id']==10019917]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5634f4ca",
   "metadata": {},
   "source": [
    "### For patients with multiple observations, let's calculate the median\n",
    "\n",
    "**Notes:** \n",
    "\n",
    "* The best way to handle these types of inconsistencies will largely depend on the exact variable, question of interest, and domain knowledge. \n",
    "<br> </br>\n",
    "* Seek consultation from mentors and domain experts! \n",
    "<br> </br>\n",
    "* In this case, I'm not considering the date or specific admission. Some patients have multiple admissions. \n",
    "<br> </br>\n",
    "* Additionally, depending how old they are, height may be expected to change over time. In those cases, we might want to consider the patient's age before just taking the median"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfb0995f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a new variable with the median height value for each patient\n",
    "height['result_value_median'] = height.groupby('subject_id')['result_value'].transform('median')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a98617ef",
   "metadata": {},
   "source": [
    "**Save the new dataset**\n",
    "\n",
    "This next step ensures that we will now only have one unique row per patient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ba5ff40",
   "metadata": {},
   "outputs": [],
   "source": [
    "# assuming we don't need to link by date\n",
    "height_clean = height[['subject_id', 'result_name', 'result_value_median']].drop_duplicates()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9331f22f",
   "metadata": {},
   "source": [
    "Let's review the distribution again - now it looks a bit more normally distributed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7d72d93",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.distplot(height_clean['result_value_median'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dbdb5d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "height_clean.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cab04a1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(height_clean)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "792882fd",
   "metadata": {},
   "source": [
    "### Weight (lbs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdbdf0a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "weight.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f31c3d57",
   "metadata": {},
   "source": [
    "### Review variable types \n",
    "\n",
    "It looks like `result_value` is a continuous variable. Let's make sure that Python also has recognized it as such."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39db656c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check variable data types\n",
    "weight.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63c76423",
   "metadata": {},
   "outputs": [],
   "source": [
    "weight['result_value'] = pd.to_numeric(weight['result_value'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33eb7d22",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check variable data types\n",
    "weight.dtypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed8487aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "weight['result_value'].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2094cd6b",
   "metadata": {},
   "source": [
    "### Distribution of values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a68a72ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.distplot(weight['result_value'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdbaf178",
   "metadata": {},
   "outputs": [],
   "source": [
    "weight.sort_values(['result_value'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb5e64d7",
   "metadata": {},
   "source": [
    "### Investigate the variance of in measurements between patients with multiple observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9b44176",
   "metadata": {},
   "outputs": [],
   "source": [
    "weight['std'] = weight.groupby('subject_id')['result_value'].transform(np.std, ddof=0)# ddof = 0 indicates that std for patients with one measurement will be displayed as 0 rather than NaN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5ef71ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "weight.sort_values('std', ascending = False)[['subject_id', 'std']].drop_duplicates().head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7493dee",
   "metadata": {},
   "outputs": [],
   "source": [
    "weight[weight['subject_id']==10019385] # did this patient lose > 110 pounds in < 1 month?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b53af5d",
   "metadata": {},
   "source": [
    "**Potential Unit Issue? 97 kg ~ 214 pounds**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f9e92a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "weight[weight['subject_id']==10021487].sort_values('chartdate')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f7ea010",
   "metadata": {},
   "source": [
    "### Notes:\n",
    "\n",
    "In the first case, my intuition is that the discrepancy in weight measurements is due to a unit conversions issue: 97 kg ~ 214 pounds\n",
    "\n",
    "In the second case, the patient has many weight measurements over the a 1.5 year period. They are not varrying too widely. It looks like there is an additional weight loss, followed by an increase in weight (U-shaped curve)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86b7e80c",
   "metadata": {},
   "source": [
    "### Notes on Missing Data\n",
    "\n",
    "This could be it's own class. \n",
    "\n",
    "Rule of thumb: Talk to your domain experts/collaborators!\n",
    "\n",
    "Some common ways to handle missing data:\n",
    "\n",
    "* Remove those patients\n",
    "* Impute mean or median\n",
    "* Forward/backward filling\n",
    "* MICE\n",
    "\n",
    "In the above cases, I used the median to impute. This may or may not be the best method depending on your use cases. Here, I wanted to keep it simple!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79ecd716",
   "metadata": {},
   "source": [
    "### Evaluating Unit Conversions\n",
    "\n",
    "Another clinical variable that often needs cleaning are laboratory test measurements.\n",
    "\n",
    "In MIMIC we can use the `labevents` and `d_labitems` tables to work with laboratory data. \n",
    "\n",
    "Importantly, `d_labitems` links to the `labevents` table. \n",
    "\n",
    "The `d_labitems` provides the definitions (or names) or the specific labs via the `itemid`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56cbe6cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's first look at the labevents table\n",
    "labevents.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc8cdc50",
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's first look at the d_labitems table\n",
    "d_labitems.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62d5c5e5",
   "metadata": {},
   "source": [
    "### There are a few ways you could think about working with the labs dataset. \n",
    "\n",
    "1. We can merge `d_labitems` and `labevents` by shared identifiers (we'll review this in the next section on merging tables)\n",
    "\n",
    "2. We can search for specific labs of interest via the `label` column\n",
    "\n",
    "3. We can also filter with specific `itemid` numbers if we know which correspond to labs of interest"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9be92774",
   "metadata": {},
   "source": [
    "First, let's try and see whether we can find the `itemid` for potassium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "061fbc9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "d_labitems[d_labitems['label'].str.contains('potassium', case = False, na=False)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c15eb0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "potassium = labevents[labevents['itemid']==50971]\n",
    "potassium.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1eb659c",
   "metadata": {},
   "outputs": [],
   "source": [
    "potassium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0771e36a",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.distplot(potassium['valuenum'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "199c0143",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# let's check the distinct units of potassium\n",
    "potassium[['valueuom']].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e236f58",
   "metadata": {},
   "source": [
    "### Let's evaluate another lab test, one which has more than one unit specification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "318115e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "d_labitems[d_labitems['itemid']==51249]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e57a9eba",
   "metadata": {},
   "outputs": [],
   "source": [
    "MCHC = labevents[(labevents['itemid']==51249)]\n",
    "MCHC['valueuom'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d1a4a9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "lab_perc = MCHC[MCHC['valueuom']=='%']\n",
    "lab_perc.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e215fdf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "lab_gdl = MCHC[MCHC['valueuom']=='g/dL']\n",
    "lab_gdl.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eccde471",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.distplot(lab_gdl['valuenum'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e8b71aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.distplot(lab_perc['valuenum'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cea31750",
   "metadata": {},
   "source": [
    "## Merging tables\n",
    "\n",
    "Tables can be merged (e.g. linked) together through shared identifiers. The MIMIC official [documentation](https://mimic.mit.edu/docs/iv/) is a great reference for learning more about the associations between tables\n",
    "\n",
    "---\n",
    "\n",
    "### First, we will merge the different lab tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "322d779d",
   "metadata": {},
   "outputs": [],
   "source": [
    "labevents_with_itemid = pd.merge(d_labitems, labevents,\n",
    "                      on = ['itemid'],\n",
    "                      how = 'inner')\n",
    "\n",
    "labevents_with_itemid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16fd5ec5",
   "metadata": {},
   "source": [
    "### Let's merge some additional MIMIC-IV tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a9a871b",
   "metadata": {},
   "outputs": [],
   "source": [
    "admissions = pd.read_csv('hosp/admissions.csv.gz')\n",
    "diagnoses = pd.read_csv('hosp/diagnoses_icd.csv.gz')\n",
    "d_icd_diagnoses = pd.read_csv('hosp/d_icd_diagnoses.csv.gz')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47aafd49",
   "metadata": {},
   "source": [
    "### Merge diagnoses (icd) tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40220c76",
   "metadata": {},
   "outputs": [],
   "source": [
    "diagnoses.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b051443",
   "metadata": {},
   "outputs": [],
   "source": [
    "d_icd_diagnoses.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83cf0657",
   "metadata": {},
   "outputs": [],
   "source": [
    "diagnoses_with_code = pd.merge(diagnoses, \n",
    "                         d_icd_diagnoses,\n",
    "                         on = ['icd_code', 'icd_version'],\n",
    "                         how = 'inner')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70035a27",
   "metadata": {},
   "outputs": [],
   "source": [
    "diagnoses_with_code"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e2c2a30",
   "metadata": {},
   "source": [
    "### Let's merge the new diagnoses table with the admissions table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a068f1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "admissions.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39a6ecaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "admissions_with_icd = pd.merge(admissions,\n",
    "                               diagnoses_with_code,\n",
    "                               on = ['subject_id', 'hadm_id'],\n",
    "                               how = 'inner')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1dc86c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "admissions_with_icd.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "659b6acc",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "admissions_with_icd.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd1bc2a3",
   "metadata": {},
   "source": [
    "### Filtering datasets (data cleaning / cohort specification)\n",
    "\n",
    "Let's say we are only interested in patients with sepsis. We could filter our dataset using keyword matching or specific ICD codes (if we have specific codes we are interested in)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2417bb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "sepsis_cohort = admissions_with_icd[admissions_with_icd['long_title'].str.contains('sepsis', case = False)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff80e49c",
   "metadata": {},
   "outputs": [],
   "source": [
    "sepsis_cohort.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b2d01b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "sepsis_cohort.nunique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f38d441a",
   "metadata": {},
   "outputs": [],
   "source": [
    "sepsis_cohort[['icd_code', 'icd_version', 'long_title']].drop_duplicates()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83b7016e",
   "metadata": {},
   "outputs": [],
   "source": [
    "sepsis_cohort_icd = admissions_with_icd[admissions_with_icd['icd_code']=='99592']\n",
    "sepsis_cohort_icd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "462dba09",
   "metadata": {},
   "source": [
    "# Introduction to Machine Learning with scikit-learn\n",
    "\n",
    "---\n",
    "\n",
    "For this demonstration, we will use a previously prepared dataset from an earlier version of MIMIC. This data comes from a [community challenge](https://physionet.org/content/challenge-2012/1.0.0/) focused on predicting mortality with the MIMIC dataset.\n",
    "\n",
    "To demonstate the use of scikit-learn for machine learning, we will use the previously pre-processed dataset from one of the challenge winners (https://github.com/alistairewj/challenge2012). Data can be accessed via quest or from Alistair Johnson's github. He also includes his [pre-processing notebook](https://github.com/alistairewj/challenge2012/blob/master/prepare-data.ipynb) which can provide another example of how to proceed with cleaning clinical data.\n",
    "\n",
    "### Hypothesis: Machine Learning can facilitate prediction of ICU mortality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fdff6f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load in our machine learning functions from scikit-learn\n",
    "from sklearn.metrics import confusion_matrix, ConfusionMatrixDisplay, classification_report, roc_curve, roc_auc_score, precision_score, recall_score, accuracy_score, auc, average_precision_score, precision_recall_curve, precision_recall_fscore_support, f1_score, log_loss\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import KFold, GridSearchCV,RandomizedSearchCV,cross_validate\n",
    "from sklearn.preprocessing import StandardScaler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76b555f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option('display.max_columns', None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b609a367",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load in pre-processed challenge dataset\n",
    "data_ml = pd.read_csv('PhysionetChallenge2012-set-b.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b8dc3f2",
   "metadata": {},
   "source": [
    "### Let's review our dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "151e032d",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_ml.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca142212",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_ml.nunique()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ad42edf",
   "metadata": {},
   "source": [
    "### **Notes**: \n",
    "\n",
    "* This dataset has 4000 unique subjects and > 100 columns.\n",
    "\n",
    "* As mentioned, this dataset is already pre-processed. I would still advocate for performing your own data quality checks. For the purposes of this exercise, my checks will be extremely limited. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4fff5a9",
   "metadata": {},
   "source": [
    "### Review missing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7927c8fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# count how many variables have no missing values; We have 10 columns with complete entries\n",
    "data_ml.dropna(axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a4333cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# count how many variables have < 100 missing values\n",
    "data_ml.columns[ np.sum(data_ml.isnull()) < 100 ]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bf1e6e9",
   "metadata": {},
   "source": [
    "### Our limited data pre-processing\n",
    "\n",
    "* We will subset the data frame to keep only the columns that have < 100 missing values. \n",
    "* To keep things simple, we will also just take the variables that are the `first` measurement\n",
    "* We will subsequently drop the patients who have missing data.\n",
    "* Note: `In-hospital_death` is our outcome variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b387d0a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we will remove `GCS_last`, 'Length_of_stay', and 'Survival' since these may leak info about our likely deceased patients\n",
    "data_ml_clean = data_ml[['In-hospital_death', 'Age', 'SAPS-I', 'SOFA', 'Weight',\n",
    "       'Height', 'GCS_first', 'Glucose_first', 'HR_first', 'Temp_first',\n",
    "       'BUN_first', 'Creatinine_first', 'HCO3_first', 'HCT_first', 'K_first', 'Mg_first',\n",
    "       'Na_first', 'Platelets_first', 'WBC_first']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b6fc4c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# drop any rows (patients) with a missing value\n",
    "data_ml_clean = data_ml_clean.dropna()\n",
    "data_ml_clean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce04895d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we have ~2,000 rows (or unique patients)\n",
    "len(data_ml_clean)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51af8f17",
   "metadata": {},
   "source": [
    "### Prepare data for modeling\n",
    "\n",
    "For continuous varibales, many machine learning algorithms work best when the data has been normalized (mean = 0, standard deviation 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e396b2aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify the columns you want to standardize (in this case, we want to standardize every column by Gender and our outcome `In-hospital_death`)\n",
    "columns_to_standardize = ['Age', 'SAPS-I', 'SOFA', 'Weight',\n",
    "       'Height', 'GCS_first', 'Glucose_first', 'HR_first', 'Temp_first',\n",
    "       'BUN_first', 'Creatinine_first', 'HCO3_first', 'HCT_first', 'K_first', 'Mg_first',\n",
    "       'Na_first', 'Platelets_first', 'WBC_first']\n",
    "\n",
    "# Create a StandardScaler instance\n",
    "scaler = StandardScaler()\n",
    "\n",
    "# Fit and transform only the specified columns\n",
    "data_ml_clean[columns_to_standardize] = scaler.fit_transform(data_ml_clean[columns_to_standardize])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8dbc88f",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_ml_clean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e83aeb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_ml_clean['Age'].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9661cce6",
   "metadata": {},
   "source": [
    "## Create Training and Test Sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "037bb296",
   "metadata": {},
   "outputs": [],
   "source": [
    "## first, we will create two dataframes. \n",
    "\n",
    "# X will contain only our predictors \n",
    "X = data_ml_clean.drop(['In-hospital_death'], axis = 1)\n",
    "\n",
    "# y contains our outcome \n",
    "y = data_ml_clean['In-hospital_death']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e1d77d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=12345, shuffle = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0c351aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train.value_counts()#/len(y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8f3f6ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_test.value_counts()#/len(y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6163ec83",
   "metadata": {},
   "source": [
    "# Model Training & Evaluation\n",
    "\n",
    "---\n",
    "\n",
    "## Logistic Regression Algorithm\n",
    "\n",
    "<img src=\"lecture_images/logistic_regression.png\" alt=\"title\" width=\"400\" align=\"left\">\n",
    "\n",
    "\n",
    "* Commonly used for classification, thus need a binary outcome variable (e.g. survival or death)\n",
    "* Good baseline model (simple!)\n",
    "* L1 (Lasso) vs L2 (Ridge) - penalty to avoid overfitting\n",
    "\n",
    "[Image Source](https://www.javatpoint.com/logistic-regression-in-machine-learning)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a2cd240",
   "metadata": {},
   "source": [
    "### Train L2 (ridge) logistic regression model using crossfold validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23729000",
   "metadata": {},
   "outputs": [],
   "source": [
    "# instantiate logistic regression model\n",
    "lr = LogisticRegression(penalty='l2',solver='liblinear', random_state = 12345)\n",
    "\n",
    "# use cross fold validation to randomly split data in k folds\n",
    "cv_results = cross_validate(lr, X_train, y_train, cv=10, \n",
    "                            scoring = ['accuracy', 'roc_auc', 'recall', 'precision'],\n",
    "                            return_train_score = False, return_estimator = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5438320b",
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7918b799",
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the mean across folders\n",
    "print('Accuracy:', cv_results['test_accuracy'].mean().round(3))\n",
    "print('AUC:', cv_results['test_roc_auc'].mean().round(3))\n",
    "print('Recall:', cv_results['test_recall'].mean().round(3))\n",
    "print('Precision:', cv_results['test_precision'].mean().round(3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "499b9bd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "##basic model performance - code prepared by Garrett Eickelberg\n",
    "\n",
    "# fit our model on the training set\n",
    "fit_lr = lr.fit(X_train, y_train)\n",
    "\n",
    "# predicted classes using default 0.5 threshold\n",
    "y_hat = fit_lr.predict(X_train) \n",
    "\n",
    "#predicted probabilities\n",
    "y_proba = fit_lr.predict_proba(X_train)[:,1] \n",
    "\n",
    "# auc\n",
    "auc = roc_auc_score(y_train, y_proba)\n",
    "\n",
    "# model loss\n",
    "loss = log_loss(y_train, y_hat)\n",
    "\n",
    "print ('the AUC is: {:0.3f}'.format(auc))\n",
    "print ('the logloss is: {:0.3f}'.format(loss))\n",
    "print(\"classification report:\\n \", classification_report(y_train,y_hat, digits=3))\n",
    "print(\"confusion matrix:\\n \")\n",
    "\n",
    "# Create the confusion matrix\n",
    "cm = confusion_matrix(y_train, y_hat)\n",
    "ConfusionMatrixDisplay(confusion_matrix=cm).plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da245b4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_proba"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e58e0774",
   "metadata": {},
   "source": [
    "## Evaluation of Training \n",
    "\n",
    "---\n",
    "\n",
    "### Receiver Operating Characteristic Curve\n",
    "\n",
    "A ROC plots the true postitve rate (sensitivity) vs false positive rate (1-specificity) at different decision thresholds. This allows us to calculate the area under the curve (AUC) to estimate model performance. An AUC of 0.5 indicates a model that performs as well as chance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acf27208",
   "metadata": {},
   "outputs": [],
   "source": [
    "fpr, tpr, thresholds = roc_curve(y_train, y_proba, pos_label=1)\n",
    "\n",
    "plt.title('ROC curve')\n",
    "ax1 = plt.plot(fpr, tpr, 'b', label = '%s AUC = %0.3f' % ('', auc), linewidth=2)\n",
    "plt.legend(loc = 'lower right')\n",
    "plt.plot([0, 1], [0, 1],'r--')\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([0, 1])\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.xlabel('False Positive Rate')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e616fdcd",
   "metadata": {},
   "source": [
    "### Precision-Recall Curve\n",
    "\n",
    "The curve plots the prcesion and recall at different decision thresholds. Recall or the true positive rate tells us how well the model correctly identifies all positive cases, whereas precision estimates the proportion of cases the model identified as positive that were truely positive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddd96acc",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_proba = fit_lr.predict_proba(X_train)[:,1]\n",
    "\n",
    "precision, recall, thresholds = precision_recall_curve(y_train, y_proba, pos_label=1, sample_weight=None)\n",
    "avg_p = average_precision_score(y_train, y_proba, pos_label=1, sample_weight=None)\n",
    "\n",
    "plt.title('Precision-Recall curve')\n",
    "ax1= plt.plot(precision, recall, 'b', label = '%s AP = %0.3f' % ('', avg_p), linewidth=2)\n",
    "plt.legend(loc = 'lower left')\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([0, 1])\n",
    "plt.ylabel('Precision')\n",
    "plt.xlabel('Recall')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ee496dd",
   "metadata": {},
   "source": [
    "### Evaluation on Test Set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9353f0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# predict on the test set\n",
    "# predicted classes using default 0.5 threshold\n",
    "y_hat = fit_lr.predict(X_test) \n",
    "\n",
    "#predicted probabilities\n",
    "y_proba = fit_lr.predict_proba(X_test)[:,1] \n",
    "\n",
    "# auc\n",
    "auc=roc_auc_score(y_test, y_proba)\n",
    "\n",
    "# model loss\n",
    "loss= log_loss(y_test, y_hat)\n",
    "\n",
    "print ('the AUC is: {:0.3f}'.format(auc))\n",
    "print ('the logloss is: {:0.3f}'.format(loss))\n",
    "print(\"classification report:\\n \", classification_report(y_test, y_hat, digits=3))\n",
    "print(\"confusion matrix:\\n \")\n",
    "\n",
    "# Create the confusion matrix\n",
    "cm = confusion_matrix(y_test, y_hat)\n",
    "ConfusionMatrixDisplay(confusion_matrix=cm).plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29377237",
   "metadata": {},
   "source": [
    "### Receiver Operating Characteristic Curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c4fee19",
   "metadata": {},
   "outputs": [],
   "source": [
    "fpr, tpr, thresholds = roc_curve(y_test, y_proba, pos_label=1)\n",
    "\n",
    "plt.title('ROC curve')\n",
    "ax1 = plt.plot(fpr, tpr, 'b', label = '%s AUC = %0.3f' % ('', auc), linewidth=2)\n",
    "plt.legend(loc = 'lower right')\n",
    "plt.plot([0, 1], [0, 1],'r--')\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([0, 1])\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.xlabel('False Positive Rate')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a733f903",
   "metadata": {},
   "source": [
    "### Precision-Recall Curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64d4e575",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_proba = fit_lr.predict_proba(X_test)[:,1]\n",
    "\n",
    "precision, recall, thresholds = precision_recall_curve(y_test, y_proba, pos_label=1, sample_weight=None)\n",
    "avg_p = average_precision_score(y_test, y_proba, pos_label=1, sample_weight=None)\n",
    "\n",
    "plt.title('Precision-Recall curve')\n",
    "ax1= plt.plot(precision, recall, 'b', label = '%s AP = %0.3f' % ('', avg_p), linewidth=2)\n",
    "plt.legend(loc = 'lower left')\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([0, 1])\n",
    "plt.ylabel('Precision')\n",
    "plt.xlabel('Recall')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "588e52f7",
   "metadata": {},
   "source": [
    "## Random Forest Algorithm\n",
    "\n",
    "<img src=\"lecture_images/decision_tree.png\" alt=\"title\" width=\"400\" align=\"left\">\n",
    "\n",
    "* Random forest is an esembling algorithm that learns from generating multiple (hundreds-thousands) of decision trees\n",
    "* Randomly selects a subset of features when growing each tree \n",
    "* This approach allows us to average the predictions of the \"crowd\"\n",
    "\n",
    "<b></b>\n",
    "* Many **hyperparameters** that can be tuned during training:\n",
    "<br></br>\n",
    "* `n_estimators`: Number of trees\n",
    "* `max_features`: Number of featues to be considered at each split \n",
    "* `max_depth`: Maximum number of levels in tree (can help regularize the tree to prevent overfitting)\n",
    "* `min_samples_split`: Minimum samples to split a node \n",
    "* `min_samples_leaf`: Minimum number of samples required at each leaf node "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "802d4a35",
   "metadata": {},
   "outputs": [],
   "source": [
    "### tuning RF hyperparameters\n",
    "# Number of trees in random forest\n",
    "n_estimators = [100, 300]\n",
    "# Number of features to consider at every split\n",
    "max_features = [2, 5]\n",
    "# Maximum number of levels in tree\n",
    "max_depth = [5,10]\n",
    "# Minimum number of samples required to split a node\n",
    "min_samples_split = [5, 10]\n",
    "# Minimum number of samples required at each leaf node\n",
    "min_samples_leaf = [10, 15]\n",
    "\n",
    "# create grid of hypterparameter settings\n",
    "# we will permute through this grid\n",
    "param_grid = {'n_estimators': n_estimators,\n",
    "               'max_features': max_features,\n",
    "               'max_depth': max_depth,\n",
    "               'min_samples_split': min_samples_split,\n",
    "               'min_samples_leaf': min_samples_leaf}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1f6a97a",
   "metadata": {},
   "source": [
    "### Hyperparameter Tuning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2e006f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a random forest classifier\n",
    "rf = RandomForestClassifier(criterion='entropy', random_state=12345)\n",
    "\n",
    "# Use random search to find the best hyperparameters\n",
    "grid_search = GridSearchCV(estimator = rf,\n",
    "                           param_grid = param_grid,\n",
    "                           cv = 5,\n",
    "                           scoring = 'roc_auc',\n",
    "                           return_train_score = False,\n",
    "                           n_jobs = -1)\n",
    "\n",
    "# Fit the random search object to the data\n",
    "grid_search.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a885232b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a variable for the best model\n",
    "best_rf = grid_search.best_estimator_\n",
    "\n",
    "# Print the best hyperparameters\n",
    "print('Best hyperparameters:',  grid_search.best_params_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b1e063c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#grid_search.cv_results_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d96bcfe",
   "metadata": {},
   "source": [
    "### Training Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8518e972",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = best_rf.predict(X_train)\n",
    "y_prob = best_rf.predict_proba(X_train)[:,1]\n",
    "\n",
    "accuracy = accuracy_score(y_train, y_pred)\n",
    "precision = precision_score(y_train, y_pred)\n",
    "recall = recall_score(y_train, y_pred)\n",
    "auroc = roc_auc_score(y_train, y_prob)\n",
    "\n",
    "print(\"Accuracy:\", accuracy)\n",
    "print(\"AUC:\", auroc)\n",
    "print(\"Precision:\", precision)\n",
    "print(\"Recall:\", recall)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9befb41a",
   "metadata": {},
   "source": [
    "### Feature Importance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e77a553",
   "metadata": {},
   "outputs": [],
   "source": [
    "# review feature importances on the training set\n",
    "feature_importances = pd.Series(best_rf.feature_importances_, index=X_train.columns).sort_values(ascending=False)\n",
    "\n",
    "# Plot a simple bar chart\n",
    "feature_importances.plot.bar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60105593",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_hat = best_rf.predict(X_train) # predicted classes using default 0.5 threshold\n",
    "y_proba = best_rf.predict_proba(X_train)[:,1] #predicted probabilities\n",
    "auc=roc_auc_score(y_train, y_proba)\n",
    "loss= log_loss(y_train, y_hat)\n",
    "\n",
    "print ('the AUC is: {:0.3f}'.format(auc))\n",
    "print ('the logloss is: {:0.3f}'.format(loss))\n",
    "print(\"classification report:\\n \", classification_report(y_train,y_hat, digits=3))\n",
    "print(\"confusion matrix:\\n \")\n",
    "# Create the confusion matrix\n",
    "cm = confusion_matrix(y_train, y_hat)\n",
    "ConfusionMatrixDisplay(confusion_matrix=cm).plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b53b08d",
   "metadata": {},
   "source": [
    "### Receiver Operating Characteristic Curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "805d2b0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fpr, tpr, thresholds = roc_curve(y_train, y_proba, pos_label=1)\n",
    "\n",
    "plt.title('ROC curve')\n",
    "ax1 = plt.plot(fpr, tpr, 'b', label = '%s AUC = %0.3f' % ('', auc), linewidth=2)\n",
    "plt.legend(loc = 'lower right')\n",
    "plt.plot([0, 1], [0, 1],'r--')\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([0, 1])\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.xlabel('False Positive Rate')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88ed8623",
   "metadata": {},
   "source": [
    "### Precision-Recall Curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90eaecf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_proba = best_rf.predict_proba(X_train)[:,1]\n",
    "\n",
    "precision, recall, thresholds = precision_recall_curve(y_train, y_proba, pos_label=1, sample_weight=None)\n",
    "avg_p = average_precision_score(y_train, y_proba, pos_label=1, sample_weight=None)\n",
    "\n",
    "plt.title('Precision-Recall curve')\n",
    "ax1= plt.plot(precision, recall, 'b', label = '%s AP = %0.3f' % ('', avg_p), linewidth=2)\n",
    "plt.legend(loc = 'lower left')\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([0, 1])\n",
    "plt.ylabel('Precision')\n",
    "plt.xlabel('Recall')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e8e3089",
   "metadata": {},
   "source": [
    "### Evaluation of Test Set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "676b3e0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_hat = best_rf.predict(X_test) # predicted classes using default 0.5 threshold\n",
    "y_proba = best_rf.predict_proba(X_test)[:,1] #predicted probabilities\n",
    "auc=roc_auc_score(y_test, y_proba)\n",
    "loss= log_loss(y_test, y_hat)\n",
    "\n",
    "print ('the AUC is: {:0.3f}'.format(auc))\n",
    "print ('the logloss is: {:0.3f}'.format(loss))\n",
    "print(\"classification report:\\n \", classification_report(y_test,y_hat, digits=3))\n",
    "print(\"confusion matrix:\\n \")\n",
    "# Create the confusion matrix\n",
    "cm = confusion_matrix(y_test, y_hat)\n",
    "ConfusionMatrixDisplay(confusion_matrix=cm).plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8139f850",
   "metadata": {},
   "source": [
    "### Receiver Operating Characteristic Curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c99806ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "fpr, tpr, thresholds = roc_curve(y_test, y_proba, pos_label=1)\n",
    "\n",
    "plt.title('ROC curve')\n",
    "ax1 = plt.plot(fpr, tpr, 'b', label = '%s AUC = %0.3f' % ('', auc), linewidth=2)\n",
    "plt.legend(loc = 'lower right')\n",
    "plt.plot([0, 1], [0, 1],'r--')\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([0, 1])\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.xlabel('False Positive Rate')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0a24bad",
   "metadata": {},
   "source": [
    "### Precision-Recall Curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a453e89f",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_proba = best_rf.predict_proba(X_test)[:,1]\n",
    "\n",
    "precision, recall, thresholds = precision_recall_curve(y_test, y_proba, pos_label=1, sample_weight=None)\n",
    "avg_p = average_precision_score(y_test, y_proba, pos_label=1, sample_weight=None)\n",
    "\n",
    "plt.title('Precision-Recall curve')\n",
    "ax1= plt.plot(precision, recall, 'b', label = '%s AP = %0.3f' % ('', avg_p), linewidth=2)\n",
    "plt.legend(loc = 'lower left')\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([0, 1])\n",
    "plt.ylabel('Precision')\n",
    "plt.xlabel('Recall')"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "kernelspec": {
   "display_name": "mimic4-env",
   "language": "python",
   "name": "mimic4-env"
  },
  "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.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}