lab_boosting.ipynb

{
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3", 
      "name": "python3", 
      "language": "python"
    }, 
    "language_info": {
      "mimetype": "text/x-python", 
      "nbconvert_exporter": "python", 
      "name": "python", 
      "file_extension": ".py", 
      "version": "3.6.8", 
      "pygments_lexer": "ipython3", 
      "codemirror_mode": {
        "version": 3, 
        "name": "ipython"
      }
    }
  }, 
  "nbformat": 4, 
  "nbformat_minor": 2, 
  "cells": [
    {
      "source": [
        "# Boosting"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "source": [
        "We will look here into the practicalities of boosted trees. "
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 1, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline\n", 
        "import numpy as np\n", 
        "import scipy as sp\n", 
        "import matplotlib as mpl\n", 
        "import matplotlib.cm as cm\n", 
        "import matplotlib.pyplot as plt\n", 
        "import pandas as pd\n", 
        "pd.set_option('display.width', 500)\n", 
        "pd.set_option('display.max_columns', 100)\n", 
        "pd.set_option('display.notebook_repr_html', True)\n", 
        "import seaborn.apionly as sns\n", 
        "sns.set_style(\"whitegrid\")\n", 
        "sns.set_context(\"poster\")"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "## Dataset\n", 
        "\n", 
        "First, the data. We will be attempting to predict the presidential election results (at the county level) from 2016, measured as 'votergap' = (trump - clinton) in percentage points, based mostly on demographic features of those counties.  Let's quick take a peak at the data:"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 2, 
      "cell_type": "code", 
      "source": [
        "elect_df = pd.read_csv(\"data/county_level_election.csv\")\n", 
        "elect_df.head()"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 3, 
      "cell_type": "code", 
      "source": [
        "from sklearn.model_selection import train_test_split"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true
      }
    }, 
    {
      "execution_count": 11, 
      "cell_type": "code", 
      "source": [
        "# split 80/20 train-test\n", 
        "X = elect_df[['population','hispanic','minority','female','unemployed','income','nodegree','bachelor','inactivity','obesity','density','cancer']]\n", 
        "response = elect_df['votergap']\n", 
        "Xtraindf, Xtestdf, ytrain, ytest = train_test_split(X,response,test_size=0.2, random_state=1983)\n"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": true
      }
    }, 
    {
      "execution_count": 12, 
      "cell_type": "code", 
      "source": [
        "plt.hist(ytrain)\n", 
        "Xtraindf.hist(column=['minority', 'population','hispanic','female']);"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "How would you describe these variables?"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 13, 
      "cell_type": "code", 
      "source": [
        "x = Xtraindf['minority'].values\n", 
        "o = np.argsort(x)\n", 
        "x = x[o]\n", 
        "y = ytrain.values\n", 
        "y = y[o]\n", 
        "plt.plot(np.log(x),y,'.')\n", 
        "xx = np.log(x).reshape(-1,1)"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "## Boosting Regression Trees\n", 
        "\n", 
        " Gradient Bossting is very state of the art, and has major connections to logistic regression, gradient descent in a functional space, and search in information space. See Schapire and Freund's MIT Press book for details (Google is a wonderful thing).\n", 
        "\n", 
        "But briefly, let us cover the idea here. The idea is that we will use a bunch of weak 'learners' (aka, models) which are fit sequentially. The first one fits the signal, the second one the first model's residual, the third the second model's residual, and so on. At each stage we upweight the places that our previous model did badly on. First let us illustrate."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 17, 
      "cell_type": "code", 
      "source": [
        "from sklearn.tree import DecisionTreeRegressor\n", 
        "from sklearn.ensemble import AdaBoostRegressor\n"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 20, 
      "cell_type": "code", 
      "source": [
        "# code from http://nbviewer.jupyter.org/github/pprett/pydata-gbrt-tutorial/blob/master/gbrt-tutorial.ipynb\n", 
        "import time\n", 
        "from IPython import display\n", 
        "plt.figure(figsize=(12, 8))\n", 
        "plt.plot(xx, y, '.');\n", 
        "counter = 0\n", 
        "estab = AdaBoostRegressor(base_estimator=DecisionTreeRegressor(max_depth=1),\n", 
        "                          n_estimators=500, learning_rate=1)\n", 
        "estab.fit(xx, y)\n", 
        "staged_predict_generator = estab.staged_predict(xx)\n", 
        "for stagepred in staged_predict_generator:\n", 
        "    counter = counter + 1\n", 
        "    if counter in [1, 2, 4, 8, 10, 50, 100, 200]:\n", 
        "        plt.plot(xx, stagepred, alpha=0.7, label=str(counter), lw=4)\n", 
        "        plt.legend();\n", 
        "        display.display(plt.gcf())\n", 
        "        display.clear_output(wait=True)\n", 
        "        time.sleep(2)"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "## Now on the whole data set"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 22, 
      "cell_type": "code", 
      "source": [
        "estab2 = AdaBoostRegressor(base_estimator=DecisionTreeRegressor(max_depth=2),\n", 
        "                          n_estimators=1000, learning_rate=0.3)\n", 
        "estab2.fit(Xtraindf, ytrain)\n", 
        "print(estab2.score(Xtraindf, ytrain))\n", 
        "print(estab2.score(Xtestdf, ytest))\n", 
        "\n", 
        "staged_scores = estab2.staged_score(Xtest, ytest)\n", 
        "i=1\n", 
        "scores=[]\n", 
        "for stagesc in staged_scores:\n", 
        "    scores=np.append(scores,stagesc)\n", 
        "    i = i+1\n", 
        "plt.plot(scores);"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "Ok, so this demonstration helps us understand some things about boosting.\n", 
        "\n", 
        "- `n_estimators` is the number of trees, and thus the stage in the fitting. It also controls the complexity for us. The more trees we have the more we fit to the tiny details.\n", 
        "- `staged_predict` gives us the prediction at each step\n", 
        "- once again `max_depth` from the underlying decision tree tells us the depth of the tree. But here it tells us the amount of features interactions we have, not just the scale of our fit. But clearly it increases the variance again.\n", 
        "\n", 
        "Ideas from decision trees remain. For example, increase `min_samples_leaf` or decrease `max_depth` to reduce variance and increase the bias."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "source": [
        "EXERCISE\n", 
        "\n", 
        "1. What do you expect to happen if you increase `max_depth` to 5?  Edit the code above to explore the result.\n", 
        "2. What do you expect to happen if you put `max_depth` back to 1 and decrease the `learning_rate` to 0.1?  Edit the code above to explore the result.\n", 
        "3. Do a little work to find some sort of 'best' values of `max_depth` and `learning_rate`.  Does this result make sense?"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "source": [
        "**Answers**\n", 
        "1. Increasing `max_depth` for each model does not really improve things.  It will make each model more complex, which may not be a good thing.  We want lots of little models that each are not overfit, so allowing any single model to have a lot of complexity (and thus overfit) could be a problem.\n", 
        "2. The model will progress much more slowly to fitting to the data well.  This may be a good thing: this is essentially a form of 'regularization' in each model and will help prevent overfitting to any one model.  It allows future models to have more 'error' to be able to explain.\n", 
        "3. See code below.\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 25, 
      "cell_type": "code", 
      "source": [
        "from collections import OrderedDict\n", 
        "from sklearn.model_selection import GridSearchCV\n", 
        "\n", 
        "param_dict3 = OrderedDict(\n", 
        "    n_estimators = [200],\n", 
        "    learning_rate = [0.2, 0.4, 0.6, 0.8]\n", 
        ")\n", 
        "\n", 
        "# your code here\n"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 164, 
      "cell_type": "code", 
      "source": [
        "best_scores"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "### Deviance"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "source": [
        "A deviance plot can be used to compare train and test errors against the number of iterations.\n", 
        "\n", 
        "Training error (deviance, related to the KL-divergence) is stored in est.train_score_\n", 
        "Test error is computed using est.staged_predict (this uses est.loss_)"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": 33, 
      "cell_type": "code", 
      "source": [
        "from sklearn.ensemble import GradientBoostingRegressor"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 38, 
      "cell_type": "code", 
      "source": [
        "def deviance_plot(est, X_test, y_test, ax=None, label='', train_color='#2c7bb6', \n", 
        "                  test_color='#d7191c', alpha=1.0, ylim=(-80, 100)):\n", 
        "    \"\"\"Deviance plot for ``est``, use ``X_test`` and ``y_test`` for test error. \"\"\"\n", 
        "    n_estimators = len(est.estimators_)\n", 
        "    test_dev = np.empty(n_estimators)\n", 
        "\n", 
        "    for i, pred in enumerate(est.staged_predict(X_test)):\n", 
        "       test_dev[i] = est.loss_(y_test, pred)\n", 
        "\n", 
        "    if ax is None:\n", 
        "        fig = plt.figure(figsize = (12, 8))\n", 
        "        ax = plt.gca()\n", 
        "        \n", 
        "    ax.plot(np.arange(n_estimators) + 1, test_dev, color=test_color, label='Test %s' % label, \n", 
        "             linewidth=2, alpha=alpha)\n", 
        "    ax.plot(np.arange(n_estimators) + 1, est.train_score_, color=train_color, \n", 
        "             label='Train %s' % label, linewidth=2, alpha=alpha)\n", 
        "    ax.set_ylabel('Error')\n", 
        "    ax.set_xlabel('n_estimators')\n", 
        "    #ax.set_ylim(ylim)\n", 
        "    return test_dev, ax"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "execution_count": 39, 
      "cell_type": "code", 
      "source": [
        "gb = GradientBoostingRegressor(n_estimators=125, max_depth=5)\n", 
        "gb.fit(Xtraindf, ytrain)\n", 
        "deviance_plot(gb, Xtest, ytest, ylim=(0,0.8));\n", 
        "plt.legend();"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "Notice the wide gap. This is an indication of overfitting!\n", 
        "\n", 
        "Unlike random forests, where we are using the randomness to our benefits, the GBRT requires careful cross-validation\n", 
        "\n", 
        "Peter Prettenhofer, who wrote sklearns GBRT implementation writes in his pydata14 talk (worth watching!)\n", 
        "\n", 
        "Hyperparameter tuning I usually follow this recipe to tune the hyperparameters:\n", 
        "\n", 
        "- Pick n_estimators as large as (computationally) possible (e.g. 3000)\n", 
        "- Tune max_depth, learning_rate, min_samples_leaf, and max_features via grid search\n", 
        "- A lower learning_rate requires a higher number of n_estimators. Thus increase n_estimators even more and tune learning_rate again holding the other parameters fixed\n", 
        "- This last point is a tradeof between number of iterations or runtime against accuracy. And keep in mind that it might lead to overfitting.\n", 
        "\n", 
        "Let me add however, that poor learners do rather well. So you might want to not cross-validate max_depth. And min_samples_per_leaf is not independent either, so if you do use cross-val, you might just use one of those."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "source": [
        "EXERCISE: use a param grid and `GridSearchCV` to try and overfit this model as little as possible. We'll only explore part of the space as this takes time.\n", 
        "\n", 
        "eg:\n", 
        "\n", 
        "```\n", 
        "param_grid = {'learning_rate': [0.1, 0.01],\n", 
        "              'min_samples_leaf': [3, 5],  ## depends on the num of training examples, or use max_depth\n", 
        "              'max_features': [0.2, 0.6]\n", 
        "              }\n", 
        "```\n", 
        "\n", 
        "And remember you;ll have to figure the number of trees"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# your code here"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "Once you find the best model, make a deviance plot with the best params."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# your code here"
      ], 
      "outputs": [], 
      "metadata": {}
    }, 
    {
      "source": [
        "Also find the most informative features"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "source": [
        "### Note: what's the relationship between residuals and the gradient?\n", 
        "\n", 
        "Pavlos showed in class that for the squared loss, taking the gradient in the \"data point functional space\", ie a N-d space for N data points with each variable being $f(x_i)$ just gives us the residuals. It turns out that the gradient descent is a more general idea, and one can use this for different losses. And the upweighting of poorly fit points in AdaBoost is simply a weighing by gradient. If the gradient (or residual) is high it means you are far away from optimum in this functional space, and if you are at 0, you have a flat gradient!\n", 
        "\n", 
        "The ideas from the general theory of gradient descent tell us this: we can slow the learning by shrinking the predictions of each tree by some small number, which is called the learning_rate (learning_rate). This \"shrinkage\" helps us not overshoot, but for a finite number of iterations also simultaneously ensures we dont overfit by being in the neighboorhood of the minimum rather than just at it! But we might need to increase the iterations some to get into the minimum area."
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "source": [
        "---"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }
  ]
}