pa_pynb_sess6_7_IntroCART.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "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.8.3"
    },
    "colab": {
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/danielbauer1979/CAS_PredMod/blob/main/pa_pynb_sess6_7_IntroCART.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Introduction to Classification and Regression Trees (CART)\n",
        "\n",
        "Dani Bauer, 2022\n",
        "\n",
        "In this tutorial, we will introduce classification and regression trees.  We will start with a basic simulated example with a single predictor that compares tress to (linear) regression.  In particular, we will demonstrate when trees do well and when the don't.  Subsequently, we will revisit our Caravan case study.\n",
        "\n",
        "As usually, let's start with loading the relevant libaries."
      ],
      "metadata": {
        "id": "_Dgf3D6YiA0C"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7vCC-SQb7VRR"
      },
      "source": [
        "import numpy as np \n",
        "import matplotlib.pyplot as plt  \n",
        "import pandas as pd \n",
        "import seaborn as sns\n",
        "import graphviz\n",
        "import pydot\n",
        "from io import StringIO  \n",
        "\n",
        "from sklearn.model_selection import train_test_split, cross_val_score\n",
        "from sklearn.linear_model import LinearRegression\n",
        "from sklearn.tree import DecisionTreeRegressor, DecisionTreeClassifier, export_graphviz\n",
        "from sklearn.metrics import mean_squared_error,confusion_matrix, classification_report"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Regression Trees\n",
        "\n",
        "## Review of Concepts and Maths\n",
        "\n",
        "### Trees\n",
        "\n",
        "As we discuss in the lecture, trees successively split the domains of predictors into two areas (greater and smaller), so that the final tree can be thought of as a model with constant predictions in *square partitions* of the feature domain.  Hence, in building a tree, we need to address two questions:\n",
        "\n",
        "1. <span style=\"color:blue\"> How to come up with the splits that determine the square partitions? </span>\n",
        "\n",
        "2. <span style=\"color:blue\"> How to come up with the predictions for each of the square partitions? </span>\n",
        "\n",
        "The second task is straightforward:  We simply take the average of all predictions in the current partition.  \n",
        "\n",
        "For 2., one relies on a *greedy algorithm*: For each node, we go over all possible variables and all possible split levels, and then choose the combination of variable and split level that mimizes the prediction error.  However, note that the tree thus will generally not minimize the squared prediction error -- finding such a partition is computationally/algorithmically infeasible...\n",
        "\n",
        "### Excursion: The Sigmoid Function\n",
        "\n",
        "In generating our dataset we will use the so-called *sigmoid function*, which will be a key ingredient when we discuss *neural networks*.  It has the advantage that depending on the choice of a parameter, it can mimic a highly linear or a highly non-linear relationship:\n",
        "\n",
        "$$\n",
        "\\sigma_{\\alpha}(x) = \\frac{1}{1+e^{-\\alpha\\,x}}\n",
        "$$\n"
      ],
      "metadata": {
        "id": "oFco9P3jzMD4"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sBGwtLkV7VRS"
      },
      "source": [
        "def sigmoid(x):\n",
        "    return(1 / (1 + np.exp(-x)))\n",
        "x = np.linspace(-5, 5, 1000)\n",
        "plt.figure(figsize=(15,6))\n",
        "plt.plot(x, sigmoid(x), color='b')\n",
        "plt.plot(x, sigmoid(0.5*x), color='r')\n",
        "plt.plot(x, sigmoid(5*x), color='g')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "So we notice that for a small choice of $\\alpha$ (red curve), we obtain a relationship that is almost linear (and, indeed, as $\\alpha$ becomes smaller and smaller, $\\sigma_{\\alpha}$ approaches a perfectly linear function).  Whereas for a larger choice of $\\alpha$ (green curve), we almost obtain a step function (and, indeed, as $\\alpha$ becomes larger and larger, $\\sigma_{\\alpha}$ approaches a step function).\n",
        "\n",
        "## A Simulated Example with One Predictor\n",
        "\n",
        "1. Let's start by generating two datasets based on the sigmoid function with different parameters taken from above. "
      ],
      "metadata": {
        "id": "wyJW2guIzZ1U"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "00OVN3o97VRT"
      },
      "source": [
        "np.random.seed(1)\n",
        "x = 3 * np.random.normal(0, 1, 150)\n",
        "eps = 0.25 * np.random.normal(0, 1, 150)\n",
        "\n",
        "y_1 = sigmoid(0.5 * x) + eps\n",
        "y_2 = sigmoid(5 * x) + eps\n",
        "\n",
        "mydata1 = pd.DataFrame({'y_1':y_1,'x':x})\n",
        "mytraindata1 = mydata1[1:100]\n",
        "mytestdata1 = mydata1[101:150]\n",
        "\n",
        "plt.scatter(mydata1.x,mydata1.y_1)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DDvNPVwD7VRT"
      },
      "source": [
        "mydata2 = pd.DataFrame({'y_2':y_2,'x':x})\n",
        "mytraindata2 = mydata2[1:100]\n",
        "mytestdata2 = mydata2[101:150]\n",
        "plt.scatter(mydata2.x,mydata2.y_2)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "2. And, for comparison, let's fit a linear regression model to both datasets and let's evaluate the out-of-sample error. "
      ],
      "metadata": {
        "id": "umVdXq4ozu7E"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5KPRQ8-T7VRT"
      },
      "source": [
        "X_1train = mytraindata1['x'].values.reshape(-1,1)\n",
        "X_1test = mytestdata1['x'].values.reshape(-1,1)\n",
        "y_1train = mytraindata1['y_1'].values\n",
        "lmfit1 = LinearRegression(fit_intercept=True)\n",
        "lmfit1.fit(X_1train,y_1train)\n",
        "yhat_OOS1 = lmfit1.predict(X_1test)\n",
        "OLS_OOS_MSE1 = sum((mytestdata1['y_1'] - yhat_OOS1)**2)/50\n",
        "OLS_OOS_MSE1"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0SEVVA_07VRU"
      },
      "source": [
        "X_2train = mytraindata2['x'].values.reshape(-1,1)\n",
        "X_2test = mytestdata2['x'].values.reshape(-1,1)\n",
        "y_2train = mytraindata2['y_2'].values\n",
        "lmfit2 = LinearRegression(fit_intercept=True)\n",
        "lmfit2.fit(X_2train,y_2train)\n",
        "yhat_OOS2 = lmfit2.predict(X_2test)\n",
        "OLS_OOS_MSE2 = sum((mytestdata2['y_2'] - yhat_OOS2)**2)/50\n",
        "OLS_OOS_MSE2"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "3. Let's fit a tree to the second dataset. "
      ],
      "metadata": {
        "id": "09U5VR1Fz0ji"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6nffD9ev7VRU"
      },
      "source": [
        "tree2 = DecisionTreeRegressor(max_leaf_nodes=3)\n",
        "tree2.fit(X_2train, y_2train)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "The following function creates images of tree models using pydot, as the package sklearn doesn't offer graphs of the trees\n"
      ],
      "metadata": {
        "id": "KpELVK_oz7DY"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9VpYt9zW7VRU"
      },
      "source": [
        "import pydot\n",
        "from IPython.display import Image\n",
        "def print_tree(estimator, features, class_names=None, filled=True):\n",
        "    tree = estimator\n",
        "    names = features\n",
        "    color = filled\n",
        "    classn = class_names\n",
        "    \n",
        "    dot_data = StringIO()\n",
        "    export_graphviz(estimator, out_file=dot_data, feature_names=features, class_names=classn, filled=filled)\n",
        "    graph = pydot.graph_from_dot_data(dot_data.getvalue())\n",
        "    return(graph)"
      ],
      "execution_count": 8,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gF0SN58g7VRU"
      },
      "source": [
        "graph, = print_tree(tree2, features=['x'])\n",
        "Image(graph.create_png())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Note: Tree algorithms are implemented differently in R and Python, so generally we don't expect exactly the same results from R and Python."
      ],
      "metadata": {
        "id": "yFWlIPQw0DIt"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "4. Let's use cross validation to see what size is optimal, and let's prune accordingly."
      ],
      "metadata": {
        "id": "NdRvpUV40KJ5"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bDmCteFY7VRV"
      },
      "source": [
        "from sklearn.model_selection import cross_val_score\n",
        "OLS_OOS_MSE2 = []\n",
        "cv_score = []\n",
        "for i in range(2,7):\n",
        "    tree2 = DecisionTreeRegressor(max_leaf_nodes=i)\n",
        "    tree2.fit(X_2train, y_2train)\n",
        "    ytreehat2 = tree2.predict(X_2test)\n",
        "    OLS_OOS_MSE2.append(sum((mytestdata2['y_2'] - ytreehat2)**2)/50)\n",
        "    cv_score.append(sum(cross_val_score(tree2, X_2train, y_2train, cv=5)))\n",
        "plt.plot(range(2,7),OLS_OOS_MSE2)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "I8yvrf5NLA03"
      },
      "source": [
        "plt.plot(range(2,7),cv_score)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Python doesn't have automatic pruning algorithm. We can prune ourselves by cross validation with hyper-parameters in the DecisionTreeRegressor function. \n",
        "More info in the  [scikit-docu](https://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeRegressor.html)"
      ],
      "metadata": {
        "id": "cPpNhInG0aAF"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "5. Let's redo for the first dataset."
      ],
      "metadata": {
        "id": "Le9_ijWW0iMk"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "i7hunbTP7VRW"
      },
      "source": [
        "tree1 = DecisionTreeRegressor()\n",
        "tree1.fit(X_1train, y_1train)\n",
        "graph, = print_tree(tree1, features=['x'])\n",
        "Image(graph.create_png())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LK9couni7VRW"
      },
      "source": [
        "OLS_OOS_MSE1 = []\n",
        "cv_score = []\n",
        "for i in range(2,7):\n",
        "    tree1 = DecisionTreeRegressor(max_leaf_nodes=i)\n",
        "    tree1.fit(X_1train, y_1train)\n",
        "    ytreehat1 = tree1.predict(X_1test)\n",
        "    OLS_OOS_MSE1.append(sum((mytestdata1['y_1'] - ytreehat1)**2)/50)\n",
        "    cv_score.append(sum(cross_val_score(tree1, X_1train, y_1train, cv=5)))\n",
        "plt.plot(range(2,7),OLS_OOS_MSE1)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vDGMecqORdJk"
      },
      "source": [
        "plt.plot(range(2,7),cv_score)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Hence, the tree performs very differently between the two datasets:  For the second dataset, the tree provided a very simple (and understandable!) model of the data, and despite the simplicity the prediction accuracy improved relative to the linear regressions.  For the first dataset, on the other hand, the tree is far more complex than the regression model but still provides worse predictions.\n",
        "\n",
        "Hence, overall, trees are well-suited for modeling non-linear relationships -- which explains their relevance as the foundational structural element for more advanced learners."
      ],
      "metadata": {
        "id": "G6X4SK1D0u-p"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Case Study I: Caravan Insurance Purchases\n",
        "\n",
        "Let's go back to the `Caravan` insurance data:"
      ],
      "metadata": {
        "id": "RnUCh9D_06IU"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "joqrQm2gSDdZ"
      },
      "source": [
        "!git clone https://github.com/danielbauer1979/CAS_PredMod.git"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jILOaSki7VRX"
      },
      "source": [
        "Caravan = pd.read_csv('CAS_PredMod/pa_data_caravan.csv', index_col=0)"
      ],
      "execution_count": 16,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "As a reminder, our GLM and GAM models did not get a single true-positive with a cutoff of 50% -- essentially the predictions were all zero.  We obtained an AUC of 0.7456 for the GAM.  The knn classifier did better.  While we cannot obtain an AUC, particularlty for $k=1$ we obtained 45 true classifications, and 10 were correct.  So that makes for a true-postitive rate of roughly 18%.\n",
        "\n",
        "### Classification Tree\n",
        "\n",
        "Let's try a classification tree, with the caveat that we have to change the default parameters. The standard value of the so-called complexity parameter `cp` is insufficient to generate sufficient splits, because a split only happens if there is sufficient heterogeneity in the nodes.  We set it to 0.001 but we can generate an even larger tree by a lower choice."
      ],
      "metadata": {
        "id": "E6c0_HmO1MjV"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IsBCRB477VRX"
      },
      "source": [
        "Caravan.Purchase = Caravan.Purchase=='Yes'\n",
        "test = Caravan.iloc[0:1000]\n",
        "train = Caravan.iloc[1000:len(Caravan)]\n",
        "X = train.drop(columns = ['Purchase'])\n",
        "y = train.Purchase\n",
        "Car_tree_first = DecisionTreeRegressor(max_leaf_nodes=4)\n",
        "Car_tree_first.fit(X, y)\n",
        "graph, = print_tree(Car_tree_first, features=X.columns)\n",
        "Image(graph.create_png())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "The issue with growing the tree is that there are few positives, leading to substantial \"note purity\" even after a few modeling steps. We have to adjust the parameters to build a larger tree:"
      ],
      "metadata": {
        "id": "OFO0ZdDw1Q3F"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "E7P66TyA7VRX"
      },
      "source": [
        "Car_tree = DecisionTreeRegressor(min_samples_split=5,min_impurity_decrease=0.0001)\n",
        "Car_tree.fit(X, y)\n",
        "graph, = print_tree(Car_tree, features=X.columns)\n",
        "Image(graph.create_png())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let's look at the top features by importance scores:"
      ],
      "metadata": {
        "id": "XRuWjoD51rAo"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zHT8Fw9B7VRX"
      },
      "source": [
        "summary_tree = pd.DataFrame({'Features':X.columns,'Importance':Car_tree.feature_importances_}) \n",
        "summary_tree.sort_values(by=['Importance'], ascending=False)[0:10]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "The tree has the following final nodes:"
      ],
      "metadata": {
        "id": "4jL3bgiL1zKp"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EOhl2SU77VRX"
      },
      "source": [
        "Car_tree.tree_.node_count"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let's look at the number of \"purchases\":"
      ],
      "metadata": {
        "id": "7D0OFUhm199v"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UIi7W3sU7VRY"
      },
      "source": [
        "yhat = Car_tree.predict(test.drop(columns = ['Purchase']))\n",
        "pred = yhat > 0.5\n",
        "np.sum(pred)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "And the confusion tavke is:"
      ],
      "metadata": {
        "id": "FkLSnvOS2KD8"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "VZ9J7aNq7VRY"
      },
      "source": [
        "table = pd.DataFrame({'Purchase':test.Purchase,'pred':pred})\n",
        "table.groupby(['Purchase','pred']).size().unstack('Purchase')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Much better!"
      ],
      "metadata": {
        "id": "FWJTIsQF2OYQ"
      }
    }
  ]
}