Where does probably fit in?

Davis Vaughan

2018-12-09

Introduction

An obvious question regarding probably might be: where does this fit in with the rest of the tidymodels ecosystem? Like the other pieces of the ecosystem, probably is designed to be modular, but plays well with other tidymodels packages.

Regarding placement in the modeling workflow, probably best fits in as a post processing step after the model has been fit, but before the model performance has been calculated.

Example

As an example, we’ll use parsnip to fit a logistic regression on some Lending Club loan data, and then use probably to investigate what happens to performance when you vary the threshold of what a “good” loan is.

library(parsnip)
library(probably)
library(dplyr)
data("lending_club")

# I think it makes more sense to have "good" as the first level
# By default it comes as the second level
lending_club <- lending_club %>%
  mutate(Class = relevel(Class, "good"))

# There are a number of columns in this data set, but we will only use a few
# for this example
lending_club <- select(lending_club, Class, annual_inc, verification_status, sub_grade)

lending_club
#> # A tibble: 9,857 x 4
#>    Class annual_inc verification_status sub_grade
#>    <fct>      <dbl> <fct>               <fct>    
#>  1 good       35000 Not_Verified        C4       
#>  2 good       72000 Verified            C1       
#>  3 good       72000 Source_Verified     D1       
#>  4 good      101000 Verified            C3       
#>  5 good       50100 Source_Verified     A4       
#>  6 good       32000 Source_Verified     B5       
#>  7 good       65000 Not_Verified        A1       
#>  8 good      188000 Not_Verified        B2       
#>  9 good       89000 Source_Verified     B3       
#> 10 good       48000 Not_Verified        C2       
#> # … with 9,847 more rows

Let’s split this into 70% training and 30% testing for something to predict on.

# 70% train, 30% test
set.seed(123)
n <- nrow(lending_club)
prop_train <- .7
train_idx <- sample(seq_len(n), floor(.7 * n))

lending_train <- lending_club[train_idx,]
lending_test  <- lending_club[-train_idx,]

Before we do anything, let’s look at the counts of what we are going to be predicting, the Class of the loan.

count(lending_train, Class)
#> # A tibble: 2 x 2
#>   Class     n
#>   <fct> <int>
#> 1 good   6553
#> 2 bad     346

Clearly there is a large imbalance here with the number of good and bad loans. This is probably a good thing for the bank, but poses an interesting issue for us because we might want to ensure we are sensitive to the bad loans and are not overwhelmed by the number of good ones. One thing that we might do is downsample the number of good loans so that the total number of them is more in line with the number of bad loans. We could do this before fitting the model using recipes::step_downsample(), but for now, let’s continue with the data unchanged.

We’ll use parsnip’s logistic_reg() to create a model specification for logistic regression, set the engine to be glm and then actually fit the model using our data and the model formula.

logi_reg <- logistic_reg()
logi_reg_glm <- logi_reg %>% set_engine("glm")

# A small model specification that defines the type of model you are 
# using and the engine
logi_reg_glm
#> Logistic Regression Model Specification (classification)
#> 
#> Computational engine: glm

# Fit the model
logi_reg_fit <- fit(
  logi_reg_glm, 
  formula = Class ~ annual_inc + verification_status + sub_grade, 
  data = lending_train
)

logi_reg_fit
#> parsnip model object
#> 
#> 
#> Call:  stats::glm(formula = formula, family = stats::binomial, data = data)
#> 
#> Coefficients:
#>                        (Intercept)                          annual_inc  
#>                         -5.685e+00                           2.233e-06  
#> verification_statusSource_Verified         verification_statusVerified  
#>                          8.465e-02                           3.911e-01  
#>                        sub_gradeA2                         sub_gradeA3  
#>                         -1.216e+01                           9.129e-02  
#>                        sub_gradeA4                         sub_gradeA5  
#>                         -3.779e-02                           1.529e+00  
#>                        sub_gradeB1                         sub_gradeB2  
#>                          1.545e+00                          -5.800e-01  
#>                        sub_gradeB3                         sub_gradeB4  
#>                          1.868e+00                           1.571e+00  
#>                        sub_gradeB5                         sub_gradeC1  
#>                          2.077e+00                           2.188e+00  
#>                        sub_gradeC2                         sub_gradeC3  
#>                          2.239e+00                           2.436e+00  
#>                        sub_gradeC4                         sub_gradeC5  
#>                          2.626e+00                           3.062e+00  
#>                        sub_gradeD1                         sub_gradeD2  
#>                          3.039e+00                           2.838e+00  
#>                        sub_gradeD3                         sub_gradeD4  
#>                          3.024e+00                           3.144e+00  
#>                        sub_gradeD5                         sub_gradeE1  
#>                          3.387e+00                           3.326e+00  
#>                        sub_gradeE2                         sub_gradeE3  
#>                          3.342e+00                           3.266e+00  
#>                        sub_gradeE4                         sub_gradeE5  
#>                          3.283e+00                           3.357e+00  
#>                        sub_gradeF1                         sub_gradeF2  
#>                          3.495e+00                           3.833e+00  
#>                        sub_gradeF3                         sub_gradeF4  
#>                          3.967e+00                           3.661e+00  
#>                        sub_gradeF5                         sub_gradeG1  
#>                          3.778e+00                           2.929e+00  
#>                        sub_gradeG2                         sub_gradeG3  
#>                          4.611e+00                          -1.231e+01  
#>                        sub_gradeG4                         sub_gradeG5  
#>                          5.020e+00                           3.489e+00  
#> 
#> Degrees of Freedom: 6898 Total (i.e. Null);  6861 Residual
#> Null Deviance:       2745 
#> Residual Deviance: 2448  AIC: 2524

