1.1 Kaggle Data Description

The files contain complete loan data for all loans issued through 2007-2015, including the current loan status (Current, Late, Fully Paid, etc.) and latest payment information. The file containing loan data through the “present” contains complete loan data for all loans issued through the previous completed calendar quarter. Additional features include credit scores, number of finance inquiries, address including zip codes, and state, and collections among others. The file is a matrix of about 890 thousand observations and 75 variables. A data dictionary is provided in a separate file.

2.1 Libraries

This will automatically install any missing libraries. Make sure to have proper proxy settings. You do not necessarily need to have all of them but only those needed for respective sections (which will be indicated).

# define used libraries
libraries_used <- 
  c("lazyeval", "readr","plyr" ,"dplyr", "readxl", "ggplot2", 
    "funModeling", "scales", "tidyverse", "corrplot", "GGally", "caret",
    "rpart", "randomForest", "pROC", "gbm", "choroplethr", "choroplethrMaps",
    "microbenchmark", "doParallel", "e1071")

# check missing libraries
libraries_missing <- 
  libraries_used[!(libraries_used %in% installed.packages()[,"Package"])]
# install missing libraries
if(length(libraries_missing)) install.packages(libraries_missing)

2.2 Session Info

sessionInfo()

## R version 3.4.1 (2017-06-30)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 15063)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] methods   stats     graphics  grDevices utils     datasets  base     
## 
## loaded via a namespace (and not attached):
##  [1] compiler_3.4.1  backports_1.1.0 bookdown_0.5    magrittr_1.5   
##  [5] rprojroot_1.2   tools_3.4.1     htmltools_0.3.6 yaml_2.1.14    
##  [9] Rcpp_0.12.12    stringi_1.1.5   rmarkdown_1.6   blogdown_0.1   
## [13] knitr_1.17      stringr_1.2.0   digest_0.6.12   evaluate_0.10.1

3.1 Library versions

Also note that often the version of the library used matters (it always matters but what we mean is, it makes a difference) and some libraries are developed more actively than others which may lead to issues. For example, note that library caret is using plyr while the tidyverse of which dplyr (successor of plyr) is part has some new concepts, e.g. the default data frame is a tibble::tibble(). It seems that caret has issues with tibbles, see e.g. “Wrong model type for classification” in regression problems in R-Caret.

library(plyr)
library(tidyverse)

## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr

## Conflicts with tidy packages ----------------------------------------------

## arrange():   dplyr, plyr
## compact():   purrr, plyr
## count():     dplyr, plyr
## failwith():  dplyr, plyr
## filter():    dplyr, stats
## id():        dplyr, plyr
## lag():       dplyr, stats
## mutate():    dplyr, plyr
## rename():    dplyr, plyr
## summarise(): dplyr, plyr
## summarize(): dplyr, plyr

library(caret)

## Loading required package: lattice

## 
## Attaching package: 'caret'

## The following object is masked from 'package:purrr':
## 
##     lift

5.1 Data Transformation

We should also check the variable types. As mentioned in section Data Import the import function uses heuristics to guess the data types and these may not always work.

• annual_inc_joint = character but should be numeric
• dates are read as character and need to be transformed

First let’s transform character to numeric. We make use of dplyr::mutate_at() function and provide a vector of columns to be mutated (transformed). In general, when using libraries from the tidyverse (these libraries are mainly authored by Hadley Wickham and other RStudio people), most functions offer a standard version as opposed to an NSE (non-standard evaluation) version which can take character values as variable names. These functions usually have a trailing underscore, e.g. dplyr::select_() as compared to the non-standard evaluation function dplyr::select(). For details, see the dplyr vignette on Non-standard evaluation. However, note that the tidyverse libraries are still changing a fair amount so best to cross-check whether used function is the most efficient one for given task or if there is a new one introduced in the mean time.

As opposed to packages like data.table tidyverse packages do not alter their input in place so we have to re-assign the result to the original object. We could also create a new object but this is memory-inefficient (in fact even the re-assignment to the original object creates a temporary copy). For more details on memory management, see chapter Memory in Hadley Wickham’s book “Advanced R”.

chr_to_num_vars <- 
  c("annual_inc_joint", "mths_since_last_major_derog", "open_acc_6m",
    "open_il_6m", "open_il_12m", "open_il_24m", "mths_since_rcnt_il",
    "total_bal_il", "il_util", "open_rv_12m", "open_rv_24m",
    "max_bal_bc", "all_util", "total_rev_hi_lim", "total_cu_tl",
    "inq_last_12m", "dti_joint", "inq_fi", "tot_cur_bal", "tot_coll_amt")

loans <-
  loans %>%
  mutate_at(.funs = funs(as.numeric), .vars = chr_to_num_vars)

Let’s have a look at the date variables to see how they need to be transformed.

chr_to_date_vars <- 
  c("issue_d", "last_pymnt_d", "last_credit_pull_d",
    "next_pymnt_d", "earliest_cr_line", "next_pymnt_d")

loans %>%
  select_(.dots = chr_to_date_vars) %>%
  str()

## Classes 'tbl_df', 'tbl' and 'data.frame':    887379 obs. of  5 variables:
##  $ issue_d           : chr  "Dec-2011" "Dec-2011" "Dec-2011" "Dec-2011" ...
##  $ last_pymnt_d      : chr  "Jan-2015" "Apr-2013" "Jun-2014" "Jan-2015" ...
##  $ last_credit_pull_d: chr  "Jan-2016" "Sep-2013" "Jan-2016" "Jan-2015" ...
##  $ next_pymnt_d      : chr  NA NA NA NA ...
##  $ earliest_cr_line  : chr  "Jan-1985" "Apr-1999" "Nov-2001" "Feb-1996" ...

head(unique(loans$next_pymnt_d))

## [1] NA         "Feb-2016" "Jan-2016" "Sep-2013" "Feb-2014" "May-2014"

for (i in chr_to_date_vars){
  print(head(unique(loans[, i])))
  }

## # A tibble: 6 x 1
##    issue_d
##      <chr>
## 1 Dec-2011
## 2 Nov-2011
## 3 Oct-2011
## 4 Sep-2011
## 5 Aug-2011
## 6 Jul-2011
## # A tibble: 6 x 1
##   last_pymnt_d
##          <chr>
## 1     Jan-2015
## 2     Apr-2013
## 3     Jun-2014
## 4     Jan-2016
## 5     Apr-2012
## 6     Nov-2012
## # A tibble: 6 x 1
##   last_credit_pull_d
##                <chr>
## 1           Jan-2016
## 2           Sep-2013
## 3           Jan-2015
## 4           Sep-2015
## 5           Dec-2014
## 6           Aug-2012
## # A tibble: 6 x 1
##   next_pymnt_d
##          <chr>
## 1         <NA>
## 2     Feb-2016
## 3     Jan-2016
## 4     Sep-2013
## 5     Feb-2014
## 6     May-2014
## # A tibble: 6 x 1
##   earliest_cr_line
##              <chr>
## 1         Jan-1985
## 2         Apr-1999
## 3         Nov-2001
## 4         Feb-1996
## 5         Jan-1996
## 6         Nov-2004
## # A tibble: 6 x 1
##   next_pymnt_d
##          <chr>
## 1         <NA>
## 2     Feb-2016
## 3     Jan-2016
## 4     Sep-2013
## 5     Feb-2014
## 6     May-2014

