Lab 8 - Gradient Boosting using xgboost

This lab introduces workflows using xgboost, a well-known, powerful implementation of gradient boosting. Because xgboost is a new library that isn't included in the Anaconda distribution, you most likely will need to install it yourself using the command below:

# pip install xgboost

If your installation was successful, you should be able to import xgboost without any error messages:

import warnings   # Turns off the future warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
import xgboost as xgb

Additionally, you'll still need our standard libraries for this lab:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import sklearn

The data used in the lab's examples comes from a public health study conducted in Arahazar Upazila, Bangladesh where researchers identified community wells that were unsafe due to high levels of arsenic. They encouraged the families using these wells to switch to getting water from the nearest safe well. Several years after their initial visit, the researchers returned and logged which households they had asked to switch had actually done so, recorded as an outcome of switch = 'yes' in the dataset.

# Well switching data
wells = pd.read_csv("https://remiller1450.github.io/data/Wells.csv")

# Train-test split
from sklearn.model_selection import train_test_split
train, test = train_test_split(wells, test_size=0.1, random_state=9)

# Separate the outcome and predictors (and drop the "association" predictor)
train_y = (train['switch'] == 'yes').astype(int)
train_X = train.drop(['switch', 'association'], axis="columns")

Part 1 - Tuning the Model Structure

Similar to random forests, there are two categories of hyperparameters we are concerned with:

1. Those impacting the base learners directly (max depth, etc.)
2. Those impacting how the base learners are incorporated into the ensemble (namely the learning rate and number of iterations).

The most commonly tuned hyperparameters for the base learners, roughly in their order of importance are:

• max_depth - controls the complexity of individual trees by restricting their maximum depth.
• min_child_weight - controls how readily a tree is allowed to create additional splits by requiring a minimum amount of data (technically, weighted observations) within a node for it to be partitioned. Larger values produce simpler trees and can help reduce overfitting.
• subsample - promotes variety among trees and helps prevent overfitting by randomly selecting a fraction of the training observations to be used when building each tree.
• colsample_bytree - promotes variety among trees and helps prevent overfitting by randomly selecting a fraction of the available predictors that may be used when constructing each tree.

Less widely used, but potentially useful hyperparameters are reg_alpha and reg_lambda:

• reg_alpha (L1 regularization) encourages some tree contributions to become exactly zero. Larger values create a simpler model and can be useful when many predictors contain little useful information.
• reg_lambda (L2 regularization) shrinks tree contributions toward zero without forcing them to be exactly zero. Larger values generally produce a smoother, more conservative model that is less likely to overfit the training data.

Standard practice is to first find a configuration for the base learners that works adequately well, then to find a learning rate that is not too aggressive or too slow, then use early stopping to determine the effective number of boosting iterations.

The xgboost library is designed for compatibility with many functions from sklearn, including pipelines and cross-validated searches:

## Simple pipeline
from sklearn.model_selection import RandomizedSearchCV
import xgboost as xgb

## Simple xgboost model
model1 = xgb.XGBClassifier(eval_metric='auc')

## Randomized search space
parms = {
    "max_depth": [2, 3, 4, 5, 6],
    "min_child_weight": [1, 5, 10, 20],
    "subsample": [0.50, 0.75, 1.00],
    "colsample_bytree": [0.50, 0.75, 1.00]
}

## Randomized search
rand_res = RandomizedSearchCV(model1, param_distributions=parms, n_iter=20, scoring="roc_auc", cv=5, n_jobs=-1).fit(train_X, train_y)
print(rand_res.best_estimator_)
print(rand_res.best_score_)
print(rand_res.best_params_)

XGBClassifier(base_score=None, booster=None, callbacks=None,
              colsample_bylevel=None, colsample_bynode=None,
              colsample_bytree=1.0, device=None, early_stopping_rounds=None,
              enable_categorical=True, eval_metric='auc', feature_types=None,
              feature_weights=None, gamma=None, grow_policy=None,
              importance_type=None, interaction_constraints=None,
              learning_rate=None, max_bin=None, max_cat_threshold=None,
              max_cat_to_onehot=None, max_delta_step=None, max_depth=2,
              max_leaves=None, min_child_weight=20, missing=nan,
              monotone_constraints=None, multi_strategy=None, n_estimators=None,
              n_jobs=None, num_parallel_tree=None, ...)
0.6484246582797105
{'subsample': 1.0, 'min_child_weight': 20, 'max_depth': 2, 'colsample_bytree': 1.0}