The output of the parsnip fit() call is a parsnip model_fit object, but the underlying print method for the glm fit is used.

Now let’s predict on our testing set, and use type = "prob" to get class probabilities back rather than hard predictions. We will use these with probably to investigate performance.

predictions <- logi_reg_fit %>%
  predict(new_data = lending_test, type = "prob")

head(predictions, n = 2)
#> # A tibble: 2 x 2
#>   .pred_good .pred_bad
#>        <dbl>     <dbl>
#> 1      0.950    0.0500
#> 2      0.933    0.0671

lending_test_pred <- bind_cols(predictions, lending_test)

lending_test_pred
#> # A tibble: 2,958 x 6
#>    .pred_good .pred_bad Class annual_inc verification_status sub_grade
#>         <dbl>     <dbl> <fct>      <dbl> <fct>               <fct>    
#>  1      0.950   0.0500  good       72000 Verified            C1       
#>  2      0.933   0.0671  good      101000 Verified            C3       
#>  3      0.914   0.0863  bad        80000 Source_Verified     C5       
#>  4      0.958   0.0421  good       56000 Not_Verified        C3       
#>  5      0.892   0.108   good       65000 Verified            D1       
#>  6      0.998   0.00207 good       39676 Not_Verified        B2       
#>  7      0.972   0.0278  good       75000 Verified            B4       
#>  8      0.961   0.0393  good       73400 Source_Verified     C2       
#>  9      0.898   0.102   good       36000 Verified            D1       
#> 10      0.900   0.100   good       36000 Source_Verified     E1       
#> # … with 2,948 more rows

With our class probabilities in hand, we can use make_two_class_pred() to convert these probabilities into hard predictions using a threshold. A threshold of 0.5 just says that if the predicted probability is above 0.5, then classify this prediction as a “good” loan, otherwise, bad.

hard_pred_0.5 <- lending_test_pred %>%
  mutate(.pred = make_two_class_pred(.pred_good, levels(Class), threshold = .5)) %>%
  select(Class, contains(".pred"))

hard_pred_0.5 %>% 
  count(.truth = Class, .pred)
#> # A tibble: 2 x 3
#>   .truth      .pred     n
#>   <fct>  <clss_prd> <int>
#> 1 good         good  2787
#> 2 bad          good   171

Hmm, with a 0.5 threshold, all loans were predicted as “good”. Perhaps this has something to do with the large class imbalance. On the other hand, the bank might want to be more stringent with what is classified as a “good” loan, and might require a probability of 0.75 as the threshold.

hard_pred_0.75 <- lending_test_pred %>%
  mutate(.pred = make_two_class_pred(.pred_good, levels(Class), threshold = .75)) %>%
  select(Class, contains(".pred"))

hard_pred_0.75 %>% 
  count(.truth = Class, .pred)
#> # A tibble: 4 x 3
#>   .truth      .pred     n
#>   <fct>  <clss_prd> <int>
#> 1 good         good  2777
#> 2 good          bad    10
#> 3 bad          good   166
#> 4 bad           bad     5

In this case, 5 of the bad loans were correctly classified as bad, but some of the good loans were also misclassified as bad now. There is a tradeoff here, which can be somewhat captured by the metrics sensitivity and specificity. Both metrics have a max value of 1.

library(yardstick)

# Currently yardstick can't deal with the class_pred objects that come from
# probably, but it will be able to soon!
hard_pred_0.5 <- mutate(hard_pred_0.5, .pred = as.factor(.pred))
hard_pred_0.75 <- mutate(hard_pred_0.75, .pred = as.factor(.pred))

