Predicting large and imbalanced data set using the R package tidymodels

Summary of the variables

First let’s check the dimension of the two sets to be aware of what we are dealing with.

dim(train)

[1] 60000   171

dim(test)

[1] 16000   171

The training set has 60000 rows and 171 variables which is moderately large dataset. Inspecting thus this data by the usual functions such as summary, str would give heavy and not easily readable outputs. the best alternative However is by extracting the most important information that is required for building any machine learning model in aggregated way, for instance, what type of variables it has, some statistics about the variable values, the missing values..etc.

map_chr(train, typeof) %>% 
  tibble() %>% 
  table()

.
character    double 
      170         1 

Strangely, all the variables but one are characters, which is in contradiction with the description of this data from the file description. To figure out what is going on we display some few rows and some few columns.

train[1:5,1:7]

# A tibble: 5 x 7
  class aa_000 ab_000 ac_000     ad_000 ae_000 af_000
  <chr>  <dbl> <chr>  <chr>      <chr>  <chr>  <chr> 
1 neg    76698 na     2130706438 280    0      0     
2 neg    33058 na     0          na     0      0     
3 neg    41040 na     228        100    0      0     
4 neg       12 0      70         66     0      10    
5 neg    60874 na     1368       458    0      0     

I think, the problem is that the missing values in the data indicated by na are not recognized as missing values, instead they are treated as characters and this what makes the function read_csv coerces every variable that has this na values to character type. To fix this problem we can either go back and set the na argument to “na”, or we set the missing values by hand as follows.

train[-1] <- train[-1] %>% 
  modify(~replace(., .=="na", NA)) %>%
  modify(., as.double)

Now let’s check again

map_chr(train, typeof) %>% 
  tibble() %>% 
  table()

.
character    double 
        1       170 

The first column excluded above is our target variable class. We should not forget to do the same transformation to the test set.

test[-1] <- test[-1] %>% 
  modify(~replace(., .=="na", NA)) %>%
  modify(., as.double)

If we try to apply the summary function on the entire variables (170), we will spent a lot of time to read the summary of each variable without much gain. Instead, we try to get an automated way to get only the information needed to build efficiently our model. To decide whether we should normalize the data or not, for instance, we display the standard deviances of all the variable in decreasing order.

Note: with tree based models we do not need neither normalize the data nor converting factors to dummies.

map_dbl(train[-1], sd, na.rm=TRUE) %>% 
  tibble(sd = .) %>% 
  arrange(-sd)

# A tibble: 170 x 1
           sd
        <dbl>
 1 794874918.
 2  97484780.
 3  42746746.
 4  40404413.
 5  40404412.
 6  40404411.
 7  11567771.
 8  10886737.
 9  10859905.
10  10859904.
# ... with 160 more rows

We have very large variability, which means that the data should be normalized for any machine learning model that uses gradient descent or based on class distances.

Another thing we can check is if some variabels have small number of unique values which can hence be converted to factor type.

map(train[-1], unique) %>% 
  lengths(.) %>% 
  sort(.) %>% 
  head(5)

cd_000 ch_000 as_000 ef_000 ab_000 
     2      3     22     29     30 

To make things simple we consider only the first two ones to be converted to factor type.

the first one is constant which is of type zero variance because its variance equals to zero, and the second one should be converted to factor type with two levels (for the two sets), but since it has large missing values we will decide about it later on. Notice that we do not apply theses transformations here because they will be combined at once with all the required transformations as what will be shown shortly.

Missing values

The best way to deal with missing values depends on their number compared to the dataset size. if we have small number then it would be easier to completely remove them from the data, if in contrast we have large number then the best choice is to impute them using one of the common methods designed for this type of issue.

dim(train[!complete.cases(train),])

[1] 59409   171

As we see almost every row contains at least one missing value in some columns. Let’s check the distribution of missing values within columns.

df <- modify(train[-1], is.na) %>% 
  colSums() %>%
  tibble(names = colnames(train[-1]),missing_values=.) %>% 
  arrange(-missing_values)
  
df

# A tibble: 170 x 2
   names  missing_values
   <chr>           <dbl>
 1 br_000          49264
 2 bq_000          48722
 3 bp_000          47740
 4 bo_000          46333
 5 ab_000          46329
 6 cr_000          46329
 7 bn_000          44009
 8 bm_000          39549
 9 bl_000          27277
10 bk_000          23034
# ... with 160 more rows

I think the best strategy is to first remove columns that have a large number of missing values then we impute the rest, thereby we reduce the number of predictors and the number of missing values at ones. The following script keep the predictors that have a number of missing values less than 10000.

names <- modify(train[-1], is.na) %>% 
  colSums() %>%
  tibble(names = colnames(train[-1]), missing_values=.) %>% 
  filter(missing_values < 10000) 
train1 <- train[c("class",names$names)]
test1 <- test[c("class",names$names)]

An important thing should be noted here is that, if we use imputation methods that use information from all other columns and/or rows to predict the current missing value, therefore the data must be first split between training and testing sets before any imputation, to abide by the crucial rule of machine learning: the test data should never be seen by the model during training process. Fortunately, our data is already split so that the imputation can be done separately. However, the imputation methods will be implemented later on by the help of the recipes package where we bundle all the pre-processing steps together. Note: the above ch_000 was removed since it did not fulfill the required threshold.

