<jupyter_start><jupyter_text>Gradient Boosting for Prediction and Inference Lesson 2In this notebook, we will walk through the details of how Gradient Boosting Works, and demonstrate their use via the scikit-learn `GradientBoostingClassifier` and `GradientBoostingRegressor`.  We will also begin to discuss how to approach setting parameters and using early stopping.<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Gradient Boosting: History and Terminology- Generally, when people refer to "Gradient Boosting" what they really mean is "Gradient Boosted Decision Trees"- "Boosting" in general, originally referred to a process where you fit a model, and then gave more weight to the examples in the training data that it got wrong, refit the model, and then added the new model to the previous one.  - AdaBoost (Freund and Schapire) was an early example of this, and using trees as the base learner was shown to be quite effective- Later, it was shown that AdaBoost was equivalent to fitting the base learner to the *gradient* of a particular loss function (exponential). - The idea was then generalized to any particular loss function, and referred to as *gradient boosting*.  Additionally, other loss functions were shown to be better in practice than the exponential loss.- Though the idea can apply to any base learner, decision trees have been most effective.  So, what we are learning about today is *gradient boosting with decision trees as the base learner* or *gradient boosted decision trees* (but usually just called *gradient boosting*)  Gradient boosting: rough pass- Consider a regression problem.  Suppose you have a model $M_1$ that already performs reasonably well.  But you would like to improve it if possible.- You use your model to make predictions on the training data and measure the *residuals*: i.e. the difference between the true answer and the prediction.- Next, you build a decision tree to try and predict the *residuals* given your predictors.  For example you want to find a tree that can distinguish between the cases where your previous model underpredicted vs overpredicted (and also distinguish between a small under/overprediction and a large under/overprediction)    - In theory, if your new tree does a great job at predicting the residuals, you can add the result of the tree to the previous model, and voila, you will have a better model.- In reality we need to be concerned with overfitting: Suppose the residuals are truly random noise.  With enough predictors, you may find an incidental pattern that works on your training data, but is worthless on new cases.- Therefore, we do a few things to try to prevent this from happening (i.e. to help *regularize*)- Most importantly, rather than directly add the new model to the old model, we first multiply it by a small fraction (say .1 or even .01 or .001).  We will do this process many times, so it doesn't (generally) hurt to take small steps.- It also doesn't matter whether we start with a "pretty good" model.  We can start with a "nothing" model (predict 0 for everything) and still iterate through this process.  Wait! This makes sense, but I still don't get a few things...- Where does the loss function come in?- Where does the gradient come in?- What if I am doing a classification problem?  Good questions!- What I just described was specific to regression.- In fact, it was specific to regression with a mean squared error loss function- When using mean squared error, the derivative of the loss function is the residual (distance from the true answer) (times a constant of 2)- $\frac{\partial}{\partial y} (y -y_t)^2 = 2(y-y_t)$- So fitting a model to the residuals is *equivalent* to fitting a model to the *gradient of the loss function*- So the notion of predicting small/large under/overpredictions can be thought of as predicting where my loss is increasing/decreasing by a small/large amount.- The latter notion will even generalize to cases where the loss function has a weirder form.  Reiteration of Gradient Boosting- In general, (for classification, regression with any loss function) we take the current model, use it to make predictions on the training data, and then find the *derivative of the loss function* for each point *at its current predicted value*.- For most loss functions, the gradient is "steeper" the further you are from the right answer, so the idea of predicting the "degree" of over/underprediction is still the right one.  Example of a regression loss function<jupyter_code>tvec1 = np.linspace(-20,20,201)
plt.figure(figsize=(8,4))
plt.subplot(1,2,1)
plt.plot(tvec1, (tvec1-2)**2)
plt.title('Loss fn when true value = 2')

plt.subplot(1,2,2)
plt.plot(tvec1, (tvec1-(-7))**2)
plt.title('Loss fn when true value = -7');<jupyter_output><empty_output><jupyter_text>Example of classification (cross-entropy) loss function<jupyter_code>tvec = np.linspace(.01,.99,99)
plt.figure(figsize=(8,4))
plt.subplot(1,2,1)
plt.plot(tvec, np.log(1/tvec))
plt.title('Loss fn when true value = 1')

