Implement XGBoost in R

XGBoost is a popular supervised machine learning algorithm that can be used for a wide variety of classification and prediction tasks. It can model linear and non-linear relationships and is highly interpretable as well.

XGBoost combines the strengths of multiple decision trees, guided by strategic optimization and regularization techniques, to deliver exceptional predictive performance and interpretability. These characteristics make it a popular choice for applications across various domains like finance, healthcare, e-commerce, recommendation and ranking systems, business intelligence, time-series forecasting, and demand forecasting.

XGBoost is a very popular algorithm in data science competitions (such as the Kaggle competitions) because of its flexibility, accuracy, speed, and the support of a thriving open source community.

XGBoost stands for Extreme Gradient Boosting. Gradient boosting is an approach where new models are created that predict the residuals or errors of prior models and are then added together in ensemble to make the final prediction. It is called gradient boosting because it uses a gradient descent algorithm to minimize the loss when adding new models. Tree-based algorithms such as decision trees, bagging, or random forest solve classification and regression problems using branching decision structures. These trees are built by randomly sampling the data set. Boosting algorithms like XGBoost use the gradient descent algorithm to build trees sequentially and minimize errors without overfitting the dataset.

A critical aspect of training machine learning models is preventing your model from overfitting the training data. Data scientists are always looking to find the right balance between minimizing train error and overfitting the training data so much that the model doesn’t generalize to new data. Regularization techniques, careful hyperparameter tuning, and monitoring both training and testing errors all help develop robust models with strong, predictive performance. The hyperparameters used to tune XGBoost are a particularly important part of training.

Key XGBoost hyperparameters include:

• The number of boosting iterations, referred to as nrounds
• The learning rate (also known as eta), which is a hyperparameter that scales the contribution of each tree in the ensemble
• The maximum depth of a tree, which is a pruning parameter designed to control the overall tree depth
• gamma, which is the minimum loss reduction required to make a further partition on a leaf node of the tree.
• subsample, which is the fraction of all observations to be randomly sampled for each tree
• lambda, which is the L2 regularization term
• alpha, which is the L1 regularization term

This tutorial will guide you through using XGBoost in R to train a model on the Pima-Indian-Diabetes data set from OpenML. This data set comprises 768 instances and traces its origin to the National Institute of Diabetes and Digestive and Kidney Diseases. The goal of the data set is to make predictions on whether a patient has diabetes using specific diagnostic measurements provided in the data set.

Let’s get started!

Prerequisites

You need an IBM Cloud account to create a watsonx.ai project.

Steps

Step 1. Set up your environment

While you can choose from several tools, this tutorial walks you through how to set up an IBM account to use a Jupyter Notebook. Jupyter Notebooks are widely used within data science to combine code, text, images, and data visualizations to formulate a well-formed analysis.

1. Log in to watsonx.ai using your IBM Cloud account.

2. Create a watsonx.ai project.

3. Create a Jupyter Notebook using the R Kernel.

This step will open a Notebook environment where you can load your data set and copy the code from this tutorial to implement a binary classification task using the gradient boosting algorithm.

Step 2. Install and load the required libraries

In this step, we’ll install the following libraries:

• data.table provides an enhanced version of data frames to perform efficient data manipulation
• xgboost is a library that provides an implementation of XGBoost
• caret used for simplifying the process of model training, tuning, and evaluation
• dpylr is a library for data manipulation
• ggplot2 is used to create a variety of static and dynamic plots, including scatter plots, bar plots, and heatmaps

Before moving forward, check in the console to be sure that all the required libraries have been installed correctly.

#Install and load the required libraries 
install.packages("data.table")
install.packages("xgboost")
install.packages("caret")
install.packages("dpylr")
install.packages("ggplot2")

#Load required libraries
library(data.table)
library(xgboost)
library(caret)
library(dpylr)
library(ggplot2)

Step 3. Load the data set

In this step, we will load the Pima-Indians-Diabetes data set that is available in the OpenML repository. The data set is accessible in the Attribute-Relation File Format (ARFF), a text file format commonly used to represent data for machine learning tasks. This format is often used for small to medium-sized data sets in educational settings. It is less common in large-scale or production environments, where other data interchange formats like CSV or Parquet are more widely used.