Randomized search is a common strategy for xgboost due to the large number of available hyperparameters and diverse range of values we might seek to explore.

Additionally, since we are only tuning the structure of the base learners in this section, we will keep the ensemble-level parameters fixed at their defaults. After identifying a reasonable tree structure, we will return to the learning rate and use early stopping to determine the effective number of boosting iterations (shown in the next section).

Question #1:

• Part A: Inspect the randomized search results data frame created as part of this section's example, and identify the values of max_depth, min_child_weight, subsample and colsample_bytree that seem to be most often associated with lower performing models. Next, create an updated parameter search space that removes these values and considers two different values of reg_alpha and reg_lambda of either 0 (no shrinkage) or 10 (relatively high shrinkage). Run 30 iterations of randomized search and report the best hyperparameter configuration you identified.
• Part B: Compare the cross-validated ROC-AUC of the model you found in Part A (the best_estimator_) and a single decision tree with a maximum depth of 12. Briefly describe how these models compare in terms of their performance. What might be leading one model to outperform the other?

Part 2 - Learning Rates and Early Stopping

To properly determine the ensemble-level hyperparameters (learning rate and number of iterations) we'll need to split off part of the training data to use as a validation set. We can then have xgboost check the performance of the current ensemble on the validation set every fixed number of boosting iterations. The example below monitors the AUC across boosting iterations, stopping after the validation metric has failed to improve for 30 consecutive boosting iterations.

from sklearn.model_selection import train_test_split
import xgboost as xgb

## Create validation set
xgb_train_X, val_X, xgb_train_y, val_y = train_test_split(train_X,train_y,test_size=0.25,random_state=0)

## Start with best structure from Part 1
model2 = rand_res.best_estimator_

## Update parameters to include early stopping, and a trial learning rate with a large number of estimators
model2.set_params(early_stopping_rounds=30, learning_rate=0.03, n_estimators=5000, eval_metric="auc")

## Fit model
model2.fit(xgb_train_X, xgb_train_y, eval_set=[(xgb_train_X, xgb_train_y),(val_X, val_y)],verbose=False)
print(model2.best_iteration)
print(model2.best_score)

25
0.6589115136595938

In practice we'd check a few different learning rates to make sure we're not building the ensemble too aggressively or passively for the hyperparamter configuration we found in Part 1. Below we'll check a slightly smaller learning rate and slightly larger learning rate:

## Model w/ smaller eta
model2.set_params(early_stopping_rounds=30, learning_rate=0.01, n_estimators=5000, eval_metric="auc")
model2.fit(xgb_train_X, xgb_train_y, eval_set=[(xgb_train_X, xgb_train_y),(val_X, val_y)],verbose=False)
print(model2.best_iteration)
print(model2.best_score)

81
0.6589962450610513

## Model w/ larger eta
model2.set_params(early_stopping_rounds=30, learning_rate=0.08, n_estimators=5000, eval_metric="auc")
model2.fit(xgb_train_X, xgb_train_y, eval_set=[(xgb_train_X, xgb_train_y),(val_X, val_y)],verbose=False)
print(model2.best_iteration)
print(model2.best_score)

9
0.6564899793968908

We could continue this process, until we settle upon what we think is the best learning rate for our application.

Additionally, it can be useful to visualize the learning curve to ensure that early stopping worked as intended. Sometimes this curve might be noisy enough that we might decide to increase the value of early_stopping_rounds:

## Extract results from the model
xgb_results = model2.evals_result()

## Visualize by iteration
plt.figure(figsize=(7,4))
plt.plot(xgb_results["validation_0"]["auc"], label="Training")
plt.plot(xgb_results["validation_1"]["auc"], label="Validation")
plt.axvline(model2.best_iteration, color="gray", label="Best Iteration")
plt.xlabel("Number of boosting iterations")
plt.ylabel("AUC")
plt.legend()
plt.show()

Question #2:

• Part A: Considering the relationship between the learning_rate and how much the model changes after each boosting iteration. Should you favor larger or smaller values of early_stopping_rounds when the learning rate is small? Briefly explain your reasoning.
• Part B: Using the results from the expanded randomized search you performed in Part A, explore learning rates of 0.01, 0.1, and 0.30 using early_stopping_rounds = 40. Which learning rate appears to produce the strongest model?
• Part C: Create a plot of validation AUC versus boosting iteration for the three learning rates you identified in Part B. Briefly describe any relationships you see between the learning rate and the shape of the curve.