plt.subplot(1,2,2)
plt.plot(tvec, np.log(1/(1-tvec)))
plt.title('Loss fn when true value = 0');<jupyter_output><empty_output><jupyter_text>Detailed steps- Let $M_0$ be the current model- Let $y_0$ be the current model's predictions on the training set    - These will be point estimates for regression    - These will be predicted probabilities for classification- Let $y_t$ be the true answers (numbers for regression, 0/1 for classification)- Let $L(y_{true}, y_{pred})$ be the loss function.- Let $L^{\prime}_{y_{true}}(y_{pred}) = \frac{\partial}{\partial y_{pred}}(L(y_{true}, y_{pred}))$- Let $z = L^{\prime}_{y_t}(y_1)$ -- i.e. the derivative of the loss function (with respect to the prediction) evaluated at each pair of (true, predicted) values.- Finally, fit a decision tree that tries to predict $z$ using your predictors $X$.    - Call the resulting tree $T_1$.- Let $M_1 = M_0 + \epsilon T_1$ (epsilon is the "learning rate")- Repeat with $M_1$ as the current model to build $T_2$  (etc. etc.) Alternate interpretation: - We are doing "gradient descent" by building trees  Let's play with some dataWe'll return to the games from the NBA, predicting the winner from team statistics (not including points).<jupyter_code>from sklearn.ensemble import GradientBoostingClassifier, GradientBoostingRegressor
from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
from sklearn.metrics import log_loss, roc_auc_score, accuracy_score
from sklearn.model_selection import train_test_split
df_nba = pd.read_csv('data/games.csv')
df_nba.dropna(inplace=True)
df_nba.shape
df_nba.sample(5)
df_nba.columns
feat_1 = ['FG_PCT_home','FT_PCT_home', 'FG_PCT_away', 'FT_PCT_away']
feat_2 = ['FG_PCT_home','FT_PCT_home', 'FG3_PCT_home', 'AST_home', 'REB_home', 
        'FG_PCT_away', 'FT_PCT_away', 'FG3_PCT_away', 'AST_away', 'REB_away']

X = df_nba.iloc[:,:-1]  # everything except winner
y = df_nba.HOME_TEAM_WINS 
X_train_full, X_test, y_train_full, y_test = train_test_split(X,y,test_size = 2000, random_state=0)
X_train_full.columns
n_train_pts = 2000
X_train_1 = X_train_full.iloc[:n_train_pts].loc[:, feat_1]
y_train_1 = y_train_full.iloc[:n_train_pts]
X_test_1 = X_test.loc[:, feat_1]
## For reference, lets redo the random forest
rf1 = RandomForestClassifier(n_estimators=1000)
rf1
rf1.fit(X_train_1, y_train_1)
prob_preds_rf1 = rf1.predict_proba(X_test_1)[:,1]
log_loss(y_test, prob_preds_rf1), roc_auc_score(y_test, prob_preds_rf1)
hard_preds_rf1 = (prob_preds_rf1 > .5).astype(int)
accuracy_score(y_test, hard_preds_rf1)
## now let's try a gb model
gb1 = GradientBoostingClassifier(max_depth=2, 
                                 n_estimators=100, learning_rate=.1)