#Load the Pima-Indians-Diabetes data set
url <- "https://www.openml.org/data/download/4965316/pima.arff"
colnames <- c("Pregnancies", "Glucose", "BloodPressure", "SkinThickness", "Insulin", "BMI", "DiabetesPedigreeFunction", "Age", "Outcome")
diabetes_data <- fread(url, skip = 16, header = FALSE, col.names = colnames)

#Display the first few rows of the data set
head(diabetes_data)

The call to head() will show you the first 5 entries in the data set:

   Pregnancies Glucose BloodPressure SkinThickness Insulin   BMI DiabetesPedigreeFunction   Age Outcome
         <int>   <int>         <int>         <int>   <int> <num>                    <num> <int>   <int>
1:           4     117            62            12       0  29.7                    0.380    30       1
2:           4     158            78             0       0  32.9                    0.803    31       1
3:           2     118            80             0       0  42.9                    0.693    21       1
4:          13     129             0            30       0  39.9                    0.569    44       1
5:           5     162           104             0       0  37.7                    0.151    52       1
6:           7     114            64             0       0  27.4                    0.732    34       1

Our outcome is a binary categorical variable so we'll create a logistic regression model. XGBoost can also create linear models as well as multi-class classification models.

To learn more the parameters of the data set, visit Pima-Indian-Diabetes data set.

Step 4. Explore the data set

Before pre-processing data, you should conduct an exploratory data analysis (EDA) to understand the data’s structure and format. This can include the types of variables, their distributions, and the overall organization of information.

We begin by examining if there are any missing values in the data set.

#Check for missing data
missing_data <- colSums(is.na(diabetes_data))
missing_data

We see that there is no missing data in any of the columns in the data set.

Next, we’ll check how the attributes are statistically correlated by plotting a correlation heatmap. Since all the attributes are numeric and there are no categorical variables, we don’t have to numerically encode the data.

#Create a correlation heatmap for the data
correlation_matrix <- abs(cor(diabetes_data))

#Now plot the data
plt <- ggplot(data = as.data.frame(as.table(correlation_matrix)), aes(x = Var1, y = Var2, fill = Freq)) +
    geom_tile() +
    scale_fill_gradientn(colors = colorRampPalette(c("purple", "blue", "red"))(50), name = "Correlation") +
    geom_text(aes(label = round(Freq, 2)), color = "white", size=4) +
    theme_minimal() +
    theme(axis.text.x = element_text(size=10, angle = 45, hjust = 1), 
        axis.text.y = element_text(size = 10),
        plot.title = element_text(size = 10),   # Adjust plot title font size
        legend.text = element_text(size = 10),
        legend.title = element_text(size=10))
print(plt)

This creates a correlation heatmap graph:

If the graph appears too compressed you can re-run the code block to regenerate it.

From the correlation heatmap we can see that the attributes are not significantly correlated with each other. Hence, we’ll apply a threshold of 0.1 to identify the parameters that have an acceptable correlation with the variable Outcome. The parameters include Pregnancies, Glucose, Insulin, BMI, DiabetesPedigreeFunction, and Age. We’ll use these six features to train our model and use Outcome as the target variable. The Outcome variable uses binary encoding, which means that it displays a value of 1 (or True) for a diabetes-positive entry and 0 (or False) for a diabetes-negative entry. This type of encoding is referred to as ordinal encoding.

Another important step in EDA is whether the data set is balanced or imbalanced; that is, if one class has more instances than another. One way to find this is to visualize the histogram distribution of the data set for different classes.

#Show the class-wise distribution of the data set
class_counts <- table(diabetes_data$Outcome)
barplot(class_counts,
        names.arg = c("Tested Negative", "Tested Positive"),
        col = c("blue", "red"),
        xlab = "Class", ylab = "Number of Instances",
        main = "Class wise Distribution of Diabetes Dataset")

This creates the following:

We can see from the graph that the negative class has nearly twice as many instances as the positive class. This is an indication that the data set is imbalanced and it might lead to the model overfitting on the negative class. This doesn’t necessarily imply that the training data will overfit, it is good to be aware of this upfront in case we need to use over-sampling (upsampling) or under-sampling (downsampling) techniques.

Now that we’re confident in our data set and in what features we’ll be using, we’ll extract the required features and split the data set into a training set and a testing set. We’ll hold out 20% of the data for testing and use the remaining 80% for training.

#Extract required features
X_data <- diabetes_data[, c("Pregnancies", "Glucose", "BMI", "Age", "Insulin", "DiabetesPedigreeFunction"), with = FALSE]
y_data <- diabetes_data$Outcome

#Split the data set
set.seed(50)

# 80/20 split
split_indices <- createDataPartition(y_data, p = 0.8, list = FALSE)
X_train <- X_data[split_indices, ]
X_test <- X_data[-split_indices, ]
y_train <- y_data[split_indices]
y_test <- y_data[-split_indices]

Step 5. Train and evaluate the model with default hyperparameters

Now that we have a training data set and a test data set, we’re ready to train our model and evaluate performances. First, we’ll just use the default hyperparameters for XGBoost and see how our model performs on our data. The call to the xgboost training function expects, at a minimum, the following parameters:

• data: our training data
• labels: the labels that we’re predicting
• nrounds: the number of rounds to train
• objective: what we want our model to predict (In our case, this is a binary outcome, in that patients either have or do not have diabetes.)

Additionally we can pass in any hyperparameters that you want to set for the training. In this case, though, we we want to see how the default settings do first. We need to convert our X_train to a data.matrix and then pass in just the y_train list.

default_model <- xgboost(data = as.matrix(X_train),
                         label = y_train,
                         booster = "gbtree",
                         objective = "binary:logistic",
                         nrounds = 100,
                         verbose = 0)

Now that we’ve trained the model, we can check its performance by passing the model and our the testing data we split out to the predict function. The accuracy is simply calculated as the ratio of right responses.

y_pred <- predict(default_model, as.matrix(X_test), type = "response") > 0.5
accuracy <- sum(y_pred == y_test) / length(y_test)

So, we can see that our model has an accuracy of 0.6732026. To improve on this, we’ll look at using two techniques: a grid search of hyperparameters and cross validation.

Step 6. Use a grid search of hyperparameters and perform cross validation

We want to see if we can improve our models performance by tuning hyperparameters but it would be inefficient to try each different combination of values across all of the different parameters we could tune. Enter grid search.

Grid search is a process that searches exhaustively through a manually specified subset of the hyperparameter space of the targeted algorithm and evaluates the cost function based on the generated hyperparameter sets. If we want to test 4 different values for learning_rate and 4 different values for min_child_weight, then we’ll have 16 different models to train and compare. The grid search process runs each of these sets of hyperparameters for us and returns the set with the highest performance based on the relevant accuracy metric.

Cross validation is a way to estimate the performance of our model with less variance than we might get from using a single train-test set split. We split the data set into different chunks, often referred to as folds. The model is trained on all but one of the folds, so for k folds, we would train on k-1. The model is then tested on the held out fold. We then repeat this test so that each fold of the data set is the held back as the test fold for testing. After running this cross validation, we’ll have k performance scores that are summarized with a mean and a standard deviation value. This is a more time-consuming way of training data but it frequently can lead to better performance and is a common practice with XGBoost.

To perform both the grid search and cross validation, we could either roll our own or we could use the implementation in the caret library. First, we’ll create a list of all the hyperparameters that we’d like to search through:

hyperparam_grid <- expand.grid(
  nrounds = seq(from = 100, to = 300, by = 100),
  eta = c(0.025, 0.05, 0.1, 0.3),
  max_depth = c(4, 5, 6),
  gamma = c(0, 1, 2),
  colsample_bytree = c(0.5, 0.75, 1.0),
  min_child_weight = c(1, 3, 5),
  subsample = 1
)