sens(hard_pred_0.5, Class, .pred)
#> # A tibble: 1 x 3
#>   .metric .estimator .estimate
#>   <chr>   <chr>          <dbl>
#> 1 sens    binary             1
spec(hard_pred_0.5, Class, .pred)
#> # A tibble: 1 x 3
#>   .metric .estimator .estimate
#>   <chr>   <chr>          <dbl>
#> 1 spec    binary             0

sens(hard_pred_0.75, Class, .pred)
#> # A tibble: 1 x 3
#>   .metric .estimator .estimate
#>   <chr>   <chr>          <dbl>
#> 1 sens    binary         0.996
spec(hard_pred_0.75, Class, .pred)
#> # A tibble: 1 x 3
#>   .metric .estimator .estimate
#>   <chr>   <chr>          <dbl>
#> 1 spec    binary        0.0292

In this example, as we increased specificity (by capturing those 5 bad loans with a higher threshold), we lowered sensitivity (by incorrectly reclassifying some of the good loans as bad). It would be nice to have some combination of these metrics to represent this tradeoff. Luckily, j_index is exactly that.

\[ j\_index = sens + spec - 1 \]

j_index has a maximum value of 1 when there are no false positives and no false negatives. It can be used as justification of whether or not an increase in the threshold value is worth it. If increasing the threshold results in more of an increase in the specificity than a decrease in the sensitivity, we can see that with j_index.

j_index(hard_pred_0.5, Class, .pred)
#> # A tibble: 1 x 3
#>   .metric .estimator .estimate
#>   <chr>   <chr>          <dbl>
#> 1 j_index binary             0
j_index(hard_pred_0.75, Class, .pred)
#> # A tibble: 1 x 3
#>   .metric .estimator .estimate
#>   <chr>   <chr>          <dbl>
#> 1 j_index binary        0.0257

Now, this is not the only way to optimize things. If you care about low false positives, you might be more interested in keeping sensitivity high, and this wouldn’t be the best way to tackle this problem. But for now, let’s see how we can use probably to optimize the j_index.

threshold_perf() will recalculate a number of metrics across varying thresholds. One of these is j_index.

threshold_data <- lending_test_pred %>%
  threshold_perf(Class, .pred_good, thresholds = seq(0.5, 1, by = 0.0025))

threshold_data %>%
  filter(.threshold %in% c(0.5, 0.6, 0.7))
#> # A tibble: 12 x 4
#>    .threshold .metric  .estimator .estimate
#>         <dbl> <chr>    <chr>          <dbl>
#>  1        0.5 sens     binary        1     
#>  2        0.6 sens     binary        0.999 
#>  3        0.7 sens     binary        0.998 
#>  4        0.5 spec     binary        0     
#>  5        0.6 spec     binary        0.0175
#>  6        0.7 spec     binary        0.0234
#>  7        0.5 j_index  binary        0     
#>  8        0.6 j_index  binary        0.0161
#>  9        0.7 j_index  binary        0.0212
#> 10        0.5 distance binary        1     
#> 11        0.6 distance binary        0.965 
#> 12        0.7 distance binary        0.954

With ggplot2, we can easily visualize this varying performance to find our optimal threshold for maximizing j_index.

library(ggplot2)

threshold_data <- threshold_data %>%
  filter(.metric != "distance") %>%
  mutate(group = case_when(
    .metric == "sens" | .metric == "spec" ~ "1",
    TRUE ~ "2"
  ))

max_j_index_threshold <- threshold_data %>%
  filter(.metric == "j_index") %>%
  filter(.estimate == max(.estimate)) %>%
  pull(.threshold)

ggplot(threshold_data, aes(x = .threshold, y = .estimate, color = .metric, alpha = group)) +
  geom_line() +
  theme_minimal() +
  scale_color_viridis_d(end = 0.9) +
  scale_alpha_manual(values = c(.4, 1), guide = "none") +
  geom_vline(xintercept = max_j_index_threshold, alpha = .6, color = "grey30") +
  labs(
    x = "'Good' Threshold\n(above this value is considered 'good')",
    y = "Metric Estimate",
    title = "Balancing performance by varying the threshold",
    subtitle = "Sensitivity or specificity alone might not be enough!\nVertical line = Max J-Index"
  )

It’s clear from this visual that the optimal threshold is very high, exactly 0.9475. This is pretty high, so again, this optimization method won’t be useful for all cases. To wrap up, here are all of the test set metrics for that threshold value.

threshold_data %>%
  filter(.threshold == max_j_index_threshold)
#> # A tibble: 3 x 5
#>   .threshold .metric .estimator .estimate group
#>        <dbl> <chr>   <chr>          <dbl> <chr>
#> 1      0.948 sens    binary         0.675 1    
#> 2      0.948 spec    binary         0.690 1    
#> 3      0.948 j_index binary         0.365 2