gb1
gb1.fit(X_train_1, y_train_1)
prob_preds_gb1 = gb1.predict_proba(X_test_1)[:,1]
log_loss(y_test, prob_preds_gb1), roc_auc_score(y_test, prob_preds_gb1)
hard_preds_gb1 = (prob_preds_gb1 > .5).astype(int)
accuracy_score(y_test, hard_preds_gb1)<jupyter_output><empty_output><jupyter_text>Now let's repeat with the bigger feature set<jupyter_code>X_train_2 = X_train_full.iloc[:n_train_pts].loc[:, feat_2]
y_train_2 = y_train_full.iloc[:n_train_pts]
X_test_2 = X_test.loc[:, feat_2]
rf2 = RandomForestClassifier(n_estimators=1000)
rf2.fit(X_train_2, y_train_2)
prob_preds_rf2 = rf2.predict_proba(X_test_2)[:,1]
log_loss(y_test, prob_preds_rf2), roc_auc_score(y_test, prob_preds_rf2)
hard_preds_rf2 = (prob_preds_rf2>.5).astype(int)
accuracy_score(y_test, hard_preds_rf2)
gb2 = GradientBoostingClassifier(max_depth=2, n_estimators=100, learning_rate=.1)
gb2.fit(X_train_2, y_train_2)
prob_preds_gb2 = gb2.predict_proba(X_test_2)[:,1]
log_loss(y_test, prob_preds_gb2), roc_auc_score(y_test, prob_preds_gb2)
hard_preds_gb2 = (prob_preds_gb2>.5).astype(int)
accuracy_score(y_test, hard_preds_gb2)<jupyter_output><empty_output><jupyter_text>Gradient Boosting vs Random Forest- Typically, Gradient Boosting can outperform Random Forests by a small, but not insignificant amount.- However, Gradient Boosting requires much more "parameter tuning" to get the best performance- For Random Forest, can usually use 1K or 2K trees and do well.  The only "major" parameter is the `max_features` and the default is usually reasonable.- For Gradient Boosting, parameter tuning is considerably more complicated.- Random Forests tend to be less "well-calibrated" than GB.  That means that the log_loss of random forests might be bad even when the AUROC or Accuracy are still reasonably good.  Major parameters to tune for Gradient Boosting- `max_depth`: How deep to build the trees.  This is a very important parameter and performance may change dramatically for different values.  Larger values are more likely to overfit, smaller values more likely to underfit.  I've worked on problems where max_depth was 10 or even higher for best performance.- `learning_rate`, `n_estimators`: These parameters are also very important and interact highly with one another (and with `max_depth`).  Typically, the smaller your "step_size" (learning_rate), the more steps you will need to take to reach maximum performance. However, unlike random forests, if you continue to build trees in boosting, you will start overfitting, and performance (measured on the test set) will get worse.The best way to deal with `learning_rate` and `n_estimators` is to use a "validation set".  That is, take some data (that is not in the training set) and evaluate the performance of the model periodically as it trains.  When performance starts getting worse on the validation set, then stop.In theory, you should not use the test set as your validation set, because you will be overstating your performance a bit.  (Since you will be choosing the moment of peak performance on that test set to stop).  In practice, this is usually not such a big deal.  It is also fine just to set aside some data out of the training set.  The validation set does not need to be huge, and it is not crucial to stop at the exact perfect moment -- you just want to be in the right neighborhood.  Sklearn early stoppingin version 0.20 (Sep 2018), sklearn added support for early stopping via `validation_fraction` and `n_iter_no_change`.  Unfortunately, it does not output the loss values on the validation set, so it is hard to truly see it in action.  Also, it does not support passing a separate validation set.  We will revisit early stopping when we explore the other boosting packages.- `validation_fraction`: fraction of data to use for validation- `n_iter_no_change`: Compare current iteration to the value $x$ trees ago to decide whether to stop or not<jupyter_code># Let's run with early stopping
# We can make the number of trees large and just rely on early stopping
# We can also make the learning rate small to be on the safe side
gb3 = GradientBoostingClassifier(max_depth=2, 
                                 n_estimators=5000, 
                                 learning_rate=.01,
                                validation_fraction=.1,
                                n_iter_no_change=20,
                                verbose=2)
gb3.fit(X_train_2, y_train_2)
prob_preds_gb3 = gb3.predict_proba(X_test_2)[:,1]
log_loss(y_test, prob_preds_gb3), roc_auc_score(y_test, prob_preds_gb3)<jupyter_output><empty_output><jupyter_text>NOTE: When comparing boosting models, I put more credence in the log_loss values then the ROC or Accuracy scores<jupyter_code>hard_preds_gb3 = (prob_preds_gb3>.5).astype(int)
accuracy_score(y_test, hard_preds_gb3)<jupyter_output><empty_output>