Next, we’ll set up the configuration for the cross validation by calling the caret::trainControl method and creating a control object. We want the training process to generate cross validation for us, so we set the method parameter to cv and the number of folds to 4. The created object will be passed to train in the next step.

tune_control <- caret::trainControl(
  method = "cv", # cross-validation
  number = 4, # with n folds
  verboseIter = FALSE, # no training log
  allowParallel = FALSE
)

Finally, we’re ready to train our model. We pass the training data, our grid of parameters, and the control object to caret::train along with the type of model that we’re looking to find optimal parameters for. caret::train supports many different model types, in our case we want an XGBoost model, so we’ll pass xgbTree. We could additionally pass which evaluation metric to use, such as "RMSE" for a linear regression or "AUC" along with many other parameters for model tuning. There are also options for parallel computing if your data set is very large or if you have large numbers of hyperparameters to evaluate.

This might take a while to run since our model will be training and evaluating 1728 different models, but at the end we’ll have a very good picture of the best hyper-parameters.

bst <- caret::train(
  x = X_train,
  y = as.factor(y_train),
  trControl = tune_control,
  tuneGrid = hyperparam_grid,
  method = "xgbTree", #  this says we want XGB
  verbose = FALSE,
  verbosity = 0
)

Once tuning is complete, we can see the best hyperparameters:

bst$bestTune

This prints out the best values found in the grid search:

   nrounds max_depth   eta gamma colsample_bytree min_child_weight subsample
37     100         4 0.025     1             0.75                1         1

We can see that hyperparameter combination 37 had the best performing combination of values. Now that we’ve found the best combination of hyperparameters, we can pass those to the xgboost method for training. One additional hyperparameter used here is scale_pos_weight, which tells XGBoost to reweight the positive values in our training data to account for the imbalance in our data set.

final_model <- xgboost(data = as.matrix(X_train),
                       label = y_train,
                       booster = "gbtree",
                       objective = "binary:logistic",
                       nrounds = bst$bestTune$nrounds,
                       max_depth=bst$bestTune$max_depth,
                       colsample_bytree=bst$bestTune$colsample_bytree,
                       min_child_weight=bst$bestTune$min_child_weight,
                       subsample=bst$bestTune$subsample,
                       eta=bst$bestTune$eta,
                       gamma=bst$bestTune$gamma,
                       scale_pos_weight = 0.5, # because our dataset is unbalanced
                       verbose = 0)

# now evaluate
y_pred <- predict(final_model, as.matrix(X_test), type = "response") > 0.5
accuracy <- sum(y_pred == y_test) / length(y_test)

Our accuracy has increased to 0.7712418, which is approximately a 10% improvement. You might try using different sets of hyperparameters to see whether you can improve on this accuracy.

One other helpful feature of XGBoost is its interpretability. We can easily see how much different features affect the output of the model overall using the xgb.plot.importance method to graph the contributions of each feature to the model prediction.

importance_matrix <- xgb.importance(colnames(X_train), model = final_model)
xgb.plot.importance(importance_matrix, rel_to_first = TRUE, xlab = "Relative importance")

The XGBoost package also provides the ability to generate Shapley value plots and other diagnostic techniques, which are all possible because of the interpretable nature of XGBoost models.

Summary

In this tutorial, we introduced XGBoost and how it uses gradient boosting to combine the strengths of multiple decision trees for strong predictive performance and interpretability.

We’ve explored the Pima-Indians-Diabetes data set and identified the features to use to train an XGBoost binary classification model. We learned about cross validation and grid search and implemented them using the caret library. Once the best hyperparameters were identified, we trained an XGBoost model using the best hyperparameters and compared the performance to our default hyperparameter version.

Try watsonx for free

Build an AI strategy for your business on one collaborative AI and data platform called IBM watsonx, which brings together new generative AI capabilities, powered by foundation models, and traditional machine learning into a powerful platform spanning the AI lifecycle. With watsonx.ai, you can train, validate, tune, and deploy models with ease and build AI applications in a fraction of the time with a fraction of the data.

Try watsonx.ai, the next-generation studio for AI builders.

Next steps

Explore more articles and tutorials about watsonx on IBM Developer.