It seems the date format is consistent and follows month-year convention. We can use the base::as.Date() function to convert this to a date format. The function requires a day as well so we simply add the 1st of each month to the existing character string via pasting together strings using base::paste0(). Note that base::paste0("bla", NA) will not return NA but the concatenated string (here: blaNA). Conveniently, base::as.Date() will return NA so we can leave it at that. Alternatively, we could include an exception handler for NA values to explicitly handle those because we saw previously that some date variables include NA values. One way to approach that would be to wrap the date conversion function into a base::ifelse() call as so ifelse(is.na(x), NA, some_date_function(). However, it seems that base::ifelse() is dropping the date class that we just created with our date function, for details, see e.g. How to prevent ifelse from turning Date objects into numeric objects. As we do not want to deal with these issues at this moment, we simply go with the default behavior of base::as.Date() as it returns NA in cases of bad input anyway.

Let’s also return the NA values of the input to remind ourselves of the number. It should coincide with the number of NA after we have converted / transformed the date variables.

meta_loans %>% 
  select(variable, q_na) %>% 
  filter(variable %in% chr_to_date_vars)

##             variable   q_na
## 1            issue_d      0
## 2   earliest_cr_line     29
## 3       last_pymnt_d  17659
## 4       next_pymnt_d 252971
## 5 last_credit_pull_d     53

Finally, this is how our custom date conversion function will look like, we call it convert_date() and it will take the string value to be converted as input. We define the date format by following the function conventions of base::as.Date(), for details see ?base::as.Date. We could have also used an anonymous function directly in the dplyr::mutate_at() call but the creation of a specific function seemed appropriate as we may use it several times.

convert_date <- function(x){
  as.Date(paste0("01-", x), format = "%d-%b-%Y")
  } 

loans <-
  loans %>%
  mutate_at(.funs = funs(convert_date), .vars = chr_to_date_vars)

num_vars <- 
  loans %>% 
  sapply(is.numeric) %>% 
  which() %>% 
  names()

meta_loans %>%
  select(variable, p_zeros, p_na, unique) %>%
  filter_(~ variable %in% num_vars) %>%
  knitr::kable()

We also see that variables mths_since_last_delinq, mths_since_last_record, mths_since_last_major_derog, dti_joint and annual_inc_joint have a large share of NA values. If we think about this in more detail, it may be reasonable to assume that NA values for the variables mths_since_last_delinq, mths_since_last_record and mths_since_last_major_derog actually indicate that there was no event/record of any missed payment so there cannot be any time value. Analogously, a missing value for annual_inc_joint and dti_joint may simply indicate that it is a single borrower or the partner has no income. Thus, the first three variables actually carry valuable information that may be lost if we ignored it. We will thus replace the missing values with zeros to make them available for modeling. It should be noted though that a zero time could indicate an event that is just happening so we have to document our assumptions carefully.

na_to_zero_vars <-
  c("mths_since_last_delinq", "mths_since_last_record",
    "mths_since_last_major_derog")

loans <- 
  loans %>%
  mutate_at(.vars = na_to_zero_vars, .funs = funs(replace(., is.na(.), 0)))

These transformations should have us covered for now. We recreate the meta table after all these changes.

meta_loans <- funModeling::df_status(loans, print_results = FALSE)
meta_loans <-
  meta_loans %>%
  mutate(uniq_rat = unique / nrow(loans))

5.2 Additional Meta Data File

Next we look at the meta descriptions provided in the additional Excel file.

knitr::kable(meta_loan_stats[,1:2])

As expected, id is “A unique LC assigned ID for the loan listing.” and member_id is “A unique LC assigned Id for the borrower member.”. The policy_code is “publicly available policy_code=1 new products not publicly available policy_code=2”. That means there could be different values but as we have seen before in this data it only takes on one value.

Finally, let’s look at variables that are either in the data and not in the meta description or vice versa. The dplyr::setdiff() function does what its name suggests, just pay attention to the order of arguments to understand which difference you actually get.

Variables in loans data but not in meta description.

dplyr::setdiff(colnames(loans), meta_loan_stats$LoanStatNew)

## [1] "verification_status_joint" "total_rev_hi_lim"

Variables in meta description but not in loans data.

dplyr::setdiff(meta_loan_stats$LoanStatNew, colnames(loans))

##  [1] "acc_open_past_24mths"           "avg_cur_bal"                   
##  [3] "bc_open_to_buy"                 "bc_util"                       
##  [5] "chargeoff_within_12_mths"       "delinq_amnt"                   
##  [7] "fico_range_high"                "fico_range_low"                
##  [9] "last_fico_range_high"           "last_fico_range_low"           
## [11] "mo_sin_old_il_acct"             "mo_sin_old_rev_tl_op"          
## [13] "mo_sin_rcnt_rev_tl_op"          "mo_sin_rcnt_tl"                
## [15] "mort_acc"                       "mths_since_recent_bc"          
## [17] "mths_since_recent_bc_dlq"       "mths_since_recent_inq"         
## [19] "mths_since_recent_revol_delinq" "num_accts_ever_120_pd"         
## [21] "num_actv_bc_tl"                 "num_actv_rev_tl"               
## [23] "num_bc_sats"                    "num_bc_tl"                     
## [25] "num_il_tl"                      "num_op_rev_tl"                 
## [27] "num_rev_accts"                  "num_rev_tl_bal_gt_0"           
## [29] "num_sats"                       "num_tl_120dpd_2m"              
## [31] "num_tl_30dpd"                   "num_tl_90g_dpd_24m"            
## [33] "num_tl_op_past_12m"             "pct_tl_nvr_dlq"                
## [35] "percent_bc_gt_75"               "pub_rec_bankruptcies"          
## [37] "tax_liens"                      "tot_hi_cred_lim"               
## [39] "total_bal_ex_mort"              "total_bc_limit"                
## [41] "total_il_high_credit_limit"     "total_rev_hi_lim  "            
## [43] "verified_status_joint"          NA

It seems for the loans variables verification_status_joint and total_rev_hi_lim there are actually equivalents in the meta data but they carry slightly different names or have leading / trailing blanks. There are quite a few variables in the meta description that are not in the data. But that should be less of a concern as we don’t need those anyway.

9.1 Grade

For a start, let’s look at grade which seems to be a rating classification scheme that Lending Club uses to assign loans into risk buckets similar to other popular rating schemes like S&P or Moodys. For more details, see Lending Club Rate Information. For now, it suffices to know that grades take values A, B, C, D, E, F, G where A represents the highest quality loan and G the lowest.

As a first relation, we investigate the distribution of loan amount over the different grades with a standard boxplot highlighting potential outliers in red. One may also select a pre-defined theme from the ggthemes package by adding a call to ggplot such as theme_economist() or even create a theme oneself. For details on themes, see ggplot2 themes and Introduction to ggthemes.

As mentioned before we want to be group size aware so let’s write a short function to be used inside the ggplot::geom_boxplot() call in combination with ggplot::stat_summary() where we can call user-defined functions using parameter fun.data. As it is often the case with plots, we need to experiment with the position of additional text elements which we achieve by scaling the y-position with some constant multiplier around the median (as the boxplot will have the median as horizontal bar within the rectangles). The mean can be added in a similar fashion, however we don’t need to specify the y-position explicitly. The count as text and the mean as point may overlap themselves and with the median horizontal bar. This is acceptable in an exploratory setting. For publication-ready plots one would have to perform some adjustments.

# see http://stackoverflow.com/questions/15660829/how-to-add-a-number-of-observations-per-group-and-use-group-mean-in-ggplot2-boxp
give_count <- 
  stat_summary(fun.data = function(x) return(c(y = median(x)*1.06,
                                               label = length(x))),
               geom = "text")

# see http://stackoverflow.com/questions/19876505/boxplot-show-the-value-of-mean
give_mean <- 
  stat_summary(fun.y = mean, colour = "darkgreen", geom = "point", 
               shape = 18, size = 3, show.legend = FALSE)

train %>%
  ggplot(aes(grade, loan_amnt)) +
  geom_boxplot(fill = "white", colour = "darkblue", 
               outlier.colour = "red", outlier.shape = 1) +
  give_count +
  give_mean +
  scale_y_continuous(labels = comma) +
  facet_wrap(~ default) +
  labs(title="Loan Amount by Grade", x = "Grade", y = "Loan Amount \n")

We can derive a few points from the plot:

• there is not a lot of difference between default and non-default
• lower quality loans tend to have a higher loan amount
• there are virtually no outliers except for grade B
• the loan amount spread (IQR) seems to be slightly higher for lower quality loans

According to Lending Club’s rate information, we would expect a strong increasing relationship between grade and interest rate so we can test this. We also consider the term of the loan as one might expect that longer term loans could have a higher interest rate.

train %>%
  ggplot(aes(grade, int_rate)) +
  geom_boxplot(fill = "white", colour = "darkblue", 
               outlier.colour = "red", outlier.shape = 1) +
  give_count +
  give_mean +
  scale_y_continuous(labels = comma) +
  labs(title="Interest Rate by Grade", x = "Grade", y = "Interest Rate \n") +
  facet_wrap(~ term)

We can derive a few points from the plot:

• interest rate increases with grade worsening
• a few loans seem to have an equally low interest rate independent of grade
• the spread of rates seems to increase with grade worsening
• there tend to be more outliers on the lower end of the rate
• The 3-year term has a much higher number of high-rated borrowers while the 5-year term has a larger number in the low-rating grades

9.2 Home Ownership

We can do the same for home ownership.

train %>%
  ggplot(aes(home_ownership, int_rate)) +
  geom_boxplot(fill = "white", colour = "darkblue", 
               outlier.colour = "red", outlier.shape = 1) +
  give_count +
  give_mean +
  scale_y_continuous(labels = comma) +
  facet_wrap(~ default) +
  labs(title="Interest Rate by Home Ownership", x = "Home Ownership", y = "Interest Rate \n")

We can derive a few points from the plot:

• there seems no immediate conclusion with respect to the impact of home ownership, e.g. we can see that median/mean interest rate is higher for people who own a house than for those who still pay a mortgage.
• interest rates are highest for values “any” and “none” which could be loans of limited data quality but there are very few data points.
• interest rate is higher for default loans, which is probably driven by other factors (e.g. grade)

9.3 Loan Amount and Income

Another interesting plot may be the relationship between loan amount and funded / invested amount. As all variables are continous, we can do that with a simple scatterplot but we will need to reshape the data to have both loan values plotted against loan amount.

In fact the reshaping here may be slightly odd as we like to display the same variable on the x-axis but different values on the y-axis facetted by their variable names. To achieve that, we gather the data three times with the same x variable but changing y variables.

funded_amnt <-
  train %>%
  transmute(loan_amnt = loan_amnt, value = funded_amnt, 
              variable = "funded_amnt")

funded_amnt_inv <-
  train %>%
  transmute(loan_amnt = loan_amnt, value = funded_amnt_inv, 
              variable = "funded_amnt_inv")

plot_data <- rbind(funded_amnt, funded_amnt_inv)
# remove unnecessary data using regex
ls()

##  [1] "categorical_vars"  "chr_to_date_vars"  "chr_to_num_vars"  
##  [4] "convert_date"      "default_vars"      "defaulted"        
##  [7] "excel_file"        "funded_amnt"       "funded_amnt_inv"  
## [10] "give_count"        "give_mean"         "i"                
## [13] "income_vars"       "libraries_missing" "libraries_used"   
## [16] "loan_amount_vars"  "loans"             "meta_browse_notes"
## [19] "meta_loan_stats"   "meta_loans"        "meta_reject_stats"
## [22] "na_to_zero_vars"   "num_vars"          "path"             
## [25] "plot_data"         "test"              "train"            
## [28] "train_index"       "vars_to_remove"

rm(list = ls()[grep("^funded", ls())])

Now let’s plot the data.

plot_data %>%
  ggplot(aes(x = loan_amnt, y = value)) +
  facet_wrap(~ variable, scales = "free_x", ncol = 3) +
  geom_point()

We can derive a few points from the plot:

• there are instances when funded amount is smaller loan amount
• there seems to be a number of loans where investment is smaller funded amount i.e. not the full loan is invested in

Let’s do the same but only for annual income versus loan amount.

train %>%
  ggplot(aes(x = annual_inc, y = loan_amnt)) +
  geom_point()

## Warning: Removed 4 rows containing missing values (geom_point).

We can derive a few points from the plot:

• there is no immediatly discernible relationship
• there are quite a few income outliers with questionable values (e.g. why would a person with annual income 9500000 request a loan amount of 24000)

9.4 Time Series

Now let’s take a look at interest rates over time but split the time series by grade to see if there are differences in interest rate development depending on the borrower grade. Again we make use of a dplyr pipe, first selecting the variables of interest, then grouping by attributes and finally summarising the metric interest rate by taking the mean for each goup. We use facet_wrap to split by attribute grade. As we are using the mean for building an aggregate representation, we should be weary of outliers and group size which we had already looked at earlier (ignoring time dimension).

train %>%
  select(int_rate, grade) %>%
  group_by(grade) %>%
  summarise(int_rate_mean = mean(int_rate, na.rm = TRUE),
            int_rate_median = median(int_rate, na.rm = TRUE),
            n = n()) %>%
  knitr::kable()

Group size are relatively large except for grade G which may have only a few hundred loans per group when grouping by grade and date. Mean and median seem fairly close at this non-granular level. Let’s plot the data.

train %>%
  select(int_rate, grade, issue_d) %>%
  group_by(grade, issue_d) %>%
  summarise(int_rate_mean = mean(int_rate, na.rm = TRUE)) %>%
  ggplot(aes(issue_d, int_rate_mean)) +
  geom_line(color= "darkblue", size = 1) +
  facet_wrap(~ grade)

We can derive a few points from the plot:

• the mean interest rate is falling or relatively constant for high-rated clients
• the mean interest rate is increasing significantly for low-rated clients

Let’s also look at loan amounts over time in the same manner.

train %>%
  select(loan_amnt, grade, issue_d) %>%
  group_by(grade, issue_d) %>%
  summarise(loan_amnt_mean = mean(loan_amnt, na.rm = TRUE)) %>%
  ggplot(aes(issue_d, loan_amnt_mean)) +
  geom_line(color= "darkblue", size = 1) +
  facet_wrap(~ grade)

We can derive a few points from the plot:

• the mean loan amount is increasing for all grades
• while high-rated clients have some mean loan amount volatility, it is much higher for low-rated clients

9.5 Geolocation Plots

Let’s remind ourselves of the geolocation variables in the data and their information power.

• zip_code (geo-info)
• addr_state (geo-info)

geo_vars <- c("zip_code", "addr_state")

meta_loans %>%
  select(variable, p_zeros, p_na, type, unique) %>%
  filter_(~ variable %in% geo_vars) %>%
  knitr::kable()

loans %>%
  select_(.dots = geo_vars) %>%
  str()

## Classes 'tbl_df', 'tbl' and 'data.frame':    887379 obs. of  2 variables:
##  $ zip_code  : chr  "860xx" "309xx" "606xx" "917xx" ...
##  $ addr_state: chr  "AZ" "GA" "IL" "CA" ...

We see that zip_code seems to be the truncated US postal code with only first three digits having a value. The addr_state seems to be the state names in a two-letter abbreviation.

We can use the choroplethr package to work with maps (alternatives may be maps among other packages). While choroplethr provides functions to work with the data, its sister package choroplethrMaps contains corresponding maps that can be used by choroplethr. One issue with choroplethr is that it attaches the package plyr which means many functions from dplyr are masked as it is the successor to plyr. Thus we do not load the package choroplethr despite using its functions frequently. However, we need to load choroplethrMaps as it has some data that is not exported from its namespace so syntax like choroplethrMaps::state.map where state.map is the rdata file does not work. An alternative might be to load the file directly by going to the underlying directory, e.g. R-3.x.x\library\choroplethrMaps\data but this is cumbersome. As discussed before, when attaching a library, we just need to make sure that important functions of other packages are not masked. In this case it’s fine.

For details on those libraries, see CRAN choroplethr and CRAN choroplethrMaps.

First we aggregate the default rate by state.

default_rate_state <- 
  train %>%
  select(default, addr_state) %>%
  group_by(addr_state) %>%
  summarise(default_rate = sum(default, na.rm = TRUE) / n())

knitr::kable(default_rate_state)

The data is already in a good format but when using choroplethr we need to adhere to some conventions to make it work out of the box. The first thing would be to bring the state names into a standard format. In fact, we can investigate the format by looking e.g. into the choroplethrMaps::state.map dataset which we later use to map our data. We first load the library and then the data via the function utils::data(). From the documentation, we also learn that state.map is a “data.frame which contains a map of all 50 US States plus the District of Columbia.” It is based on a shapefile which is often used in geospatial visualizations and “taken from the US Census 2010 Cartographic Boundary shapefiles page”. We are interested in the region variable which seems to hold the state names.

library(choroplethrMaps)
utils::data(state.map)
str(state.map)

## 'data.frame':    50763 obs. of  12 variables:
##  $ long      : num  -113 -113 -112 -112 -111 ...
##  $ lat       : num  37 37 37 37 37 ...
##  $ order     : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ hole      : logi  FALSE FALSE FALSE FALSE FALSE FALSE ...
##  $ piece     : Factor w/ 66 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ group     : Factor w/ 226 levels "0.1","1.1","10.1",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ id        : chr  "0" "0" "0" "0" ...
##  $ GEO_ID    : Factor w/ 51 levels "0400000US01",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ STATE     : Factor w/ 51 levels "01","02","04",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ region    : chr  "arizona" "arizona" "arizona" "arizona" ...
##  $ LSAD      : Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ...
##  $ CENSUSAREA: num  113594 113594 113594 113594 113594 ...

unique(state.map$region)

##  [1] "arizona"              "arkansas"             "louisiana"           
##  [4] "minnesota"            "mississippi"          "montana"             
##  [7] "new mexico"           "north dakota"         "oklahoma"            
## [10] "pennsylvania"         "tennessee"            "virginia"            
## [13] "california"           "delaware"             "west virginia"       
## [16] "wisconsin"            "wyoming"              "alabama"             
## [19] "alaska"               "florida"              "idaho"               
## [22] "kansas"               "maryland"             "colorado"            
## [25] "new jersey"           "north carolina"       "south carolina"      
## [28] "washington"           "vermont"              "utah"                
## [31] "iowa"                 "kentucky"             "maine"               
## [34] "massachusetts"        "connecticut"          "michigan"            
## [37] "missouri"             "nebraska"             "nevada"              
## [40] "new hampshire"        "new york"             "ohio"                
## [43] "oregon"               "rhode island"         "south dakota"        
## [46] "district of columbia" "texas"                "georgia"             
## [49] "hawaii"               "illinois"             "indiana"

So it seems we need small case long form state names which implies we need a mapping from the short names in our data. Some internet research may give us the site American National Standards Institute (ANSI) Codes for States where we find a mapping under section FIPS Codes for the States and the District of Columbia. We could try using some parsing tool to extract the table from the web page but for now we take a short cut and simply copy paste the data into a tab-delimited txt file and read the file with readr::read_tsv(). We then cross-check if all abbreviations used in our loans data are present in the newly created mapping table.

states <- readr::read_tsv("./data//us_states.txt")

## Parsed with column specification:
## cols(
##   Name = col_character(),
##   `FIPS State Numeric Code` = col_integer(),
##   `Official USPS Code` = col_character()
## )

str(states)

## Classes 'tbl_df', 'tbl' and 'data.frame':    51 obs. of  3 variables:
##  $ Name                   : chr  "Alabama" "Alaska" "Arizona" "Arkansas" ...
##  $ FIPS State Numeric Code: int  1 2 4 5 6 8 9 10 11 12 ...
##  $ Official USPS Code     : chr  "AL" "AK" "AZ" "AR" ...
##  - attr(*, "spec")=List of 2
##   ..$ cols   :List of 3
##   .. ..$ Name                   : list()
##   .. .. ..- attr(*, "class")= chr  "collector_character" "collector"
##   .. ..$ FIPS State Numeric Code: list()
##   .. .. ..- attr(*, "class")= chr  "collector_integer" "collector"
##   .. ..$ Official USPS Code     : list()
##   .. .. ..- attr(*, "class")= chr  "collector_character" "collector"
##   ..$ default: list()
##   .. ..- attr(*, "class")= chr  "collector_guess" "collector"
##   ..- attr(*, "class")= chr "col_spec"

# give some better variable names (spaces never helpful)
names(states) <- c("name", "fips_code", "usps_code")
# see if all states are covered
dplyr::setdiff(default_rate_state$addr_state, states$usps_code)

## character(0)

Comparing the unique abbreviations in the loans data with the usps codes reveals that we are all covered so let’s bring the data together using dplyr::left_join(). For details on dplyr join syntax, see Two-table verbs. To go with the choroplethr conventions, we also need to

• have state names in lower case
• rename the mapping variable to region and the numeric metric to value.
• only select above two variables in the function call

Note that we are re-assigning the data to itself to avoid creating a new table.

default_rate_state <-
  default_rate_state %>%
  left_join(states[, c("usps_code", "name")], 
            by = c("addr_state" = "usps_code")) %>%
  rename(region = name, value = default_rate) %>%
  mutate(region = tolower(region)) %>%
  select(region, value)

Finally, we can plot the data using the default colors from choroplethr. As we are having data on state level, we use function choroplethr::state_choropleth().

choroplethr::state_choropleth(df = default_rate_state, 
                              title = "Default rate by State")

We can derive a few points from the plot:

• default rate varies across states
• we may want to adjust the default binning
• we are only looking at default rate but not defaulted exposure (there may be states with many defaults but the defaulted amount is small)

We could do the same exercise but on a more granular level as we have zip codes available. They are truncated so we would have to find a way to work with them. They also have a high number of unique values so the map legend may be cramped if we were looking at all of US.

12.1 Model options

There are a couple of different model options available for a supervised binary classification problem such as determining whether a loan will default or not. In general, we may distinguish the following model approaches (among others):

• logistic regression
• linear discriminant analysis
• k-nearest neighbors (kNN)
• trees
• random forests
• boosting
• support vector machines (SVM)

In addition, there are multiple ways to assess model performance such as model parameters, cross-validation, bootstrap and multiple feature selection methods such as ridge regression, the lasso, principal component anaylsis, etc.

We will look at some of those models and methods and see how they compare. For a more detailed treatment and further background, an accessible text is “An Introduction to Statistical Learning” or for a more advanced treatment “The Elements of Statistical Learning: Data Mining, Inference, and Prediction”.

12.2 Imbalanced data

Before we can think about any modeling, we have to realize that the dataset is highly imbalanced, i.e. the target variable has a very low proportion of defaults (namely 0.0249963). This can lead to several issues in many algorithms, e.g.

• biased accuracy
• loss functions attempt to optimize quantities such as error rate, not taking the data distribution into consideration
• errors have the same cost (not pertaining to imbalanced data only)

More details can be found in

• Learning from Imbalanced Classes
• Practical Guide to deal with Imbalanced Classification Problems in R
• The caret package: Subsampling For Class Imbalances

There are a couple of approaches to deal with the problem which can be divided into (taken from Learning from Imbalanced Classes)

• do nothing (not a good idea in general)
• balance training set
  • oversample minority class
  • undersample majority class
  • synthesize new minority classes
• throw away minority examples and switch to an anomaly detection framework
• at the algorithm level, or after it:
  • adjust the class weight (misclassification costs)
  • adjust the decision threshold
  • modify an existing algorithm to be more sensitive to rare classes
• construct an entirely new algorithm to perform well on imbalanced data

We will focus on balancing the training set. There are basically three options available

• undersample majority class (randomly or informed)
• oversample minority class (randomly or informed)
• do a mixture of both (e.g. SMOTE or ROSE method)

Undersampling the majority class may lose information but via decreasing dataset also lead to more computational efficiency. We also tried SMOTE and ROSE but functions DMwR::SMOTE() and ROSE::ROSE() seem to be very picky about their input and so far we failed to run them. For details on SMOTE, see SMOTE: Synthetic Minority Over-sampling Technique and for details on ROSE, see ROSE: A Package for Binary Imbalanced Learning. As pointed out here a random sampling for either under- or oversampling is not the best idea. It is recommended to use cross-validation and perform over- or under-sampling on each fold independently to get an honest estimate of model performance. In fact, function caret::train() also allows the re-balancing of data during the training of a model via caret::trainControl(sampling = "..."), see The caret package 11.2: Subsampling During Resampling for details. For now, we go with undersampling which still leaves a fair amount of training observations (at least for non-deep learning approaches).

train_down <- 
  caret::downSample(x = train[, !(names(train) %in% c("default"))], 
                    y = as.factor(train$default), yname = "default")

base::prop.table(table(train_down$default))

## 
## FALSE  TRUE 
##   0.5   0.5

base::table(train_down$default)

## 
## FALSE  TRUE 
## 17745 17745

12.3 A note on modeling libraries

There are many modeling libraries in R (in fact it is one of the major reasons for its popularity) and they can be very roughly distinguished into

• general modeling libraries comprising many different model functions (sometimes wrappers around specialized model libraries)
• specialized modeling libraries focusing on only one or a certain class of models

We will be using general modeling libraries and in particular the (standard) library shipping with R stats and a popular library caret. Another library often used is rms. The advantage of the general libraries is their easier access via a common API, usually a fairly comprehensive documentation and a relatively large user base which means high chance that someone already did what you like to do.

names(caret::getModelInfo())

##   [1] "ada"                 "AdaBag"              "AdaBoost.M1"        
##   [4] "adaboost"            "amdai"               "ANFIS"              
##   [7] "avNNet"              "awnb"                "awtan"              
##  [10] "bag"                 "bagEarth"            "bagEarthGCV"        
##  [13] "bagFDA"              "bagFDAGCV"           "bam"                
##  [16] "bartMachine"         "bayesglm"            "binda"              
##  [19] "blackboost"          "blasso"              "blassoAveraged"     
##  [22] "bridge"              "brnn"                "BstLm"              
##  [25] "bstSm"               "bstTree"             "C5.0"               
##  [28] "C5.0Cost"            "C5.0Rules"           "C5.0Tree"           
##  [31] "cforest"             "chaid"               "CSimca"             
##  [34] "ctree"               "ctree2"              "cubist"             
##  [37] "dda"                 "deepboost"           "DENFIS"             
##  [40] "dnn"                 "dwdLinear"           "dwdPoly"            
##  [43] "dwdRadial"           "earth"               "elm"                
##  [46] "enet"                "evtree"              "extraTrees"         
##  [49] "fda"                 "FH.GBML"             "FIR.DM"             
##  [52] "foba"                "FRBCS.CHI"           "FRBCS.W"            
##  [55] "FS.HGD"              "gam"                 "gamboost"           
##  [58] "gamLoess"            "gamSpline"           "gaussprLinear"      
##  [61] "gaussprPoly"         "gaussprRadial"       "gbm_h2o"            
##  [64] "gbm"                 "gcvEarth"            "GFS.FR.MOGUL"       
##  [67] "GFS.GCCL"            "GFS.LT.RS"           "GFS.THRIFT"         
##  [70] "glm.nb"              "glm"                 "glmboost"           
##  [73] "glmnet_h2o"          "glmnet"              "glmStepAIC"         
##  [76] "gpls"                "hda"                 "hdda"               
##  [79] "hdrda"               "HYFIS"               "icr"                
##  [82] "J48"                 "JRip"                "kernelpls"          
##  [85] "kknn"                "knn"                 "krlsPoly"           
##  [88] "krlsRadial"          "lars"                "lars2"              
##  [91] "lasso"               "lda"                 "lda2"               
##  [94] "leapBackward"        "leapForward"         "leapSeq"            
##  [97] "Linda"               "lm"                  "lmStepAIC"          
## [100] "LMT"                 "loclda"              "logicBag"           
## [103] "LogitBoost"          "logreg"              "lssvmLinear"        
## [106] "lssvmPoly"           "lssvmRadial"         "lvq"                
## [109] "M5"                  "M5Rules"             "manb"               
## [112] "mda"                 "Mlda"                "mlp"                
## [115] "mlpKerasDecay"       "mlpKerasDecayCost"   "mlpKerasDropout"    
## [118] "mlpKerasDropoutCost" "mlpML"               "mlpSGD"             
## [121] "mlpWeightDecay"      "mlpWeightDecayML"    "monmlp"             
## [124] "msaenet"             "multinom"            "mxnet"              
## [127] "mxnetAdam"           "naive_bayes"         "nb"                 
## [130] "nbDiscrete"          "nbSearch"            "neuralnet"          
## [133] "nnet"                "nnls"                "nodeHarvest"        
## [136] "null"                "OneR"                "ordinalNet"         
## [139] "ORFlog"              "ORFpls"              "ORFridge"           
## [142] "ORFsvm"              "ownn"                "pam"                
## [145] "parRF"               "PART"                "partDSA"            
## [148] "pcaNNet"             "pcr"                 "pda"                
## [151] "pda2"                "penalized"           "PenalizedLDA"       
## [154] "plr"                 "pls"                 "plsRglm"            
## [157] "polr"                "ppr"                 "PRIM"               
## [160] "protoclass"          "pythonKnnReg"        "qda"                
## [163] "QdaCov"              "qrf"                 "qrnn"               
## [166] "randomGLM"           "ranger"              "rbf"                
## [169] "rbfDDA"              "Rborist"             "rda"                
## [172] "regLogistic"         "relaxo"              "rf"                 
## [175] "rFerns"              "RFlda"               "rfRules"            
## [178] "ridge"               "rlda"                "rlm"                
## [181] "rmda"                "rocc"                "rotationForest"     
## [184] "rotationForestCp"    "rpart"               "rpart1SE"           
## [187] "rpart2"              "rpartCost"           "rpartScore"         
## [190] "rqlasso"             "rqnc"                "RRF"                
## [193] "RRFglobal"           "rrlda"               "RSimca"             
## [196] "rvmLinear"           "rvmPoly"             "rvmRadial"          
## [199] "SBC"                 "sda"                 "sdwd"               
## [202] "simpls"              "SLAVE"               "slda"               
## [205] "smda"                "snn"                 "sparseLDA"          
## [208] "spikeslab"           "spls"                "stepLDA"            
## [211] "stepQDA"             "superpc"             "svmBoundrangeString"
## [214] "svmExpoString"       "svmLinear"           "svmLinear2"         
## [217] "svmLinear3"          "svmLinearWeights"    "svmLinearWeights2"  
## [220] "svmPoly"             "svmRadial"           "svmRadialCost"      
## [223] "svmRadialSigma"      "svmRadialWeights"    "svmSpectrumString"  
## [226] "tan"                 "tanSearch"           "treebag"            
## [229] "vbmpRadial"          "vglmAdjCat"          "vglmContRatio"      
## [232] "vglmCumulative"      "widekernelpls"       "WM"                 
## [235] "wsrf"                "xgbLinear"           "xgbTree"            
## [238] "xyf"

12.4 Data preparation for modeling

A number of models require data preparation in advance so here we are going to perform some pre-processing so that all models can work with the data. Note that a few model functions also allow the pre-processing during the training step (e.g. scale and center) but it is computationally more efficient to do at least some pre-processing before.

vars_to_mutate <-
  train_down %>%
  select(which(sapply(.,is.character))) %>%
  names()

vars_to_mutate

##  [1] "term"                "grade"               "emp_length"         
##  [4] "home_ownership"      "verification_status" "pymnt_plan"         
##  [7] "purpose"             "zip_code"            "addr_state"         
## [10] "initial_list_status" "application_type"

train_down <-
  train_down %>%
  mutate_at(.funs = make.names, .vars = vars_to_mutate)
  
test <-
  test %>%
  mutate_at(.funs = make.names, .vars = vars_to_mutate)

vars_to_remove <- c("zip_code")
train_down <- train_down %>% select(-one_of(vars_to_remove))

# train
dummies_train <-
  dummyVars("~.", data = train_down[, !(names(train_down) %in% c("default"))], 
            fullRank = FALSE)

train_down_dummy <-
  train_down %>%
  select(-which(sapply(.,is.character))) %>%
  cbind(predict(dummies_train, newdata = train_down))

# test
dummies_test <-
  dummyVars("~.", data = test[, dummies_train$vars], fullRank = FALSE)

test_dummy <-
  test %>%
  select(one_of(colnames(train_down))) %>%
  select(-which(sapply(.,is.character))) %>%
  cbind(predict(dummies_test, newdata = test))

12.5 Logistic regression

While linear regression is very often used for (continuous) quantitative variables, logistic regression is its counterpart for (discrete) qualitative responses also referred to as categorical. Rather than modeling the outcome directly as default or not, logistic regression models the probability that the response belongs to a particular category. The logistic function will ensure that probabilties are within the range [0, 1]. The model coefficients are estimated via the maximum likelihood method.

The logistic regression model is part of the larger family of generalized linear model function and implemented in many R packages. We will use the stats::glm() function that usually comes with the standard install. A strong contender for any regression modeling is the caret package that we have already used before and will use again.

We will start with a very simple model of one explanatory variable, namely grade, which intuitively should be a good predictor of default as it is a credit rating and therefore an evaluation of the borrower. Note that a call to stats::glm() will not return all model statistics by default but only some attributes of the created model object. We can inspect the object’s attributes with the base::attributes() function. The function base::summary() is a generic function used to produce result summaries of the results of various model fitting functions. The function invokes particular methods which depend on the class of the first argument. In the case of a glm object, it will return model statistics such as

• deviance residuals
• coefficients
• standard error
• z statistic
• p-value
• AIC

and some more. We may also use the summary function on selected model objects, e.g. the coefficients only via summary(model_object)$coef.

Finally, before running the model, the documentation of stats::glm() tells us that if the response is of type factor the first level denotes failure and all others success. We can check and find that this is already the case here.

levels(train_down$default)

## [1] "FALSE" "TRUE"

model_glm_1 <- glm(formula = default ~ grade, 
                   family = binomial(link = "logit"), 
                   data = train_down, na.action = na.exclude)
class(model_glm_1)

## [1] "glm" "lm"

attributes(model_glm_1)

## $names
##  [1] "coefficients"      "residuals"         "fitted.values"    
##  [4] "effects"           "R"                 "rank"             
##  [7] "qr"                "family"            "linear.predictors"
## [10] "deviance"          "aic"               "null.deviance"    
## [13] "iter"              "weights"           "prior.weights"    
## [16] "df.residual"       "df.null"           "y"                
## [19] "converged"         "boundary"          "model"            
## [22] "call"              "formula"           "terms"            
## [25] "data"              "offset"            "control"          
## [28] "method"            "contrasts"         "xlevels"          
## 
## $class
## [1] "glm" "lm"

summary(model_glm_1)

## 
## Call:
## glm(formula = default ~ grade, family = binomial(link = "logit"), 
##     data = train_down, na.action = na.exclude)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8377  -1.1985  -0.0327   1.1564   1.7406  
## 
## Coefficients:
##             Estimate Std. Error z value            Pr(>|z|)    
## (Intercept) -1.26665    0.03911  -32.39 <0.0000000000000002 ***
## gradeB       0.77173    0.04519   17.08 <0.0000000000000002 ***
## gradeC       1.31620    0.04389   29.99 <0.0000000000000002 ***
## gradeD       1.67669    0.04600   36.45 <0.0000000000000002 ***
## gradeE       1.99273    0.05132   38.83 <0.0000000000000002 ***
## gradeF       2.26212    0.06856   33.00 <0.0000000000000002 ***
## gradeG       2.75092    0.12634   21.77 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 49200  on 35489  degrees of freedom
## Residual deviance: 46080  on 35483  degrees of freedom
## AIC: 46094
## 
## Number of Fisher Scoring iterations: 4

summary(model_glm_1)$coef

##               Estimate Std. Error   z value      Pr(>|z|)
## (Intercept) -1.2666450 0.03910963 -32.38704 4.178600e-230
## gradeB       0.7717315 0.04518690  17.07866  2.139941e-65
## gradeC       1.3161994 0.04388804  29.98993 1.327728e-197
## gradeD       1.6766852 0.04600412  36.44641 7.841098e-291
## gradeE       1.9927260 0.05131843  38.83061  0.000000e+00
## gradeF       2.2621234 0.06855512  32.99715 8.925306e-239
## gradeG       2.7509198 0.12633659  21.77453 4.046038e-105

There are a couple of functions that work natively with lm / glm objects created via stats::glm() (taken from Steph Locke: Logistic Regressions):

• coefficients: Extract coefficients
• summary: Output a basic summary
• fitted: Return the predicted values for the training data
• predict: Predict some values for new data
• plot: Produce some basic diagnostic plots
• residuals: Return the errors on predicted values for the training data

Also note that categorical variables will get transformed into dummy variables. A standardization is not necessarily needed, see e.g. Is standardization needed before fitting logistic regression. There are a number of model diagnostics one may perform, see e.g. Diagnostics for logistic regression and pay attention to Frank Harrell’s answer.

12.6 Model evaluation illustration

For the sake of illustration, let’s go once through the model evaluation steps for this very simple model. This workflow should be similar even when the model gets more complex. There are many perspectives one can take when evaluating a model and many packages that provide helper functions to do so. For this illustration, we focus on a fairly standard approach looking at some statistics and explaining them along the way. We will also be plotting an ROC curve to establish its basic functioning.

To evaluate the model performance, we have to use the fitted model to predict the target variable default on unseen (test) data. The function stats::preditc.glm() does exactly that. We set type = "response" which returns the predicted probabilities and thus allows us to evaluate against a chosen threshold. This can be useful in cases where the cost of missclassification is not equal and one outcome should be penalized heavier than another. Also note, that we choose the positive outcome to be TRUE, i.e. a true positive would be correctly identifiying a default. We use the function caret::ConfusionMatrix() to plot a confusion matrix (also called error matrix). The function also returns a number of model statistics (however note that the test data is highly imbalanced as we have not applied any resampling as in the training data - so a number of statistics are meaningless). Here is what the function reports:

• Accuracy: (TP + TN) / (TP + FP + TN + FN)
• 95% CI: 95% confidence interval for accuracy
• No Information Rate: accuracy in case of random guessing (in case of binary outcome the proportion of the majority class)
• P-Value [Acc > NIR]: p-value that accuracy is larger No Information Rate
• Kappa: (observed accuracy - expected accuracy) / (1 - expected accuracy), see Cohen’s kappa in plain English
• Mcnemar’s Test P-Value
• Sensitivity / recall: true positive rate = TP / (TP + FN)
• Specificity: true negative rate = TN / (TN + FP)
• Pos Pred Value: (sensitivity * prevalence) / ((sensitivity * prevalence) + ((1 - specificity) * (1 - prevalence)))
• Neg Pred Value: (specificity * (1-prevalence)) / (((1-sensitivity) * prevalence) + ((specificity) * (1 - prevalence)))
• Prevalence: (TP + FN) / (TP + FP + TN + FN)
• Detection Rate: TP / (TP + FP + TN + FN)
• Detection Prevalence: (TP + FP) / (TP + FP + TN + FN)
• Balanced Accuracy: (sensitivity + specificity) / 2

model_glm_1_pred <- 
  predict.glm(object = model_glm_1, newdata = test, type = "response")
model_pred_t <- function(pred, t) ifelse(pred > t, TRUE, FALSE)
caret::confusionMatrix(data = model_pred_t(model_glm_1_pred, 0.5), 
                       reference = test$default,
                       positive = "TRUE")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction FALSE  TRUE
##      FALSE 79712  1003
##      TRUE  93327  3433
##                                              
##                Accuracy : 0.4685             
##                  95% CI : (0.4662, 0.4708)   
##     No Information Rate : 0.975              
##     P-Value [Acc > NIR] : 1                  
##                                              
##                   Kappa : 0.0211             
##  Mcnemar's Test P-Value : <0.0000000000000002
##                                              
##             Sensitivity : 0.77390            
##             Specificity : 0.46066            
##          Pos Pred Value : 0.03548            
##          Neg Pred Value : 0.98757            
##              Prevalence : 0.02500            
##          Detection Rate : 0.01934            
##    Detection Prevalence : 0.54520            
##       Balanced Accuracy : 0.61728            
##                                              
##        'Positive' Class : TRUE               
## 

Looking at the model statistics, we find a mixed picture:

• we get a fair amount of true defaults right (sensitivity)
• we get a large amount of non-defaults wrong (specificity)
• the Kappa (which should consider class distributions) is very low

We can visualize the classification (predicted probabilities) versus the actual class for the test set and the given threshold. We make use of a plotting function inspired by Illustrated Guide to ROC and AUC with its source code residing on Github. We have slightly adjusted the code to fit our needs.

plot_pred_type_distribution <- function(df, threshold) {
  v <- rep(NA, nrow(df))
  v <- ifelse(df$pred >= threshold & df$actual == 1, "TP", v)
  v <- ifelse(df$pred >= threshold & df$actual == 0, "FP", v)
  v <- ifelse(df$pred < threshold & df$actual == 1, "FN", v)
  v <- ifelse(df$pred < threshold & df$actual == 0, "TN", v)
  
  df$pred_type <- v
  
  ggplot(data = df, aes(x = actual, y = pred)) + 
    geom_violin(fill = rgb(1, 1 ,1, alpha = 0.6), color = NA) + 
    geom_jitter(aes(color = pred_type), alpha = 0.6) +
    geom_hline(yintercept = threshold, color = "red", alpha = 0.6) +
    scale_color_discrete(name = "type") +
    labs(title = sprintf("Threshold at %.2f", threshold))
  }

plot_pred_type_distribution(
  df = tibble::tibble(pred = model_glm_1_pred, actual = test$default), 
  threshold = 0.5)

We can derive a few things from the plot:

• We can clearly see the imbalance of the test dataset with many more observations being FALSE, i.e. non-defaulted
• We can see that predictions tend to cluster more around the center (i.e. the threshold) rather than the top or bottom
• We see that there are only seven distinct (unique) prediction probabilities, namely 0.512386, 0.2198321, 0.3787367, 0.6010975, 0.6739447, 0.8152174, 0.7301686 which corresponds to the number of unique levels in our predictor grade, namely B, C, A, E, F, D, G
• We can see the trade-off of moving the threshold down or up leading to higher sensitivity (lower specificity) or higher specificity (lower sensitivity), respectively.

roc_glm_1 <- pROC::roc(response = test$default, predictor = model_glm_1_pred)
roc_glm_1

## 
## Call:
## roc.default(response = test$default, predictor = model_glm_1_pred)
## 
## Data: model_glm_1_pred in 173039 controls (test$default FALSE) < 4436 cases (test$default TRUE).
## Area under the curve: 0.6641

pROC::plot.roc(x = roc_glm_1, legacy.axes = FALSE, xlim = c(1, 0), asp = NA,
               col = "green", print.auc = FALSE, print.auc.y = .4)

legend(x = "bottomright", legend=c("glm_1 AUC = 0.664"), 
       col = c("green"), lty = 1, cex = 1.0)

12.7 Variable selection

Many articles and book chapters have been written on the topic of variable selection (feature selection) and some authors believe it is one of, if not the most, important topic in model buidling. A few good resources giving an overview, can be found below:

• An Introduction to Statistical Learning - Chapter 6: Linear Model Selection and Regularization
• An Introduction to Variable and Feature Selection
• Variable selection methods in regression: Ignorable problem, outing notable solution

The different approaches can be categorized as follows:

• subset selection
• best subset selection
• stepwise selection (forward/backward/hybrid)
• shrinkage
• ridge regression (l₂-norm)
• lasso (l₁-norm)
• dimension reduction
• principal component analysis (PCA)

There are many trade-offs to be made when choosing one approach over the other but the two most important are statistical soundness and computational efficiency. We will not go into further details of variable selection but rather focus on the model building. We thus simply detrmine a predictor space that seems sensible given prior analysis.

Let’s remind ourselves of variables left for modeling.

full_vars <- colnames(train_down)
full_vars

##  [1] "term"                        "int_rate"                   
##  [3] "grade"                       "emp_length"                 
##  [5] "home_ownership"              "annual_inc"                 
##  [7] "verification_status"         "issue_d"                    
##  [9] "pymnt_plan"                  "purpose"                    
## [11] "addr_state"                  "dti"                        
## [13] "delinq_2yrs"                 "earliest_cr_line"           
## [15] "inq_last_6mths"              "open_acc"                   
## [17] "pub_rec"                     "revol_bal"                  
## [19] "revol_util"                  "initial_list_status"        
## [21] "out_prncp_inv"               "total_rec_late_fee"         
## [23] "collection_recovery_fee"     "last_pymnt_d"               
## [25] "last_pymnt_amnt"             "last_credit_pull_d"         
## [27] "collections_12_mths_ex_med"  "mths_since_last_major_derog"
## [29] "application_type"            "acc_now_delinq"             
## [31] "tot_coll_amt"                "tot_cur_bal"                
## [33] "default"

We determine a subspace and report which variables we ignore.

model_vars <-
  c("term", "grade", "emp_length", "home_ownership", "annual_inc",
    "purpose", "pymnt_plan", "out_prncp_inv", "delinq_2yrs", "default")

ignored_vars <- dplyr::setdiff(full_vars, model_vars)
ignored_vars

##  [1] "int_rate"                    "verification_status"        
##  [3] "issue_d"                     "addr_state"                 
##  [5] "dti"                         "earliest_cr_line"           
##  [7] "inq_last_6mths"              "open_acc"                   
##  [9] "pub_rec"                     "revol_bal"                  
## [11] "revol_util"                  "initial_list_status"        
## [13] "total_rec_late_fee"          "collection_recovery_fee"    
## [15] "last_pymnt_d"                "last_pymnt_amnt"            
## [17] "last_credit_pull_d"          "collections_12_mths_ex_med" 
## [19] "mths_since_last_major_derog" "application_type"           
## [21] "acc_now_delinq"              "tot_coll_amt"               
## [23] "tot_cur_bal"

12.8 A more complex logistic model with caret

The interested reader should check the comprehensive caret documentation on how the API works and for modeling in particular read model training and tuning. We only give a high-level overview here.

The caret package has several functions that attempt to streamline the model building and evaluation process. The train function can be used to

• evaluate, using resampling, the effect of model tuning parameters on performance
• choose the “optimal” model across these parameters
• estimate model performance from a training set

First, a specific model must be chosen. User-defined models can also be created. Second, the corresponding model parameters need to be passed to the function call. The workhorse function is caret::train() but an important input is parameter trControl which also allows specifying a particular resampling approach, e.g. k-fold cross-validation (once or repeated), leave-one-out cross-validation and bootstrap (simple estimation or the 632 rule), etc.

ctrl <- 
  trainControl(method = "repeatedcv", 
               number = 10,
               repeats = 5,
               classProbs = TRUE,
               summaryFunction = twoClassSummary,
               savePredictions = TRUE,
               verboseIter = FALSE)

model_glm_2 <-
  train_down %>%
  select(model_vars) %>%
  mutate(default = as.factor(ifelse(default == TRUE, "yes", "no"))) %>%
  filter(complete.cases(.)) %>%
  train(default ~ ., 
        data = ., 
        method = "glm", 
        family = "binomial",
        metric = "ROC",
        trControl = ctrl)

model_glm_2

## Generalized Linear Model 
## 
## 35488 samples
##     9 predictor
##     2 classes: 'no', 'yes' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 31938, 31938, 31940, 31939, 31940, 31940, ... 
## Resampling results:
## 
##   ROC        Sens       Spec    
##   0.6964721  0.6291356  0.662019

attributes(model_glm_2)

## $names
##  [1] "method"       "modelInfo"    "modelType"    "results"     
##  [5] "pred"         "bestTune"     "call"         "dots"        
##  [9] "metric"       "control"      "finalModel"   "preProcess"  
## [13] "trainingData" "resample"     "resampledCM"  "perfNames"   
## [17] "maximize"     "yLimits"      "times"        "levels"      
## [21] "terms"        "coefnames"    "contrasts"    "xlevels"     
## 
## $class
## [1] "train"         "train.formula"

summary(model_glm_2)

## 
## Call:
## NULL
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3163  -1.0774  -0.2782   1.0635   5.6347  
## 
## Coefficients:
##                                 Estimate     Std. Error z value
## (Intercept)                -1.7381657126   0.1420955332 -12.232
## termX60.months             -0.5453038253   0.0298754606 -18.253
## gradeB                      0.8159277950   0.0461745543  17.671
## gradeC                      1.3693629935   0.0459015145  29.833
## gradeD                      1.7227789216   0.0490499846  35.123
## gradeE                      2.0514999392   0.0562719050  36.457
## gradeF                      2.3642399559   0.0741417266  31.888
## gradeG                      2.8775473195   0.1315002573  21.882
## emp_lengthX..1.year         0.0719719375   0.0634309986   1.135
## emp_lengthX1.year          -0.0231801001   0.0663542839  -0.349
## emp_lengthX10..years       -0.1447186450   0.0541082748  -2.675
## emp_lengthX2.years         -0.0155518922   0.0628511985  -0.247
## emp_lengthX3.years          0.0019357572   0.0640494674   0.030
## emp_lengthX4.years         -0.0416200714   0.0681959839  -0.610
## emp_lengthX5.years         -0.1208526942   0.0671199278  -1.801
## emp_lengthX6.years          0.0685014427   0.0713774190   0.960
## emp_lengthX7.years          0.0169740882   0.0708191548   0.240
## emp_lengthX8.years         -0.0231357344   0.0711988793  -0.325
## emp_lengthX9.years         -0.0454158473   0.0751319778  -0.604
## home_ownershipNONE         -0.5795772695   1.3066383779  -0.444
## home_ownershipOTHER        12.3599134059  91.3562921237   0.135
## home_ownershipOWN           0.1059366445   0.0397823985   2.663
## home_ownershipRENT          0.1738380528   0.0253881630   6.847
## annual_inc                 -0.0000018620   0.0000002708  -6.875
## purposecredit_card          0.0940714394   0.1290711042   0.729
## purposedebt_consolidation   0.1957294997   0.1273405178   1.537
## purposeeducational          1.4894997798   0.4601526184   3.237
## purposehome_improvement     0.2605554222   0.1355011685   1.923
## purposehouse                0.0010870105   0.2087251602   0.005
## purposemajor_purchase       0.1111783795   0.1507309963   0.738
## purposemedical              0.1525955150   0.1664861916   0.917
## purposemoving               0.4431876570   0.1857656345   2.386
## purposeother                0.0979564862   0.1356860122   0.722
## purposerenewable_energy     0.0910160004   0.4132422786   0.220
## purposesmall_business       0.2405034075   0.1572708880   1.529
## purposevacation             0.1203054032   0.1926193435   0.625
## purposewedding             -0.4075063208   0.3202319627  -1.273
## pymnt_plany                11.9809020840 154.7613613910   0.077
## out_prncp_inv               0.0000520918   0.0000017015  30.614
## delinq_2yrs                 0.0865718494   0.0125320528   6.908
##                                       Pr(>|z|)    
## (Intercept)               < 0.0000000000000002 ***
## termX60.months            < 0.0000000000000002 ***
## gradeB                    < 0.0000000000000002 ***
## gradeC                    < 0.0000000000000002 ***
## gradeD                    < 0.0000000000000002 ***
## gradeE                    < 0.0000000000000002 ***
## gradeF                    < 0.0000000000000002 ***
## gradeG                    < 0.0000000000000002 ***
## emp_lengthX..1.year                    0.25652    
## emp_lengthX1.year                      0.72684    
## emp_lengthX10..years                   0.00748 ** 
## emp_lengthX2.years                     0.80457    
## emp_lengthX3.years                     0.97589    
## emp_lengthX4.years                     0.54166    
## emp_lengthX5.years                     0.07177 .  
## emp_lengthX6.years                     0.33720    
## emp_lengthX7.years                     0.81058    
## emp_lengthX8.years                     0.74522    
## emp_lengthX9.years                     0.54552    
## home_ownershipNONE                     0.65736    
## home_ownershipOTHER                    0.89238    
## home_ownershipOWN                      0.00775 ** 
## home_ownershipRENT            0.00000000000753 ***
## annual_inc                    0.00000000000621 ***
## purposecredit_card                     0.46610    
## purposedebt_consolidation              0.12428    
## purposeeducational                     0.00121 ** 
## purposehome_improvement                0.05449 .  
## purposehouse                           0.99584    
## purposemajor_purchase                  0.46076    
## purposemedical                         0.35937    
## purposemoving                          0.01705 *  
## purposeother                           0.47033    
## purposerenewable_energy                0.82568    
## purposesmall_business                  0.12621    
## purposevacation                        0.53225    
## purposewedding                         0.20318    
## pymnt_plany                            0.93829    
## out_prncp_inv             < 0.0000000000000002 ***
## delinq_2yrs                   0.00000000000491 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 49197  on 35487  degrees of freedom
## Residual deviance: 44806  on 35448  degrees of freedom
## AIC: 44886
## 
## Number of Fisher Scoring iterations: 11

predictors(model_glm_2)

##  [1] "termX60.months"            "gradeB"                   
##  [3] "gradeC"                    "gradeD"                   
##  [5] "gradeE"                    "gradeF"                   
##  [7] "gradeG"                    "emp_lengthX..1.year"      
##  [9] "emp_lengthX1.year"         "emp_lengthX10..years"     
## [11] "emp_lengthX2.years"        "emp_lengthX3.years"       
## [13] "emp_lengthX4.years"        "emp_lengthX5.years"       
## [15] "emp_lengthX6.years"        "emp_lengthX7.years"       
## [17] "emp_lengthX8.years"        "emp_lengthX9.years"       
## [19] "home_ownershipNONE"        "home_ownershipOTHER"      
## [21] "home_ownershipOWN"         "home_ownershipRENT"       
## [23] "annual_inc"                "purposecredit_card"       
## [25] "purposedebt_consolidation" "purposeeducational"       
## [27] "purposehome_improvement"   "purposehouse"             
## [29] "purposemajor_purchase"     "purposemedical"           
## [31] "purposemoving"             "purposeother"             
## [33] "purposerenewable_energy"   "purposesmall_business"    
## [35] "purposevacation"           "purposewedding"           
## [37] "pymnt_plany"               "out_prncp_inv"            
## [39] "delinq_2yrs"

varImp(model_glm_2)

## glm variable importance
## 
##   only 20 most important variables shown (out of 39)
## 
##                           Overall
## gradeE                    100.000
## gradeD                     96.340
## gradeF                     87.466
## out_prncp_inv              83.972
## gradeC                     81.827
## gradeG                     60.017
## termX60.months             50.059
## gradeB                     48.462
## delinq_2yrs                18.937
## annual_inc                 18.846
## home_ownershipRENT         18.770
## purposeeducational          8.866
## emp_lengthX10..years        7.323
## home_ownershipOWN           7.291
## purposemoving               6.531
## purposehome_improvement     5.261
## emp_lengthX5.years          4.925
## purposedebt_consolidation   4.202
## purposesmall_business       4.181
## purposewedding              3.477

plot(varImp(model_glm_2))

model_glm_2_pred <- 
  predict(model_glm_2, 
          newdata = test %>% 
            select(model_vars) %>%
            mutate(default = as.factor(ifelse(default == TRUE, 
                                              "yes", "no"))) %>%
            filter(complete.cases(.)), 
          type = "prob")

head(model_glm_2_pred, 3)

##          no       yes
## 1 0.7041954 0.2958046
## 2 0.8783103 0.1216897
## 3 0.6509819 0.3490181

caret::confusionMatrix(
  data = ifelse(model_glm_2_pred[, "yes"] > 0.5, "yes", "no"), 
  reference = as.factor(ifelse(test[complete.cases(test[, model_vars]), 
                                    "default"] == TRUE, "yes", "no")),
  positive = "yes")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction     no    yes
##        no  109156   1482
##        yes  63876   2953
##                                              
##                Accuracy : 0.6317             
##                  95% CI : (0.6295, 0.634)    
##     No Information Rate : 0.975              
##     P-Value [Acc > NIR] : 1                  
##                                              
##                   Kappa : 0.0378             
##  Mcnemar's Test P-Value : <0.0000000000000002
##                                              
##             Sensitivity : 0.66584            
##             Specificity : 0.63084            
##          Pos Pred Value : 0.04419            
##          Neg Pred Value : 0.98660            
##              Prevalence : 0.02499            
##          Detection Rate : 0.01664            
##    Detection Prevalence : 0.37657            
##       Balanced Accuracy : 0.64834            
##                                              
##        'Positive' Class : yes                
## 

temp <- as.factor(ifelse(test[complete.cases(test[, model_vars]),
                                             "default"] == TRUE, "yes", "no"))

roc_glm_2 <- 
  pROC::roc(response = temp, 
            predictor = model_glm_2_pred[, "yes"])

roc_glm_2

## 
## Call:
## roc.default(response = temp, predictor = model_glm_2_pred[, "yes"])
## 
## Data: model_glm_2_pred[, "yes"] in 173032 controls (temp no) < 4435 cases (temp yes).
## Area under the curve: 0.6995

#data = as.data.frame(unclass(train[complete.cases(train), ]))

pROC::plot.roc(x = roc_glm_1, legacy.axes = FALSE, xlim = c(1, 0), asp = NA, 
               col = "green")

pROC::plot.roc(x = roc_glm_2, legacy.axes = FALSE, xlim = c(1, 0), asp = NA,
               add = TRUE, col = "blue")

legend(x = "bottomright", legend=c("glm_1 AUC = 0.664", "glm_2 AUC = 0.7"), 
       col = c("green", "blue"), lty = 1, cex = 1.0)

12.9 A note on computational efficiency (parallelization via cluster)

Some models (and their tuning parameter calibration via resampling) will require significant computational resources both in terms of memory and CPUs. There are a few ways to increase computational efficiency which can be roughly divided into

• decrease training set (without loosing information)
• compute on several cores (parallelization via cluster)
• shift matrix computations to GPU (e.g. via CUDA)
• use compute instances from cloud offerings such as Amazon Web Services, googleCloud or Microsoft Azure

We will only look at parallelization here as it is supported by R and caret out of the box. For details, see The caret package 9: Parallel Processing, caret ml parallel and HOME of the R parallel WIKI! by Tobias Kind. There are a couple of libraries supporting parallel processing in R but one should be aware of the support for different operating systems. On Windows, one should use doParallel which is a parallel backend for the foreach package. Note that the multicore functionality only runs tasks on a single computer, not a cluster of computers. However, you can use the snow functionality to execute on a cluster, using Unix-like operating systems, Windows, or even a combination. It is pointless to use doParallel and parallel on a machine with only one processor with a single core. To get a speed improvement, it must run on a machine with multiple processors, multiple cores, or both. Remember: unless doParallel::registerDoParallel() is called, foreach will not run in parallel. Simply loading the doParallel package is not enough. Note that we can set the parameter caret::trainControl(allowParallel = TRUE) but it is actually the default setting, i.e. caret::train() will automatically run in parallel if it detects a registerd cluster, irrespective of library is used to build it.

To illustrate efficiency gains, we will benchmark an example model training using the previous glm model. We use library (and function with same name) microbenchmark and define to execute the operation several times (to get a stable time estimate) via parameter times. Note that parallizing usually means creating copies of the original data for each worker (core) which requires according memory. Since we only want to illustrate the logic here, we downsize the training data to avoid failures due to insufficient memory.

benchmark_train_data <- 
  train_down[seq(1, nrow(train_down), 2), model_vars] %>%
  mutate(default = as.factor(ifelse(default == TRUE, "yes", "no"))) %>%
  filter(complete.cases(.))
format(utils::object.size(benchmark_train_data), units = "Mb")

## [1] "1.3 Mb"

We see the data has ~ 1 Mb in size.

ctrl <- 
  trainControl(method = "repeatedcv", 
               number = 10,
               repeats = 5,
               classProbs = TRUE,
               summaryFunction = twoClassSummary,
               verboseIter = FALSE,
               allowParallel = TRUE) #default

library(microbenchmark)

glm_nopar <-
  microbenchmark(glm_nopar =
      train(default ~ .,
            data = benchmark_train_data,
            method = "glm",
            family = "binomial",
            metric = "ROC",
            trControl = ctrl),
    times = 5
    )

## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading

## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading

## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading

## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading

## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading

## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading

## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading

## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading

## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading

## Warning in predict.lm(object, newdata, se.fit, scale = 1, type =
## ifelse(type == : prediction from a rank-deficient fit may be misleading

We can look at the resulting run time.

glm_nopar

## Unit: seconds
##       expr      min       lq     mean   median       uq     max neval
##  glm_nopar 27.51427 27.90744 28.27665 28.38591 28.71571 28.8599     5

Now, we will setup a cluster of cores and repeat the benchmarking. We use parallel::detectCores() to establish how many cores are available and use some amount smaller 100% of cores to leave some resources for the system. Usually, we would write a log of the distributed calculation via parallel::makeCluster(outfile = ). It will be opened in append mode as all workers write to it. However, we found issues with running parallel::makeCluster(outfile = ) so we simply run the cluster directly from registerDoParallel() which has no option to define a log file.

# doParallel would also load parallel
library(parallel)
library(doParallel)

## Warning: package 'doParallel' was built under R version 3.4.2

## Loading required package: foreach

## 
## Attaching package: 'foreach'

## The following objects are masked from 'package:purrr':
## 
##     accumulate, when

## Loading required package: iterators

cores_2_use <- floor(0.8 * detectCores())
# cl <- makeCluster(cores_2_use, outfile = "./logs/parallel_log.txt")
# registerDoParallel(cl, cores_2_use)
registerDoParallel(cores = cores_2_use)
getDoParWorkers()

## [1] 6

We can inspect the open connections via base::showConnections().

knitr::kable(base::showConnections())

Now we tried running in parallel but the setup may produce an error on Windows.

Error in e$fun(obj, substitute(ex), parent.frame(), e$data) : unable to find variable "optimismBoot"

This seems to be a bug as noted in SO question Caret train function - unable to find variable “optimismBoot”. The bug seems to be solved in caret development version which can be installed from GitHub via devtools::install_github('topepo/caret/pkg/caret') (note that building caret from source requires RBuildTools which are part of RTools which need to be installed separately, e.g. from Building R for Windows - RStudio may manage the install for you). The issue is also tracked on GitHub under optimismBoot error #706 and may already be resolved in the caret version you are using.

We do not run the install of the development version here again but you should in case you do not have it and wish to test the parallel execution. Before you do, make sure to close the cluster with below code (not executed now as we continue).

parallel::stopCluster(cl)
foreach::registerDoSEQ()

glm_par <-
  microbenchmark(glm_par =
    train(default ~ .,
            data = benchmark_train_data,
            method = "glm",
            family = "binomial",
            metric = "ROC",
            trControl = ctrl),
    times = 5
    )

We make sure to stop the cluster and force R back to sequential execution by using foreach::registerDoSEQ(). In case you used the setup via parallel::makeCluster(), you need to call parallel::stopCluster() first. Otherwise, memory and CPU may be occupied with legacy clusters.

#parallel::stopCluster(cl)
foreach::registerDoSEQ()

Note that we ran into issues when using too much system resources on a Windows 7 workstation as documented in the SO question caret train binary glm fails on parallel cluster via doParallel. Below error was produced.

Error in serialize(data, node$con) : error writing to connection

However, on a Windows 10 machine with fewer cores and development version of caret the issue could not be reproduced.

We have also tried running caret::train() parallel on Linux (Ubuntu) via the doMC library and found that below setting runs just fine and also faster than a sequential execution. We only include the code here for info without execution.

library(doMC)

cores_2_use <- floor(0.8 * parallel::detectCores())
registerDoMC(cores_2_use)

microbenchmark(
  glm_par =
    train(y ~ .,
          data = df,
          method = "glm",
          family = "binomial",
          metric = "ROC",
          trControl = ctrl),
  times = 5)

Finally, let’s compare parallel versus non-parallel execution time. We can see that the parallel setup performed significantly better but we would have expected even more as we are using a multiple of the one-core setup.

glm_nopar

## Unit: seconds
##       expr      min       lq     mean   median       uq     max neval
##  glm_nopar 27.51427 27.90744 28.27665 28.38591 28.71571 28.8599     5

glm_par

## Unit: seconds
##     expr      min       lq     mean  median       uq      max neval
##  glm_par 13.65229 13.79392 17.34424 13.8397 13.89238 31.54293     5

12.10 A note on tuning parameters in caret

When it comes to more complex models in general, often tuning parameters can be used to improve the model. There are different ways to provide the range of parameters to take. In particular, one may use

• pre-defined number (e.g. if domain knowledge already provides and idea)
• use grid search (default in caret::train())
• use random search

In general, the more tuning parameters are available, the higher the computational cost of selecting them. We will show some applications of options below. For further information, see The caret package 5: Model Training and Tuning

12.11 Trees

Trees stratify or segment the predictor space into a number of simple regions. In order to make a prediction for a given observation, we typically use the mean or the mode of the training observations in the region to which it belongs. Trees are easy to interpret and built but as pointed out in “An Introduction to Statistical Learning” they may not compete with the best supervised learning approaches:

Tree-based methods are simple and useful for interpretation. However, they typically are not competitive with the best supervised learning approaches, such as those seen in Chapters 6 and 7, in terms of prediction accuracy. Hence in this chapter we also introduce bagging, random forests, and boosting. Each of these approaches involves producing multiple trees which are then combined to yield a single consensus prediction. We will see that combining a large number of trees can often result in dramatic improvements in prediction accuracy, at the expense of some loss in interpretation.

Following this guidance, we will also look into some alternatives / enhancements of the basic tree model used in CART. Some further resources on trees are found in

• An Introduction to Statistical Learning
• A Complete Tutorial on Tree Based Modeling from Scratch (in R & Python)

As usual, there are many ways to implement tree models in R but we will stay within the caret package and use the API that was already explained earlier.

ctrl <- 
  trainControl(method = "repeatedcv", 
               number = 10,
               repeats = 5,
               classProbs = TRUE,
               summaryFunction = twoClassSummary,
               verboseIter = FALSE,
               allowParallel = TRUE)

# library(rpart)

model_rpart <-
  train_down %>%
  select(model_vars) %>%
  mutate(default = as.factor(ifelse(default == TRUE, "yes", "no"))) %>%
  filter(complete.cases(.)) %>%
  train(default ~ .,
        data = .,
        method = 'rpart',
        metric = "ROC",
        preProc = c("center", "scale"),
        trControl = ctrl)

## Warning in preProcess.default(thresh = 0.95, k = 5, freqCut = 19, uniqueCut
## = 10, : These variables have zero variances: home_ownershipNONE

model_rpart

## CART 
## 
## 35488 samples
##     9 predictor
##     2 classes: 'no', 'yes' 
## 
## Pre-processing: centered (39), scaled (39) 
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 31939, 31938, 31940, 31939, 31940, 31939, ... 
## Resampling results across tuning parameters:
## 
##   cp          ROC        Sens       Spec     
##   0.01445641  0.7018554  0.5572830  0.7674466
##   0.03111086  0.6461547  0.3900115  0.8799842
##   0.26027166  0.5707699  0.6019078  0.5396320
## 
## ROC was used to select the optimal model using  the largest value.
## The final value used for the model was cp = 0.01445641.

ggplot(model_rpart)

plot(model_rpart$finalModel, uniform = TRUE, margin = 0.2)
graphics::text(model_rpart$finalModel)

ctrl <- 
  trainControl(method = "repeatedcv", 
               number = 10,
               repeats = 5,
               classProbs = TRUE,
               summaryFunction = twoClassSummary,
               verboseIter = FALSE,
               allowParallel = TRUE,
               search = "random")

model_rpart <-
  train_down %>%
  select(model_vars) %>%
  mutate(default = as.factor(ifelse(default == TRUE, "yes", "no"))) %>%
  filter(complete.cases(.)) %>%
  train(default ~ .,
        data = .,
        method = 'rpart',
        metric = "ROC",
        preProc = c("center", "scale"),
        trControl = ctrl)

model_rpart

## CART 
## 
## 35488 samples
##     9 predictor
##     2 classes: 'no', 'yes' 
## 
## Pre-processing: centered (39), scaled (39) 
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 31939, 31940, 31940, 31938, 31940, 31940, ... 
## Resampling results across tuning parameters:
## 
##   cp            ROC        Sens       Spec     
##   0.0002254410  0.7418953  0.5806143  0.7904981
##   0.0002536211  0.7433551  0.5725220  0.8003501
##   0.0005072423  0.7462419  0.5532501  0.8225663
## 
## ROC was used to select the optimal model using  the largest value.
## The final value used for the model was cp = 0.0005072423.

plot(model_rpart$finalModel, uniform = TRUE, margin = 0.1)
graphics::text(model_rpart$finalModel, cex = 0.5)

model_rpart_pred <- 
  predict(model_rpart, 
          newdata = test %>% 
            select(model_vars) %>%
            mutate(default = as.factor(ifelse(default == TRUE, 
                                              "yes", "no"))) %>%
            filter(complete.cases(.)), 
          type = "prob")

caret::confusionMatrix(
  data = ifelse(model_rpart_pred[, "yes"] > 0.5, "yes", "no"), 
  reference = as.factor(ifelse(test[complete.cases(test[, model_vars]), 
                                    "default"] == TRUE, "yes", "no")),
  positive = "yes")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    no   yes
##        no  91226   682
##        yes 81806  3753
##                                              
##                Accuracy : 0.5352             
##                  95% CI : (0.5329, 0.5375)   
##     No Information Rate : 0.975              
##     P-Value [Acc > NIR] : 1                  
##                                              
##                   Kappa : 0.0377             
##  Mcnemar's Test P-Value : <0.0000000000000002
##                                              
##             Sensitivity : 0.84622            
##             Specificity : 0.52722            
##          Pos Pred Value : 0.04386            
##          Neg Pred Value : 0.99258            
##              Prevalence : 0.02499            
##          Detection Rate : 0.02115            
##    Detection Prevalence : 0.48211            
##       Balanced Accuracy : 0.68672            
##                                              
##        'Positive' Class : yes                
## 

roc_rpart <- 
  pROC::roc(response = temp, 
            predictor = model_rpart_pred[, "yes"])

roc_rpart

## 
## Call:
## roc.default(response = temp, predictor = model_rpart_pred[, "yes"])
## 
## Data: model_rpart_pred[, "yes"] in 173032 controls (temp no) < 4435 cases (temp yes).
## Area under the curve: 0.7463

pROC::plot.roc(x = roc_glm_1, legacy.axes = FALSE, xlim = c(1, 0), asp = NA,
               col = "green")

pROC::plot.roc(x = roc_glm_2, legacy.axes = FALSE, xlim = c(1, 0), asp = NA,
               add = TRUE, col = "blue")

pROC::plot.roc(x = roc_rpart, legacy.axes = FALSE, xlim = c(1, 0), asp = NA,
               add = TRUE, col = "orange")


legend(x = "bottomright", legend=c("glm_1 AUC = 0.664", "glm_2 AUC = 0.7",
                                   "rpart AUC = 0.74"), 
       col = c("green", "blue", "orange"), lty = 1, cex = 1.0)

# library(randomForest)

ctrl <- 
  trainControl(method = "repeatedcv", 
               number = 5,
               repeats = 1,
               classProbs = TRUE,
               summaryFunction = twoClassSummary,
               verboseIter = FALSE,
               allowParallel = TRUE)

model_rf <-
  train_down %>%
  select(model_vars) %>%
  mutate(default = as.factor(ifelse(default == TRUE, "yes", "no"))) %>%
  filter(complete.cases(.)) %>%
  train(default ~ .,
        data = .,
        method = 'rf',
        ntree = 10,
        metric = "ROC",
        preProc = c("center", "scale"),
        trControl = ctrl)

model_rf

## Random Forest 
## 
## 35488 samples
##     9 predictor
##     2 classes: 'no', 'yes' 
## 
## Pre-processing: centered (39), scaled (39) 
## Resampling: Cross-Validated (5 fold, repeated 1 times) 
## Summary of sample sizes: 28391, 28390, 28390, 28391, 28390 
## Resampling results across tuning parameters:
## 
##   mtry  ROC        Sens       Spec     
##    2    0.6877044  0.5636517  0.7134646
##   20    0.7041407  0.6133559  0.6806064
##   39    0.6982071  0.6123415  0.6655027
## 
## ROC was used to select the optimal model using  the largest value.
## The final value used for the model was mtry = 20.

plot(model_rf$finalModel)

model_rf_pred <- 
  predict(model_rf, 
          newdata = test %>% 
            select(model_vars) %>%
            mutate(default = as.factor(ifelse(default == TRUE, 
                                              "yes", "no"))) %>%
            filter(complete.cases(.)), 
          type = "prob")

caret::confusionMatrix(
  data = ifelse(model_rf_pred[, "yes"] > 0.5, "yes", "no"), 
  reference = as.factor(ifelse(test[complete.cases(test[, model_vars]), 
                                    "default"] == TRUE, "yes", "no")),
  positive = "yes")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction     no    yes
##        no  114693   1723
##        yes  58339   2712
##                                              
##                Accuracy : 0.6616             
##                  95% CI : (0.6594, 0.6638)   
##     No Information Rate : 0.975              
##     P-Value [Acc > NIR] : 1                  
##                                              
##                   Kappa : 0.038              
##  Mcnemar's Test P-Value : <0.0000000000000002
##                                              
##             Sensitivity : 0.61150            
##             Specificity : 0.66284            
##          Pos Pred Value : 0.04442            
##          Neg Pred Value : 0.98520            
##              Prevalence : 0.02499            
##          Detection Rate : 0.01528            
##    Detection Prevalence : 0.34401            
##       Balanced Accuracy : 0.63717            
##                                              
##        'Positive' Class : yes                
## 

roc_rf <- 
  pROC::roc(response = temp, 
            predictor = model_rf_pred[, "yes"])

roc_rf

## 
## Call:
## roc.default(response = temp, predictor = model_rf_pred[, "yes"])
## 
## Data: model_rf_pred[, "yes"] in 173032 controls (temp no) < 4435 cases (temp yes).
## Area under the curve: 0.7021

# library(gbm)

ctrl <- 
  trainControl(method = "repeatedcv", 
               number = 5,
               repeats = 1,
               classProbs = TRUE,
               summaryFunction = twoClassSummary,
               verboseIter = FALSE,
               allowParallel = TRUE)

model_gbm_1 <- 
   train_down %>%
   select(model_vars) %>%
   mutate(default = as.factor(ifelse(default == TRUE, "yes", "no"))) %>%
   filter(complete.cases(.)) %>%
   train(default ~ ., 
         data = ., 
         method = "gbm",
         metric = "ROC",
         trControl = ctrl,
         preProc = c("center", "scale"),
         ## This last option is actually one
         ## for gbm() that passes through
         verbose = FALSE)

model_gbm_1

## Stochastic Gradient Boosting 
## 
## 35488 samples
##     9 predictor
##     2 classes: 'no', 'yes' 
## 
## Pre-processing: centered (39), scaled (39) 
## Resampling: Cross-Validated (5 fold, repeated 1 times) 
## Summary of sample sizes: 28390, 28391, 28390, 28390, 28391 
## Resampling results across tuning parameters:
## 
##   interaction.depth  n.trees  ROC        Sens       Spec     
##   1                   50      0.7218075  0.4989575  0.7923127
##   1                  100      0.7347236  0.5438152  0.7776594
##   1                  150      0.7376056  0.5922232  0.7565241
##   2                   50      0.7351958  0.5467456  0.7754611
##   2                  100      0.7416798  0.6105945  0.7484081
##   2                  150      0.7480446  0.6190476  0.7499302
##   3                   50      0.7387262  0.6076078  0.7507755
##   3                  100      0.7488759  0.6164553  0.7538746
##   3                  150      0.7515752  0.6082840  0.7630048
## 
## Tuning parameter 'shrinkage' was held constant at a value of 0.1
## 
## Tuning parameter 'n.minobsinnode' was held constant at a value of 10
## ROC was used to select the optimal model using  the largest value.
## The final values used for the model were n.trees = 150,
##  interaction.depth = 3, shrinkage = 0.1 and n.minobsinnode = 10.

ggplot(model_gbm_1)

## Warning: Ignoring unknown aesthetics: shape

model_gbm_1_pred <- 
  predict(model_gbm_1, 
          newdata = test %>% 
            select(model_vars) %>%
            mutate(default = as.factor(ifelse(default == TRUE, 
                                              "yes", "no"))) %>%
            filter(complete.cases(.)), 
          type = "prob")

caret::confusionMatrix(
  data = ifelse(model_gbm_1_pred[, "yes"] > 0.5, "yes", "no"), 
  reference = as.factor(ifelse(test[complete.cases(test[, model_vars]), 
                                    "default"] == TRUE, "yes", "no")),
  positive = "yes")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction     no    yes
##        no  104238   1033
##        yes  68794   3402
##                                              
##                Accuracy : 0.6065             
##                  95% CI : (0.6043, 0.6088)   
##     No Information Rate : 0.975              
##     P-Value [Acc > NIR] : 1                  
##                                              
##                   Kappa : 0.0438             
##  Mcnemar's Test P-Value : <0.0000000000000002
##                                              
##             Sensitivity : 0.76708            
##             Specificity : 0.60242            
##          Pos Pred Value : 0.04712            
##          Neg Pred Value : 0.99019            
##              Prevalence : 0.02499            
##          Detection Rate : 0.01917            
##    Detection Prevalence : 0.40681            
##       Balanced Accuracy : 0.68475            
##                                              
##        'Positive' Class : yes                
## 

roc_gbm_1 <- 
  pROC::roc(response = temp, 
            predictor = model_gbm_1_pred[, "yes"])

roc_gbm_1

## 
## Call:
## roc.default(response = temp, predictor = model_gbm_1_pred[, "yes"])
## 
## Data: model_gbm_1_pred[, "yes"] in 173032 controls (temp no) < 4435 cases (temp yes).
## Area under the curve: 0.7516

pROC::plot.roc(x = roc_glm_1, legacy.axes = FALSE, xlim = c(1, 0), asp = NA,
               col = "green")

pROC::plot.roc(x = roc_glm_2, legacy.axes = FALSE, xlim = c(1, 0), asp = NA,
               add = TRUE, col = "blue")

pROC::plot.roc(x = roc_rpart, legacy.axes = FALSE, xlim = c(1, 0), asp = NA,
               add = TRUE, col = "orange")

pROC::plot.roc(x = roc_gbm_1, legacy.axes = FALSE, xlim = c(1, 0), asp = NA,
               add = TRUE, col = "brown")

legend(x = "bottomright", legend=c("glm_1 AUC = 0.664", "glm_2 AUC = 0.7",
                                   "rpart AUC = 0.74", "gbm AUC = 0.753"), 
       col = c("green", "blue", "orange", "brown"), lty = 1, cex = 1.0)