Part 3 - Evaluation Metrics

Our previous examples have used AUC as the evaluation metric without much additional consideration. However, from a model building perspective it is often useful to rely upon log-loss (also called cross-entropy loss) to tune the learning rate due to its greater sensitivity to small improvements in predicted probabilities.

$$ \text{log-loss} =-\frac{1}{n}\sum_{i=1}^{n} \left[y_i \log(\hat p_i) + (1-y_i)\log(1-\hat p_i) \right] $$

where:

• $y_i \in \{0,1\}$ is the true class label
• $\hat p_i$ is the predicted probability of the positive class

Unlike classification accuracy or AUC, log-loss rewards models for producing probabilities that are both correct and confident, while heavily penalizing confident mistakes. Learning curves that use log-loss tend to be smoother and more reliable, making it the most common evaluation metric for early stopping even if it is not the final evaluation metric for the application.

Consider the following learning curves:

## log-loss
model2.set_params(learning_rate=0.05,n_estimators=5000,early_stopping_rounds=30,eval_metric="logloss")
model2.fit(xgb_train_X, xgb_train_y,eval_set=[(val_X, val_y)],verbose=False)
logloss_curve = model2.evals_result()["validation_0"]["logloss"]


## AUC
model2.set_params(learning_rate=0.05,n_estimators=5000,early_stopping_rounds=30,eval_metric="auc")
model2.fit(xgb_train_X, xgb_train_y,eval_set=[(val_X, val_y)],verbose=False)
auc_curve = model2.evals_result()["validation_0"]["auc"]

## Compare learning curves
fig, ax = plt.subplots(1, 2, figsize=(10,4))

ax[0].plot(logloss_curve)
ax[0].set_title("Validation Log-Loss")
ax[0].set_xlabel("Iteration")

ax[1].plot(auc_curve)
ax[1].set_title("Validation AUC")
ax[1].set_xlabel("Iteration")
plt.show()

Question #3:

• Part A: Suppose an additional boosting iteration causes the predicted probability of a positive observation to increase from 0.7 to 0.8. Calculate how this observation's contribution to the log-loss changes.
• Part B: In the same scenario described in Part A, what must happen for the change in predicted probability to improve the AUC? Hint: As we briefly discussed after our last quiz, changes in AUC depend upon the relative ordering of observations by predicted probability. On the quiz we saw two models that ordered the observations in the same way, and therefore had the same AUC despite assigning different predicted probabilities to each observation.
• Part C: Based upon the learning curves shown in this example and your responses to Parts A and B, explain why a practitioner might prefer log-loss when monitoring early stopping, even if AUC is going to be the final evaluation metric calculated on the test set?
• Part D: Based upon the learning curves shown in this example and your responses to Parts A, B, and C, would you recommend setting early_stopping_rounds to a higher or lower value when monitoring log-loss in comparison with monitoring AUC? Briefly explain.

Part 4 - Application

In this section you'll revisit the Wisconsin Breast Cancer data set. Recall that the goal in this application was to accurately classify tissue samples from tumors as either benign or malignant.

## Read the data
wbc = pd.read_csv("https://remiller1450.github.io/data/wisc_bc.csv")

The purpose of this section is to provide an opportunity for you to practice the xgboost workflow more independently.

Question #4:

• Part A: Split the data into three segments, a training set consisting of 80% of the observations and a test set consisting of the remaining 20%. Then, split the training set once more so that 25% of it is used for validation.
• Part B: Use either a grid search or randomized search to find a reasonable combination of max_depth, min_child_weight, subsample, colsample_bytree, reg_alpha and reg_lambda for the default learning rate and number of iterations. Your search should consider at least a few sensibly chosen values of each hyperparameter, noting that boosting tends to work best with relatively shallow trees.
• Part C: Use early stopping and your validation set to find an optimal learning rate and number of boosting iterations using the logloss metric. Display the learning curve showing training and validation performance by boosting iteration for the learning rate you selected. You should try at least three different learning rates and provide a brief justification for your choice of early_stopping_rounds.
• Part D: Calculate the classification accuracy and the ROC-AUC of the best xgboost model (from Part C) on the test set.
• Part E: Develop an effective random forest model using the training data and the guidelines from Lab 6. Calculate the classification accuracy and the ROC-AUC on the test set and briefly comment on how it compares to your xgboost model results.
• Part F: Briefly explain why a validation set is important for the approach we used to develop our xgboost model. That is, why didn't we just use the test set to monitor performance?