imbalanced data

Another important issue that we face when predicting this data is the imbalanced problem.

prop.table(table(train1$class))

       neg        pos 
0.98333333 0.01666667 

This data is highly imbalanced, which tends to make even the worst machine learning model gives very high accuracy rate. In other words, if we do not use any model and predict every class as the largest class label (in our case negative) the accuracy rate will be approximately equal to the proportion of the largest class (in our case 98%), which is very big misleading result. Moreover, this misleading result can be catastrophic if we are more interested to predict the small class (in our case positive) such as detecting fraudulent credit cards. If you would like to get more detail about how to deal with imbalanced data please check this article.

model training

Everything now is ready to run our model with the default values.

model_rf <- rf_wf %>% 
  fit(data = train1)

We can extract the summary of this model as follows

model_rf %>% pull_workflow_fit()

parsnip model object

Fit time:  55.7s 
Ranger result

Call:
 ranger::ranger(formula = ..y ~ ., data = data, num.trees = ~100,      num.threads = ~3, seed = ~123, verbose = FALSE, probability = TRUE) 

Type:                             Probability estimation 
Number of trees:                  100 
Sample size:                      118000 
Number of independent variables:  95 
Mtry:                             9 
Target node size:                 10 
Variable importance mode:         none 
Splitrule:                        gini 
OOB prediction error (Brier s.):  0.003998112 

This model has created 100 trees and has chosen randomly 9 predictors with each tree. with these settings thus we do obtain very low oob error rate which is 0.4% (accuracy rate 99.6% ). However, be cautious with such high accuracy rate, since, in practice, This result may highly related to an overfitting problem. Last thing I want to mention about this output, by looking at the confusion matrix, is the fact that we have now balanced data.

model evaluation

The best way to evaluate our model is by using the testing set. Notice that the yardstick provides bunch of metrics to use, but let’s use the most popular one for classification problems accuracy.

model_rf %>% 
  predict( new_data = test1) %>% 
  bind_cols(test1["class"]) %>% 
  accuracy(truth= as.factor(class), .pred_class) 

# A tibble: 1 x 3
  .metric  .estimator .estimate
  <chr>    <chr>          <dbl>
1 accuracy binary         0.990

with this model we get high accuracy which is very closer to the previous one. However, we should not forget that we are dealing with imbalanced data, and even though we have used subsampling methods (like smote method used here), they do not completely solve this issue, they can only minimize it at certain level and this is the reason why we have many of these methods. Therefore, it is better to use the confusion matrix from the caret package since it gives more information.

caret::confusionMatrix(as.factor(test1$class), predict(model_rf, new_data = test1)$.pred_class)

Confusion Matrix and Statistics

          Reference
Prediction   neg   pos
       neg 15532    93
       pos    64   311
                                          
               Accuracy : 0.9902          
                 95% CI : (0.9885, 0.9917)
    No Information Rate : 0.9748          
    P-Value [Acc > NIR] : < 2e-16         
                                          
                  Kappa : 0.7934          
                                          
 Mcnemar's Test P-Value : 0.02544         
                                          
            Sensitivity : 0.9959          
            Specificity : 0.7698          
         Pos Pred Value : 0.9940          
         Neg Pred Value : 0.8293          
             Prevalence : 0.9748          
         Detection Rate : 0.9708          
   Detection Prevalence : 0.9766          
      Balanced Accuracy : 0.8828          
                                          
       'Positive' Class : neg             
                                          

As said shortly, the specificity rate related to the minor class 78% is very low compared to the major class 99%, and You can think of this as a partial overfitting towards the major class. So if we are more interested to the minor class (which is often the case) then we have go back to our model and try tuning our model or try another subsampling method.

Model tuning:

For model tuning we try other values for some arguments rather than the default vaues. and leave the tuning for some others to the dials package. So let’s try the following argument values:

• num.trees = 100. The default is 500.
• num.threads = 3. The default is 1.

And tune the following:

• mtry = tune(). The default is square root of the number of the variables.
• min_n = tune(). The default is 1.

First, we define the model with these new arguments.

model_tune <- rand_forest(trees= 100, mtry=tune(), min_n = tune()) %>%
  set_engine("ranger", num.threads=3, seed=123) %>% 
  set_mode("classification")

Since in grid search the two arguments mtry and min_n are data dependent, then we should at least specify their ranges.

grid <- grid_regular(mtry(range = c(9,15)), min_n(range = c(5,40)), levels = 3)
grid

# A tibble: 9 x 2
   mtry min_n
  <int> <int>
1     9     5
2    12     5
3    15     5
4     9    22
5    12    22
6    15    22
7     9    40
8    12    40
9    15    40

By setting the levels equal to 3 we get 9 combinations and hence 9 models will be trained. The above recipe has steps that should not be repeated many times when tuning the model, we apply therefore the recipe to the training data in order to get the transformed data, and do not forget to apply the recipe to the testing data.

