eicu_decision_trees.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "05-prediction.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "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.3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/MIT-LCP/bidmc-datathon/blob/master/eicu_decision_trees.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "T3wdKZCPklNq"
      },
      "source": [
        "# eICU Collaborative Research Database\n",
        "\n",
        "# Notebook 5: Prediction with decision trees\n",
        "\n",
        "This notebook explores how a decision trees can be trained to predict in-hospital mortality of patients.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "rG3HrM7GkwCH"
      },
      "source": [
        "## Load libraries and connect to the database"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "s-MoFA6NkkbZ",
        "colab": {}
      },
      "source": [
        "# Import libraries\n",
        "import os\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# model building\n",
        "from sklearn import ensemble, impute, metrics, preprocessing, tree\n",
        "from sklearn.model_selection import cross_val_score, train_test_split\n",
        "from sklearn.pipeline import Pipeline\n",
        "\n",
        "# Make pandas dataframes prettier\n",
        "from IPython.display import display, HTML, Image\n",
        "plt.rcParams.update({'font.size': 20})\n",
        "%matplotlib inline\n",
        "plt.style.use('ggplot')\n",
        "\n",
        "# Access data using Google BigQuery.\n",
        "from google.colab import auth\n",
        "from google.cloud import bigquery"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "jyBV_Q9DkyD3",
        "colab": {}
      },
      "source": [
        "# authenticate\n",
        "auth.authenticate_user()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "cF1udJKhkzYq",
        "colab": {}
      },
      "source": [
        "# Set up environment variables\n",
        "project_id='bidmc-datathon'\n",
        "os.environ[\"GOOGLE_CLOUD_PROJECT\"]=project_id"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "xGurBAQIUDTt"
      },
      "source": [
        "To make our lives easier, we'll also install and import a set of helper functions from the `datathon2` package. We will be using the following functions from the package:\n",
        "- `plot_model_pred_2d`: to visualize our data, helping to display a class split assigned by a tree vs the true class.\n",
        "- `run_query()`: to run an SQL query against our BigQuery database and assign the results to a dataframe. \n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "GDEewAlvk0oT",
        "colab": {}
      },
      "source": [
        "!pip install pip --upgrade\n",
        "!pip install glowyr\n",
        "!pip install pandas_profiling --upgrade"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "JM6O5GPAUI89",
        "colab": {}
      },
      "source": [
        "import glowyr as dtn\n",
        "import pydotplus\n",
        "from tableone import TableOne\n",
        "from pandas_profiling import ProfileReport"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "hq_09Hh-y17k"
      },
      "source": [
        "In this notebook we'll be looking at tree models, so we'll now install a package for visualizing these models."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "jBMOwgwszGOw",
        "colab": {}
      },
      "source": [
        "!apt-get install graphviz -y"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "LgcRCqxCk3HC"
      },
      "source": [
        "## Load the patient cohort\n",
        "\n",
        "Let's extract a cohort of patients admitted to the ICU from the emergency department.  We link demographics data from the `patient` table to severity of illness score data in the `apachepatientresult` table. We exclude readmissions and neurological patients to help create a population suitable for our demonstration."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "ReCl7-aek1-k",
        "colab": {}
      },
      "source": [
        "# Link the patient, apachepatientresult, and apacheapsvar tables on patientunitstayid\n",
        "# using an inner join.\n",
        "\n",
        "%%bigquery cohort\n",
        "\n",
        "SELECT p.unitadmitsource, p.gender, p.age, p.unittype, p.unitstaytype, \n",
        "    a.actualhospitalmortality, a.acutePhysiologyScore, a.apacheScore\n",
        "FROM `physionet-data.eicu_crd_demo.patient` p\n",
        "INNER JOIN `physionet-data.eicu_crd_demo.apachepatientresult` a\n",
        "ON p.patientunitstayid = a.patientunitstayid\n",
        "WHERE a.apacheversion LIKE 'IVa'\n",
        "AND LOWER(p.unitadmitsource) LIKE \"%emergency%\"\n",
        "AND LOWER(p.unitstaytype) LIKE \"admit%\"\n",
        "AND LOWER(p.unittype) NOT LIKE \"%neuro%\";"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "yxLctVBpk9sO",
        "colab": {}
      },
      "source": [
        "cohort.head()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "NPlwRV2buYb1"
      },
      "source": [
        "## Preparing the data for analysis\n",
        "\n",
        "Before continuing, we want to review our data, paying attention to factors such as:\n",
        "- data types (for example, are values recorded as characters or numerical values?) \n",
        "- missing data\n",
        "- distribution of values"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "v3OJ4LDvueKu",
        "colab": {}
      },
      "source": [
        "# dataset info\n",
        "print(cohort.info())"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uZJXNgUjTR2h",
        "colab_type": "text"
      },
      "source": [
        "We can use a tool from `pandas_profiling` called \"ProfileReport\" to get an idea of the quality of the data."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TQNvhIWsTRVt",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "profile = ProfileReport(cohort, title='Pandas Profiling Report', html={'style':{'full_width':True}})"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "47R7EPMMTX8X",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# create an interactive widget which allows exploration of the data\n",
        "profile.to_widgets()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "s4wQ6o_RvLph",
        "colab": {}
      },
      "source": [
        "# Encode the categorical data\n",
        "encoder = preprocessing.LabelEncoder()\n",
        "cohort['gender_code'] = encoder.fit_transform(cohort['gender'])\n",
        "cohort['actualhospitalmortality_code'] = encoder.fit_transform(cohort['actualhospitalmortality'])"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "_1LYcNUdjQA5"
      },
      "source": [
        "In the eICU Collaborative Research Database, ages >89 years have been removed to comply with data sharing regulations. We will need to decide how to handle these ages. For simplicity, we will assign an age of 91.5 years to these patients."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "4ogi_ns-ylnP",
        "colab": {}
      },
      "source": [
        "# Handle the deidentified ages\n",
        "cohort['age'] = pd.to_numeric(cohort['age'], downcast='integer', errors='coerce')\n",
        "cohort['age'] = cohort['age'].fillna(value=91.5)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "77M0QJQ5wcPQ",
        "colab": {}
      },
      "source": [
        "# Preview the encoded data\n",
        "cohort[['gender','gender_code']].head()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "GqvwTNPN3KZz",
        "colab": {}
      },
      "source": [
        "# Check the outcome variable\n",
        "cohort['actualhospitalmortality_code'].unique()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "OdGX1qWdkTgY"
      },
      "source": [
        "Now let's use the [tableone package](https://doi.org/10.1093/jamiaopen/ooy012\n",
        ") to review our dataset."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "gIIsthy1WK3i",
        "colab": {}
      },
      "source": [
        "# View summary statistics\n",
        "pd.set_option('display.max_rows', 500)\n",
        "TableOne(cohort,groupby='actualhospitalmortality')"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "IGtKlTG1gvRf"
      },
      "source": [
        "From these summary statistics, we can see that the average age is higher in the group of patients who do not survive. What other differences do you see?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "ze7y5J4Ioz8u"
      },
      "source": [
        "## Creating our train and test sets\n",
        "\n",
        "We only focus on two variables for our analysis, age and acute physiology score. Limiting ourselves to two variables will make it easier to visualize our models."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "i5zXkn_AlDJW",
        "colab": {}
      },
      "source": [
        "features = ['age','acutePhysiologyScore']\n",
        "outcome = 'actualhospitalmortality_code'\n",
        "\n",
        "X = cohort[features]\n",
        "y = cohort[outcome]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "IHhIgDUwocmA",
        "colab": {}
      },
      "source": [
        "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=10)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "NvQWkuY6nkZ8",
        "colab": {}
      },
      "source": [
        "# Review the number of cases in each set\n",
        "print(\"Train data: {}\".format(len(X_train)))\n",
        "print(\"Test data: {}\".format(len(X_test)))"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "b2waK5qBqanC"
      },
      "source": [
        "## Decision trees\n",
        "\n",
        "Let's build the simplest tree model we can think of: a classification tree with only one split. Decision trees of this form are commonly referred to under the umbrella term Classification and Regression Trees (CART) [1]. \n",
        "\n",
        "While we will only be looking at classification here, regression isn't too different. After grouping the data (which is essentially what a decision tree does), classification involves assigning all members of the group to the majority class of that group during training. Regression is the same, except you would assign the average value, not the majority. \n",
        "\n",
        "In the case of a decision tree with one split, often called a \"stump\", the model will partition the data into two groups, and assign classes for those two groups based on majority vote. There are many parameters available for the DecisionTreeClassifier class; by specifying max_depth=1 we will build a decision tree with only one split - i.e. of depth 1.\n",
        "\n",
        "[1] L. Breiman, J. Friedman, R. Olshen, and C. Stone. Classification and Regression Trees. Wadsworth, Belmont, CA, 1984."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "RlG3N3OYBqAm",
        "colab": {}
      },
      "source": [
        "# specify max_depth=1 so we train a stump, i.e. a tree with only 1 split\n",
        "mdl = tree.DecisionTreeClassifier(max_depth=1)\n",
        "\n",
        "# fit the model to the data - trying to predict y from X\n",
        "mdl = mdl.fit(X_train,y_train)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "8RlioUw8B_0O"
      },
      "source": [
        "Our model is so simple that we can look at the full decision tree."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "G2t9Nz8pBqEb",
        "colab": {}
      },
      "source": [
        "graph = dtn.create_graph(mdl,feature_names=features)\n",
        "Image(graph.create_png())"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "E-iPwWWKCGY9"
      },
      "source": [
        "Here we see three nodes: a node at the top, a node in the lower left, and a node in the lower right.\n",
        "\n",
        "The top node is the root of the tree: it contains all the data. Let's read this node bottom to top:\n",
        "- `value = [384, 44]`:  Current class balance. There are 384 observations of class 0 and 44 observations of class 1.\n",
        "- `samples = 428`:  Number of samples assessed at this node.\n",
        "- `gini = 0.184`: Gini impurity, a measure of \"impurity\". The higher the value, the bigger the mix of classes. A 50/50 split of two classes would result in an index of 0.5.\n",
        "- `acutePhysiologyScore <=78.5`: Decision rule learned by the node. In this case, patients with a score of <= 78.5 are moved into the left node and >78.5 to the right. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "KS0UcZqUeJKz"
      },
      "source": [
        "The gini impurity is actually used by the algorithm to determine a split. The model evaluates every feature (in our case, age and score) at every possible split (46, 47, 48..) to find the point with the lowest gini impurity in two resulting nodes. \n",
        "\n",
        "The approach is referred to as \"greedy\" because we are choosing the optimal split given our current state. Let's take a closer look at our decision boundary."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "uXl22sNTtpHa",
        "colab": {}
      },
      "source": [
        "# look at the regions in a 2d plot\n",
        "# based on scikit-learn tutorial plot_iris.html\n",
        "plt.figure(figsize=[10,8])\n",
        "dtn.plot_model_pred_2d(mdl, X_train, y_train, \n",
        "                       title=\"Decision tree (depth 1)\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "25zSX-inCNOJ"
      },
      "source": [
        "In this plot we can see the decision boundary on the y-axis, separating the predicted classes. The true classes are indicated at each point. Where the background and point colours are mismatched, there has been misclassification. Of course we are using a very simple model. Let's see what happens when we increase the depth."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "ZuO62CL3CSGm",
        "colab": {}
      },
      "source": [
        "mdl = tree.DecisionTreeClassifier(max_depth=5)\n",
        "mdl = mdl.fit(X_train,y_train)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "A88Vi83LCSJ6",
        "colab": {}
      },
      "source": [
        "plt.figure(figsize=[10,8])\n",
        "dtn.plot_model_pred_2d(mdl, X_train, y_train,\n",
        "                      title=\"Decision tree (depth 5)\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "B88XlKDtCYmn"
      },
      "source": [
        "Now our tree is more complicated! We can see a few vertical boundaries as well as the horizontal one from before. Some of these we may like, but some appear unnatural. Let's look at the tree itself."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "V1VLrOJJCcWo",
        "colab": {}
      },
      "source": [
        "graph = dtn.create_graph(mdl,feature_names=features)\n",
        "Image(graph.create_png())"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "Ton_EnvFqHIO"
      },
      "source": [
        "Looking at the tree, we can see that there are some very specific rules. Consider our patient aged 65 years with an acute physiology score of 87. From the top of the tree, we would work our way down:\n",
        "\n",
        "- acutePhysiologyScore <= 78.5? No.\n",
        "- acutePhysiologyScore <= 106.5? Yes.\n",
        "- age <= 75.5? Yes\n",
        "- age <= 66. Yes.\n",
        "- age <= 62.5? No. \n",
        "\n",
        "This leads us to our single node with a gini impurity of 0. Having an entire rule based upon this one observation seems silly, but it is perfectly logical as at the moment. The only objective the algorithm cares about is minimizing the gini impurity. \n",
        "\n",
        "We are at risk of overfitting our data! This is where \"pruning\" comes in."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "VvsNIjCDDIo_",
        "colab": {}
      },
      "source": [
        "# let's prune the model and look again\n",
        "mdl = dtn.prune(mdl, min_samples_leaf = 10)\n",
        "graph = dtn.create_graph(mdl,feature_names=features)\n",
        "Image(graph.create_png())  "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "8pRzzV2VvdxP"
      },
      "source": [
        "Above, we can see that our second tree is (1) smaller in depth, and (2) never splits a node with <= 10 samples. We can look at the decision surface for this tree:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "5LyGDz-Cr-mU",
        "colab": {}
      },
      "source": [
        "plt.figure(figsize=[10,8])\n",
        "dtn.plot_model_pred_2d(mdl, X_train, y_train, title=\"Pruned decision tree\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "xAnqmD_Dv_dh"
      },
      "source": [
        "Our pruned decision tree has a much more intuitive boundary, but does make some errors. We have reduced our performance in an effort to simplify the tree. This is the classic machine learning problem of trading off complexity with error.\n",
        "\n",
        "Note that, in order to do this, we \"invented\" the minimum samples per leaf node of 10. Why 10? Why not 5? Why not 20? The answer is: it depends on the dataset. Heuristically choosing these parameters can be time consuming, and we will see later on how gradient boosting elegantly handles this task."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "2EFINpj-wD7H"
      },
      "source": [
        "## Decision trees have high \"variance\"\n",
        "\n",
        "Before we move on to boosting, it will be useful to demonstrate how decision trees have high \"variance\". In this context, variance refers to a property of some models to have a wide range of performance given random samples of data. Let's take a look at randomly slicing the data we have too see what that means."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "JT7fuuj6vjKB",
        "colab": {}
      },
      "source": [
        "np.random.seed(123)\n",
        "\n",
        "fig = plt.figure(figsize=[12,3])\n",
        "\n",
        "for i in range(3):\n",
        "    ax = fig.add_subplot(1,3,i+1)\n",
        "\n",
        "    # generate indices in a random order\n",
        "    idx = np.random.permutation(X_train.shape[0])\n",
        "    \n",
        "    # only use the first 50\n",
        "    idx = idx[:50]\n",
        "    X_temp = X_train.iloc[idx]\n",
        "    y_temp = y_train.values[idx]\n",
        "    \n",
        "    # initialize the model\n",
        "    mdl = tree.DecisionTreeClassifier(max_depth=5)\n",
        "    \n",
        "    # train the model using the dataset\n",
        "    mdl = mdl.fit(X_temp, y_temp)\n",
        "    txt = 'Random sample {}'.format(i)\n",
        "    dtn.plot_model_pred_2d(mdl, X_temp, y_temp, title=txt)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "j6VTIDr-yRRZ"
      },
      "source": [
        "Above we can see that we are using random subsets of data, and as a result, our decision boundary can change quite a bit. As you could guess, we actually don't want a model that randomly works well and randomly works poorly, so you may wonder why this is useful. \n",
        "\n",
        "The trick is that by combining many of instances of \"high variance\" classifiers (decision trees), we can end up with a single classifier with low variance. There is an old joke: two farmers and a statistician go hunting. They see a deer: the first farmer shoots, and misses to the left. The next farmer shoots, and misses to the right. The statistician yells \"We got it!!\".\n",
        "\n",
        "While it doesn't quite hold in real life, it turns out that this principle does hold for decision trees. Combining them in the right way ends up building powerful models."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "iWnKvx6myf9Z"
      },
      "source": [
        "## Boosting\n",
        "\n",
        "The premise of boosting is the combination of many weak learners to form a single \"strong\" learner. In a nutshell, boosting involves building a models iteratively. At each step we focus on the data on which we performed poorly. \n",
        "\n",
        "In our context, we'll use decision trees, so the first step would be to build a tree using the data. Next, we'd look at the data that we misclassified, and re-weight the data so that we really wanted to classify those observations correctly, at a cost of maybe getting some of the other data wrong this time. Let's see how this works in practice."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "YJWxu0bTwRzD",
        "colab": {}
      },
      "source": [
        "# build the model\n",
        "clf = tree.DecisionTreeClassifier(max_depth=1)\n",
        "mdl = ensemble.AdaBoostClassifier(base_estimator=clf,n_estimators=6)\n",
        "mdl = mdl.fit(X_train,y_train)\n",
        "\n",
        "# plot each individual decision tree\n",
        "fig = plt.figure(figsize=[12,6])\n",
        "for i, estimator in enumerate(mdl.estimators_):\n",
        "    ax = fig.add_subplot(2,3,i+1)\n",
        "    txt = 'Tree {}'.format(i+1)\n",
        "    dtn.plot_model_pred_2d(estimator, X_train, y_train, title=txt)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "5zNfvDjTzh2U"
      },
      "source": [
        "Looking at our example above, we can see that the first iteration builds the exact same simple decision tree as we had seen earlier. This makes sense. It is using the entire dataset with no special weighting. \n",
        "\n",
        "In the next iteration we can see the model shift. It misclassified several observations in class 1, and now these are the most important observations. Consequently, it picks the boundary that, while prioritizing correctly classifies these observations, still tries to best classify the rest of the data too. \n",
        "\n",
        "The iteration process continues, until the model is apparently creating boundaries to capture just one or two observations (see, for example, Tree 6 on the bottom right). \n",
        "\n",
        "One important point is that each tree is weighted by its global error. So, for example, Tree 6 would carry less weight in the final model. It is clear that we wouldn't want Tree 6 to carry the same importance as Tree 1, when Tree 1 is doing so much better overall. It turns out that weighting each tree by the inverse of its error is a pretty good way to do this.\n",
        "\n",
        "Let's look at final model's decision surface.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "3pVG5ytfzp_B",
        "colab": {}
      },
      "source": [
        "# plot the final prediction\n",
        "plt.figure(figsize=[9,5])\n",
        "txt = 'Boosted tree (final decision surface)'\n",
        "dtn.plot_model_pred_2d(mdl, X_train, y_train, title=txt)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "YRGRFjRgz26h"
      },
      "source": [
        "And that's AdaBoost! There are a few tricks we have glossed over here, but you understand the general principle. Now we'll move on to a different approach. With boosting, we iteratively changed the dataset to have new trees focus on the \"difficult\" observations. The next approach we discuss is similar as it also involves using changed versions of our dataset to build new trees."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "EFNDNsIpfP7j"
      },
      "source": [
        "## Bagging\n",
        "\n",
        "Bootstrap aggregation, or \"Bagging\", is another form of *ensemble learning* where we aim to build a single good model by combining many models together. With AdaBoost, we modified the data to focus on hard to classify observations. We can imagine this as a form of resampling the data for each new tree. For example, say we have three observations: A, B, and C, `[A, B, C]`. If we correctly classify observations `[A, B]`, but incorrectly classify `C`, then AdaBoost involves building a new tree that focuses on `C`. Equivalently, we could say AdaBoost builds a new tree using the dataset `[A, B, C, C, C]`, where we have *intentionally* repeated observation `C` 3 times so that the algorithm thinks it is 3 times as important as the other observations. Makes sense?\n",
        "\n",
        "Bagging involves the same approach, except we don't selectively choose which observations to focus on, but rather we *randomly select subsets of data each time*. As you can see, while this is a similar process to AdaBoost, the concept is quite different. Whereas before we aimed to iteratively improve our overall model with new trees, we now build trees on what we hope are independent datasets.\n",
        "\n",
        "Let's take a step back, and think about a practical example. Say we wanted a good model of heart disease. If we saw researchers build a model from a dataset of patients from their hospital, we would be happy. If they then acquired a new dataset from new patients, and built a new model, we'd be inclined to feel that the combination of the two models would be better than any one individually. This exact scenario is what bagging aims to replicate, except instead of actually going out and collecting new datasets, we instead use bootstrapping to create new sets of data from our current dataset. If you are unfamiliar with bootstrapping, you can treat it as \"magic\" for now (and if you are familiar with the bootstrap, you already know that it is magic).\n",
        "\n",
        "Let's take a look at a simple bootstrap model."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "JrXAspvrzv8x",
        "colab": {}
      },
      "source": [
        "np.random.seed(321)\n",
        "clf = tree.DecisionTreeClassifier(max_depth=5)\n",
        "mdl = ensemble.BaggingClassifier(base_estimator=clf, n_estimators=6)\n",
        "mdl = mdl.fit(X_train, y_train)\n",
        "\n",
        "fig = plt.figure(figsize=[12,6])\n",
        "for i, estimator in enumerate(mdl.estimators_):    \n",
        "    ax = fig.add_subplot(2,3,i+1)\n",
        "    txt = 'Tree {}'.format(i+1)\n",
        "    dtn.plot_model_pred_2d(estimator, X_train, y_train, \n",
        "                           title=txt)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "s3kKUPORfW9F"
      },
      "source": [
        "We can see that each individual tree is quite variable. This is a result of using a random set of data to train the classifier."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "w_D7_-0HfVMy",
        "colab": {}
      },
      "source": [
        "# plot the final prediction\n",
        "plt.figure(figsize=[8,5])\n",
        "txt = 'Bagged tree (final decision surface)'\n",
        "dtn.plot_model_pred_2d(mdl, X_train, y_train, title=txt)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "AOFnG0r6faLS"
      },
      "source": [
        "Not bad! Of course, since this is a simple dataset, we are not seeing that many dramatic changes between different models. Don't worry, we'll quantitatively evaluate them later. \n",
        "\n",
        "Next up, a minor addition creates one of the most popular models in machine learning."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "aiqrVfYtfcYk"
      },
      "source": [
        "## Random Forest\n",
        "\n",
        "In the previous example, we used bagging to randomly resample our data to generate \"new\" datasets. The Random Forest takes this one step further: instead of just resampling our data, we also select only a fraction of the features to include. \n",
        "\n",
        "It turns out that this subselection tends to improve the performance of our models. The odds of an individual being very good or very bad is higher (i.e. the variance of the trees is increased), and this ends up giving us a final model with better overall performance (lower bias).\n",
        "\n",
        "Let's train the model."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "u27LS36_fglG",
        "colab": {}
      },
      "source": [
        "np.random.seed(321)\n",
        "mdl = ensemble.RandomForestClassifier(max_depth=5, n_estimators=6, max_features=1)\n",
        "mdl = mdl.fit(X_train,y_train)\n",
        "\n",
        "fig = plt.figure(figsize=[12,6])\n",
        "for i, estimator in enumerate(mdl.estimators_):    \n",
        "    ax = fig.add_subplot(2,3,i+1)\n",
        "    txt = 'Tree {}'.format(i+1)\n",
        "    dtn.plot_model_pred_2d(estimator, X_train, y_train, title=txt)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "5aG0PI8lruGN",
        "colab": {}
      },
      "source": [
        "plt.figure(figsize=[9,5])\n",
        "txt = 'Random forest (final decision surface)'\n",
        "dtn.plot_model_pred_2d(mdl, X_train, y_train, title=txt)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "2KmJuztXfjzm"
      },
      "source": [
        "Again, the visualization doesn't *really* show us the power of Random Forests, but we'll quantitatively evaluate them soon enough.\n",
        "\n",
        "Last, and not least, we move on to gradient boosting."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "LTP8zFIofl2v"
      },
      "source": [
        "## Gradient Boosting\n",
        "\n",
        "Gradient boosting, our last topic, elegantly combines concepts from the previous methods. As a \"boosting\" method, gradient boosting involves iteratively building trees, aiming to improve upon misclassifications of the previous tree. Gradient boosting also borrows the concept of sub-sampling the variables (just like Random Forests), which can help to prevent overfitting.\n",
        "\n",
        "While it is hard to express in this non-technical tutorial, the biggest innovation in gradient boosting is that it provides a unifying mathematical framework for boosting models. The approach explicitly casts the problem of building a tree as an optimization problem, defining mathematical functions for how well a tree is performing (which we had before) *and* how complex a tree is. In this light, one can actually treat AdaBoost as a \"special case\" of gradient boosting, where the loss function is chosen to be the exponential loss.\n",
        "\n",
        "Let's build a gradient boosting model."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "L_QVZ9oNfnqk",
        "colab": {}
      },
      "source": [
        "np.random.seed(321)\n",
        "mdl = ensemble.GradientBoostingClassifier(n_estimators=10)\n",
        "mdl = mdl.fit(X_train, y_train)\n",
        "\n",
        "plt.figure(figsize=[9,5])\n",
        "txt = 'Gradient boosted tree (final decision surface)'\n",
        "dtn.plot_model_pred_2d(mdl, X_train, y_train, title=txt)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "tcCzP4gAsd7L"
      },
      "source": [
        "## Comparing model performance\n",
        "\n",
        "We've now learned the basics of the various tree methods and have visualized most of them. Let's finish by comparing the performance of our models on our held-out test data. Our goal, remember, is to predict whether or not a patient will survive their hospital stay using the patient's age and acute physiology score computed on the first day of their ICU stay."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "tQST4TQAtHmU",
        "colab": {}
      },
      "source": [
        "clf = dict()\n",
        "clf['Decision Tree'] = tree.DecisionTreeClassifier(criterion='entropy', splitter='best').fit(X_train,y_train)\n",
        "clf['Gradient Boosting'] = ensemble.GradientBoostingClassifier(n_estimators=10).fit(X_train, y_train)\n",
        "clf['Random Forest'] = ensemble.RandomForestClassifier(n_estimators=10).fit(X_train, y_train)\n",
        "clf['Bagging'] =  ensemble.BaggingClassifier(n_estimators=10).fit(X_train, y_train)\n",
        "clf['AdaBoost'] =  ensemble.AdaBoostClassifier(n_estimators=10).fit(X_train, y_train)\n",
        "\n",
        "fig = plt.figure(figsize=[10,10])\n",
        "\n",
        "print('AUROC\\tModel')\n",
        "for i, curr_mdl in enumerate(clf):    \n",
        "    yhat = clf[curr_mdl].predict_proba(X_test)[:,1]\n",
        "    score = metrics.roc_auc_score(y_test, yhat)\n",
        "    print('{:0.3f}\\t{}'.format(score, curr_mdl))\n",
        "    ax = fig.add_subplot(3,2,i+1)\n",
        "    dtn. plot_model_pred_2d(clf[curr_mdl], X_test, y_test, title=curr_mdl)\n",
        "    "
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "osr6iM6ltLAP"
      },
      "source": [
        "Here we can see that quantitatively, gradient boosting has produced the highest discrimination among all the models (~0.91). You'll see that some of the models appear to have simpler decision surfaces, which tends to result in improved generalization on a held-out test set (though not always!).\n",
        "\n",
        "To make appropriate comparisons, we should calculate 95% confidence intervals on these performance estimates. This can be done a number of ways. A simple but effective approach is to use bootstrapping, a resampling technique. In bootstrapping, we generate multiple datasets from the test set (allowing the same data point to be sampled multiple times). Using these datasets, we can then estimate the confidence intervals."
      ]
    }
  ]
}