Building model using GBM

Let’s first just build a model using the gradient boosting machine (GBM) with arbitrarily chosen hyper-parameter. GBM is a machine learning algorithm for regression and classification problems. It builds a model in the form of an ensemble of (weak) prediction models like decision trees. For more details refer to this, this and this.

gbm.1 <- h2o.gbm(  
  x                   = x,
  y                   = y,
  training_frame      = train_frame,
  validation_frame    = valid_frame,
  ntrees              = 100,
  max_depth           = 3,
  learn_rate          = 0.2)

In the above example we used train_frame to train the model and valid_frame to evaluate the performance. We can see the performance in tho following:

 h2o.performance(gbm.1, valid = T)

## H2OBinomialMetrics: gbm
## ** Reported on validation data. **
## 
## MSE:  0.05118193
## RMSE:  0.2262342
## LogLoss:  0.2170323
## Mean Per-Class Error:  0.126074
## AUC:  0.8941919
## Gini:  0.7883838
## 
## Confusion Matrix for F1-optimal threshold:
##        False. True.    Error     Rate
## False.    689    24 0.033661  =24/713
## True.      26    93 0.218487  =26/119
## Totals    715   117 0.060096  =50/832
## 
## Maximum Metrics: Maximum metrics at their respective thresholds
##                         metric threshold    value idx
## 1                       max f1  0.226842 0.788136 113
## 2                       max f2  0.156576 0.789902 132
## 3                 max f0point5  0.530581 0.843373  70
## 4                 max accuracy  0.226842 0.939904 113
## 5                max precision  0.999967 1.000000   0
## 6                   max recall  0.000327 1.000000 399
## 7              max specificity  0.999967 1.000000   0
## 8             max absolute_mcc  0.226842 0.753161 113
## 9   max min_per_class_accuracy  0.068766 0.857143 178
## 10 max mean_per_class_accuracy  0.156576 0.878811 132
## 
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`

Each algorithm has number of hyper-parameters to be set, like in case of GBM number of trees or learning rate, etc. As always one of the important question is which values to choose for the hyper-parameters. There are different ways to find the optimal values for hyper-parameters, like grid search, random search or other general optimization algorithm.

We can use Bayesian optimization to find out an approximation of the optimized hyper-parameters.

Bayesian Optimization

H2O supports grid search and random search. Grid search is kind of brute force (it increases exponentially with the number of parameters) and random search as its name suggests searches the hyper-parameter domain randomly.

Bayesian optimization (BO) is used for automatic tuning of hyper-parameters. In this approach the information obtained from previous function evaluations (i.e. building and evaluating models using various hyper-parameters) is used to approximate best next hyper-parameters to be tried out. BO builds a probabilistic model for the objective function, which here is the performance of the model. For more information have a look here and here.

I wished H2O guys have implemented Bayesian optimization. But no worries. There are bunch of Bayesian optimization algorithms one can use, in R, Python, C++, etc. On the other hand the heavy lifting part, which is function evaluation (i.e. building the models), will be carried out by H2O.

On the other hand there are some discussions whether BO worth the effort, e.g. see here. Any ways I think at least it worth mentioning how we can combine BO with H2O.

In the following it shall be shown how we can use rBaysianOptimisation package in R for finding the optimal hyper-parameters.

The function h2o_bayes is defined as a wrapper for h2o.gbm function, which builds and evaluates a GBM model. All the parameters that needs to be optimized are defined as an argument for this function. In order to score each model validation data set is used. However, for small data set (as in this example) usually it is better to use cross-validation (instead of validation data set) to decrease variation in the model score.

#============================
# Define the wrapper function
#============================
h2o_bayes <- function(
  max_depth, learn_rate, sample_rate, 
  col_sample_rate, balance_classes){
  bal.cl <- as.logical(balance_classes)
  gbm <- h2o.gbm(  
    x                   = x,
    y                   = y,
    training_frame      = train_frame,
    validation_frame    = valid_frame,
    #nfolds              = 3,
    ntrees              = 900,
    max_depth           = max_depth,
    learn_rate          = learn_rate,
    sample_rate         = sample_rate,
    col_sample_rate     = col_sample_rate,
    score_tree_interval = 5,
    stopping_rounds     = 2,
    stopping_metric     = "logloss",
    stopping_tolerance  = 0.005,
    balance_classes     = bal.cl)
    
  score <- h2o.auc(gbm, valid = T)
  list(Score = score,
       Pred  = 0)
}

#============================
# Find optimal values for the 
# parameters in the given range. 
#============================
OPT_Res <- BayesianOptimization(
  h2o_bayes,
  bounds = list(
    max_depth   = c(2L, 8L), 
    learn_rate  = c(1e-4, 0.2),
    sample_rate = c(0.4, 1), 
    col_sample_rate = c(0.4, 1), 
    balance_classes = c(0L, 1L)),
  init_points = 10,  n_iter = 10,
  acq = "ucb", kappa = 2.576, eps = 0.0,
  verbose = FALSE)

## 
##  Best Parameters Found: 
## Round = 10   max_depth = 6.0000  learn_rate = 0.1223 sample_rate = 0.7886    col_sample_rate = 0.4054    balance_classes = 0.0000    Value = 0.9001

Building the model using optimal parameters

After searching for the best hyper-parameters using Bayesian optimization algorithm, we use the values that were found above to train a GBM model.

gbm <- h2o.gbm(
  x                   = x,
  y                   = y,
  training_frame      = train_frame,
  validation_frame    = valid_frame,
  ntrees              = 900,
  max_depth           = OPT_Res$Best_Par["max_depth"],
  learn_rate          = OPT_Res$Best_Par["learn_rate"],
  sample_rate         = OPT_Res$Best_Par["sample_rate"],
  col_sample_rate     = OPT_Res$Best_Par["col_sample_rate"],
  balance_classes     = as.logical(OPT_Res$Best_Par["balance_classes"]),
  score_tree_interval = 5,
  stopping_rounds     = 2,
  stopping_metric     = "logloss",
  stopping_tolerance  = 0.005,
  model_id         = "my_awesome_GBM")

Train a model using only important co-variables

In the previous section we used all the columns to train a GBM model. Many times removing irrelevant data from the training data will improve the results, for various reason. Even if it doesn’t improve the result, it will reduce the amount of resources needed for calculation, e.g. CPU time and memory.

Now we remove the variables that are not important. However, note that the variable importance is obtained from previous calculation.

var.imp <- h2o.varimp(gbm)[h2o.varimp(gbm)$scaled_importance > 0.05, "variable"]
# The value of 0.05 is arbitrary, you might want to use other values.

setdiff(x, var.imp)

## [1] "Account Length" "Area Code"      "Phone"          "Day Calls"     
## [5] "Eve Calls"      "Night Calls"

gbm_varImp <- h2o.gbm(
  x                   = var.imp,
  y                   = y,
  training_frame      = train_frame,
  validation_frame    = valid_frame,
  ntrees              = 900,
  max_depth           = OPT_Res$Best_Par["max_depth"],
  learn_rate          = OPT_Res$Best_Par["learn_rate"],
  sample_rate         = OPT_Res$Best_Par["sample_rate"],
  col_sample_rate     = OPT_Res$Best_Par["col_sample_rate"],
  balance_classes     = as.logical(OPT_Res$Best_Par["balance_classes"]),
  score_tree_interval = 5,
  stopping_rounds     = 2,
  stopping_metric     = "logloss",
  stopping_tolerance  = 0.005,
  model_id         = "my_awesome_GBM_varImp")

Compare the two model

The two models are compared in the plots below, in terms of rmse and auc. It is can be seen that if we use only the important variables the performance of the model is not reduced significantly, if any.

shist <- gbm@model$scoring_history[, c("duration", "validation_rmse", "validation_auc")]
shist$algorithm <- "GBM" 
scoring_history <- shist

shist <- gbm_varImp@model$scoring_history[, c("duration", "validation_rmse", "validation_auc")]
shist$algorithm <- "GBM with var.imp." 
scoring_history <- rbind(scoring_history,shist)

scoring_history$duration <- as.numeric(
  gsub("sec", "", scoring_history$duration))

scoring_history <- melt(scoring_history, id = c("duration", "algorithm"))

ggplot(data = scoring_history, 
       aes(x     = duration, 
           y     = value, 
           color = algorithm,
           group = algorithm)) + 
  geom_line() + geom_point() +
  facet_grid(. ~ variable, scales = "free",shrink = TRUE,space = "free")

Choose the best model based on AUC measure (you could choose mse, logloss or other measures which you might find more appropriate in your case).

AUC_gbm        <- h2o.performance(gbm, valid = T)@metrics$AUC
AUC_gbm_varImp <- h2o.performance(gbm_varImp, valid = T)@metrics$AUC
if(AUC_gbm > AUC_gbm_varImp){
  bestModel <- gbm
}else{
  bestModel <- gbm_varImp
}
cat("The best model is '", bestModel@model_id, 
    "' with AUC of ", max(AUC_gbm, AUC_gbm_varImp), 
    " vs ",  min(AUC_gbm, AUC_gbm_varImp), "\n" )

## The best model is ' my_awesome_GBM_varImp ' with AUC of  0.8956946  vs  0.8955178

Variable Importance

In many situations not all the co-variables are equally important. Some variable are more important than the others. Using variable importance we can focus more on important variables.

h2o.varimp_plot(bestModel)

Variable Importance per Customer

The above mentioned variables were measured based on the whole training data set and is only valid on average. However the importance of each variable might be different for each individual customer.

We can use the trained model to approximate the variable importance based on each individual customer.

In this approach a prediction is done on the data set. And then we leave a column out and will do the prediction again and measure the difference in probabilities. The difference in probabilities can indicate the importance of that variable for each customer.

This have a little bit of calculation in R side, like calculating differences and ranking each row based on these differences. Hence, in order to speed up things I have implemented a C++ code (get_high_rank_values.cpp) which also needed to be sourced. Using Rcpp package one can simply use C++ code in R script.

require(Rcpp)
Rcpp::sourceCpp("get_high_rank_values.cpp")
source("rowVarImp.R")

df <- rowVarImp.h2o(model=bestModel, cols=var.imp, data_frame=test_frame, n=3)

Here only first 3 important variables for each customer are calculated. The first 3 columns show the differences in the probabilities by not including the variable. The next 3 columns show their values and the last 3 columns the variables name.

Knowing these information the business unit can take measures based on each individual customer (not averaged on all customer like the general variable importance).

head(df)

Partial Dependence Plots

In the previous section we talked about importance of each variable, but we didn’t answer the question of how changing variables would change the outcome.
However, it is very important to know if we change value of a variable what will be the impact on the outcome.

This kind of question can be answered by the so called partial dependence plots. Partial dependence plots depicts, indeed, the marginal effect of the variables. In the following the partial dependence plots are shown for the 5 top important variables.

In order to calculate partial dependence we can use partialDependencePlot.R. However, since at least version 3.11.0.3717 h2o has h2o.partialPlot function to calculate one dimensional partial dependence. I like my own implementation better because of its functionality, however h2o.partialPlot is much faster.

cols <- var.imp[1:5]
pdps <- h2o.partialPlot(bestModel, data_frame, cols, nbins = 100, plot = F)
mean.orig <- mean(h2o.predict(bestModel, data_frame)[,3])

lapply(1:length(cols), function(x){
  x_ <- pdps[[x]] 
  colnames(x_) <- c("variable", "diff")
  x_$diff <- x_$diff - mean.orig
  if(class(x_$variable) == 'factor' || class(x_$variable) == 'character'){
    # order based on diff
    x_ <- within(x_, variable <- factor(variable, levels=variable[order(-diff)]))
    p <- ggplot(data = x_, aes(x = variable, y = diff)) +
      geom_bar(stat="identity" )
  }else{
    p <- ggplot(data = x_, aes(x = variable, y = diff)) +
      geom_line()
  }
  p + ggtitle(cols[x])
})

The above plots show only impact of single variable if others are kept constant. But in many cases, specially if there is interaction between variables, we are interested to know the impact on changing two or more variables. For this also we can use partialDependencePlot.R. As an example, in the following we can see the impact of “Day Charge” for each state (the colored thin lines), and also on average (the thick black line). It can be seen that there are quite variances among different states.

source("partialDependencePlot.R")
cols <- c("State", "Day Charge")
pdp1 <- partialDependencePlot(bestModel, cols, data_frame)

pdp2 <- partialDependencePlot(bestModel, "Day Charge", data_frame)

pdp1 <- melt(pdp1)
pdp2 <- melt(pdp2)
pdp2$State <- "AVG"

colnames(pdp1) <- c("State", "Day.Charge", "diff")
colnames(pdp2) <- c("Day.Charge", "diff", "State")

ggplot(data = pdp1, 
       aes(x     = Day.Charge, 
           y     = diff,   
           color = State,
           group = State)) + 
  geom_line() +
  geom_line(data = pdp2, size=3, color='black')

A Simplified Single Decision Tree

Alright, so far we did quite a lot. We found important variables, impact of changing each variables on the outcome. However, we would like to have the big picture. A high level insight about the model.

For that, we can build a single decision tree. But we already know that if we use single decision tree to build a model it wouldn’t be that accurate.

In order to improve the accuracy of the single decision tree a little bit, we use the decision tree to approximate the GBM model (instead of training it on the original data). We train a decision tree with that data but the target variable is replaced by the predictions of GBM. The premise behind this, according to what a friend of me Armin Monecke told once, is that using predictions of GBM (or other complex models) kind of removes the noise from the data, and makes it easier for less complex model (like single decision trees) to train.

Additionally, here, we use only rows that are predicted with high confidence, i.e. probabilities close to 0 and 1.

Note: since we sample the data points, we need to take care of prior distribution (e.g. sampling with the same distribution of target column or specifying it for example using a prior parameter in the algorithm).

denoise.h2o.df <- function(model, data_frame, destination_frame, lower.q=0.2, higher.q=0.8){
  # Removes the data points that cannot be classified significantly.
  # That means the data points for which the prediction is far away from 0 and 1, 
  # like those close to 0.5
  #
  # Args:
  #   model: A H2O model.
  #   data_frame: A H2O data frame.
  #   destination_frame: ID for the result frame.
  #   lower.q: The lower quantile, from which the data points will be ignored.
  #   higher.q: The higher quantile, to which the data points will be ignored.
  #
  # Returns:
  #   The resulting H2O data frame.
  
  pred   <- h2o.predict(model, data_frame)
  predCol <- names(pred)[3]
  quant  <- h2o.quantile(pred[,predCol], probs = c(lower.q, higher.q))
  lower  <- min(quant)
  higher <- max(quant)
  data_pred <- h2o.cbind(data_frame, pred)
  data_denoised <- data_pred[
    data_pred[,predCol] < lower | 
      data_pred[,predCol] > higher,]
  data_denoised <- h2o.assign(data_denoised, destination_frame) 
  return(data_denoised)  
}

denoised_train_df <- 
  denoise.h2o.df(
    model             = bestModel, 
    data_frame        = train_frame, 
    destination_frame = "denoised_train_frame", 
    lower.q           = 0.1,
    higher.q          = 0.9)

denoised_test_df <- 
  denoise.h2o.df(
    model             = bestModel, 
    data_frame        = test_frame, 
    destination_frame = "denoised_test_frame", 
    lower.q           = 0.1,
    higher.q          = 0.9)

train_df <- as.data.frame(denoised_train_df)
test_df  <- as.data.frame(test_frame)
test_df1 <- as.data.frame(denoised_test_df)

# Ignore the columns that is not in used in training 
tmp.df <- train_df[, -which(names(train_df)  %in% c("Churn.","False.","True.", "Phone" ))]
colnames(tmp.df)[colnames(tmp.df) == "predict"] <- "Churn."

counts <- as.data.frame(h2o.table(data_frame[,"Churn?"]))

priorDist <- counts$Count[1] / sum(counts$Count)
priorDist <- c(priorDist, 1 - priorDist)
               
tree <- rpart(
    formula  = Churn. ~ ., 
    data     = tmp.df,
    maxdepth = 5, 
    cp       = 0.02,
    parms    = list(prior = priorDist, 
                    split = "information"))

fancyRpartPlot(tree)

id <- which(!(test_df$State %in% levels(train_df$State)))
test_df$State[id] <- NA

Pred1 <- predict(tree, test_df)
test_df$Prediction1 <- Pred1[,"True."]
ROC1 <- roc(test_df$Churn., test_df$Prediction1)
ROC1$auc

## Area under the curve: 0.6869