train2 <- prep(data_rec) %>% 
  juice()
test2 <- prep(data_rec) %>% 
  bake(test1)

To tune our model we use cross validation technique. since we have large data set we use only 3 folds.

set.seed(111)
fold <- vfold_cv(train2, v = 3, strata = class)

Now we bundle our recipe with the specified model.

tune_wf <- workflow() %>% 
  add_model(model_tune) %>%
  add_formula(class~.)

To fit these models across the folds we use the tune_grid function instead of fit.

tune_rf <- tune_wf %>% 
  tune_grid(resamples = fold, grid = grid)

For classification problems this function uses two metrics: accuracy and area under the ROC curve. SO we can extract the metric values as follows.

results <- tune_rf %>% collect_metrics()

To get the best model we have to choose one of the two metrics, so let’s go ahead with the accuracy rate.

best_param <- 
  tune_rf %>% select_best(metric = "accuracy")
best_param

# A tibble: 1 x 3
   mtry min_n .config
  <int> <int> <chr>  
1    15     5 Model3 

we can finalize the workflow with the new parameter values.

tune_wf2 <- tune_wf %>% 
  finalize_workflow(best_param)
tune_wf2

== Workflow =======================
Preprocessor: Formula
Model: rand_forest()

-- Preprocessor -------------------
class ~ .

-- Model --------------------------
Random Forest Model Specification (classification)

Main Arguments:
  mtry = 15
  trees = 100
  min_n = 5

Engine-Specific Arguments:
  num.threads = 3
  seed = 123

Computational engine: ranger 

Now we fit the model with the best parameter values to the entire training data.

best_model <- tune_wf2 %>% 
  fit(train2)
best_model

== Workflow [trained] =============
Preprocessor: Formula
Model: rand_forest()

-- Preprocessor -------------------
class ~ .

-- Model --------------------------
Ranger result

Call:
 ranger::ranger(formula = ..y ~ ., data = data, mtry = ~15L, num.trees = ~100,      min.node.size = ~5L, num.threads = ~3, seed = ~123, verbose = FALSE,      probability = TRUE) 

Type:                             Probability estimation 
Number of trees:                  100 
Sample size:                      118000 
Number of independent variables:  95 
Mtry:                             15 
Target node size:                 5 
Variable importance mode:         none 
Splitrule:                        gini 
OOB prediction error (Brier s.):  0.00359659 

Let’s get the confusion matrix

caret::confusionMatrix(as.factor(test2$class), predict(best_model, new_data = test2)$.pred_class)

Confusion Matrix and Statistics

          Reference
Prediction   neg   pos
       neg 15538    87
       pos    67   308
                                          
               Accuracy : 0.9904          
                 95% CI : (0.9887, 0.9918)
    No Information Rate : 0.9753          
    P-Value [Acc > NIR] : <2e-16          
                                          
                  Kappa : 0.7951          
                                          
 Mcnemar's Test P-Value : 0.1258          
                                          
            Sensitivity : 0.9957          
            Specificity : 0.7797          
         Pos Pred Value : 0.9944          
         Neg Pred Value : 0.8213          
             Prevalence : 0.9753          
         Detection Rate : 0.9711          
   Detection Prevalence : 0.9766          
      Balanced Accuracy : 0.8877          
                                          
       'Positive' Class : neg             
                                          

As we see we do not get any improvement for the specificity rate. so let’s try another subsampling method, say Rose method.

rf_rose <- rand_forest(trees = 100, mtry=9, min_n = 5) %>% 
  set_engine("ranger", num.threads=3, seed = 123) %>%
  set_mode("classification")
data_rec2 <- recipe(class~., data=train1) %>% 
  step_medianimpute(all_predictors() , seed_val = 111) %>% 
  step_zv(all_predictors()) %>% 
  step_corr(all_predictors(), threshold = 0.8) %>% 
  step_normalize(all_predictors()) %>%
  step_rose(class) 
rf_rose_wf <- workflow() %>% 
  add_model(rf_rose) %>% 
  add_recipe(data_rec2)
model_rose_rf <- rf_rose_wf %>% 
  fit(data = train1)
caret::confusionMatrix(as.factor(test1$class), predict(model_rose_rf, new_data = test1)$.pred_class)

Confusion Matrix and Statistics

          Reference
Prediction   neg   pos
       neg 15522   103
       pos   140   235
                                          
               Accuracy : 0.9848          
                 95% CI : (0.9828, 0.9867)
    No Information Rate : 0.9789          
    P-Value [Acc > NIR] : 2.437e-08       
                                          
                  Kappa : 0.6514          
                                          
 Mcnemar's Test P-Value : 0.02092         
                                          
            Sensitivity : 0.9911          
            Specificity : 0.6953          
         Pos Pred Value : 0.9934          
         Neg Pred Value : 0.6267          
             Prevalence : 0.9789          
         Detection Rate : 0.9701          
   Detection Prevalence : 0.9766          
      Balanced Accuracy : 0.8432          
                                          
       'Positive' Class : neg             
                                          

The rose method is much worse than smote method since the specificity rate has doped down to 69%.