Feature Selection Blog

1.0 Feature Selection Introduction

Feature selection is an important topic that requires in-depth knowledge of the problem domain. Having the right features can help the model perform better. For example, removing the highly correlated attributes can lead to a better model and improve prediction.

In this post, I will explore finding highly correlated variables, Recursive Feature Elimination, stepwise elimination, and Boruta feature selection.

1.1 Plot Histograms

First, I will plot a histogram of the variables in the Pima Indians Diabetes dataset. Various variable exhibit skew distributions and should be considered for transformations.

library("mlbench")

package 㤼㸱mlbench㤼㸲 was built under R version 3.6.3

library("caret")

package 㤼㸱caret㤼㸲 was built under R version 3.6.2Loading required package: lattice

Attaching package: 㤼㸱caret㤼㸲

The following object is masked from 㤼㸱package:purrr㤼㸲:

    lift

# load the data
data(PimaIndiansDiabetes)

1.2 Correlations

Below we will transform the response variable from a factor to numeric. Changing the variable allows us to see any correlations with other variables. I also created a highly correlated variable by combining glucose and mass. As expected glucose and mass are highly correlated with the combination variable HCorrelated.

library("corrplot")

package 㤼㸱corrplot㤼㸲 was built under R version 3.6.1corrplot 0.84 loaded

indians = PimaIndiansDiabetes
indians$diabetes = ifelse(indians$diabetes =="pos",1,0)
indians$HCorrelated =indians$glucose*indians$mass

cor_mx = cor(indians  ,use="pairwise.complete.obs", method = "pearson")
corrplot(cor_mx, method = "color", 
         type = "upper", order = "original", number.cex = .7,
         addCoef.col = "black", # Add coefficient of correlation
         tl.col = "black", tl.srt = 90, # Text label color and rotation
                  # hide correlation coefficient on the principal diagonal
         diag = TRUE)

NA
NA
NA
NA

2.0 Recursive Feature Elimination

2.1 Run the Feature Selection

Recursive Feature Elimination(RFE) builds models with different subsets of a dataset to identify a feature that might not be required. Caret provides a rfe function that facilitates this process.

Below we will load the Pima Indians Diabetes and fit the rfe function. Control was implemented using random forest cross-validated with kfold of 10. The final plot indicates that eight variables have an accuracy of 77.73.

# define the control using a random forest selection function
control = rfeControl(functions=rfFuncs, method="cv", number=10, repeats = 1) # method="cv" , leave out repeats to speed up or method = "repeatedcv", and leave out repeats 
# run the RFE algorithm
set.seed(143)
#
#   NOTE THAT THE PrimaIndians [,9] is a Factor Variable with neg pos levels (factor variable not numeric)- We choose the number of variables we want to see.
#
results = rfe(PimaIndiansDiabetes[,1:8], PimaIndiansDiabetes[,9], sizes=c(1:8), rfeControl=control,
              verbose=FALSE)
# summarize the results
print(results)

Recursive feature selection

Outer resampling method: Cross-Validated (10 fold) 

Resampling performance over subset size:

The top 5 variables (out of 8):
   glucose, mass, age, pregnant, insulin

results

Recursive feature selection

Outer resampling method: Cross-Validated (10 fold) 

Resampling performance over subset size:

The top 5 variables (out of 8):
   glucose, mass, age, pregnant, insulin

2.2 List the predictors in the Order of Choice

# list the chosen features
predictors(results)

[1] "glucose"  "mass"     "age"      "pregnant" "insulin"  "pedigree" "triceps"  "pressure"

2.3 Plot The Results

# plot the results
plot(results, type=c("g", "o"))

Another method relies on fitting a random forest model and identifying variable importance. In this method, variable importance can vary by model.

3.0 Variable Importance

glimpse(PimaIndiansDiabetes)

Rows: 768
Columns: 9
$ pregnant <dbl> 6, 1, 8, 1, 0, 5, 3, 10, 2, 8, 4, 10, 10, 1, 5, 7, 0, 7, 1, 1, 3, 8, 7, 9, 11, 10, 7...
$ glucose  <dbl> 148, 85, 183, 89, 137, 116, 78, 115, 197, 125, 110, 168, 139, 189, 166, 100, 118, 10...
$ pressure <dbl> 72, 66, 64, 66, 40, 74, 50, 0, 70, 96, 92, 74, 80, 60, 72, 0, 84, 74, 30, 70, 88, 84...
$ triceps  <dbl> 35, 29, 0, 23, 35, 0, 32, 0, 45, 0, 0, 0, 0, 23, 19, 0, 47, 0, 38, 30, 41, 0, 0, 35,...
$ insulin  <dbl> 0, 0, 0, 94, 168, 0, 88, 0, 543, 0, 0, 0, 0, 846, 175, 0, 230, 0, 83, 96, 235, 0, 0,...
$ mass     <dbl> 33.6, 26.6, 23.3, 28.1, 43.1, 25.6, 31.0, 35.3, 30.5, 0.0, 37.6, 38.0, 27.1, 30.1, 2...
$ pedigree <dbl> 0.627, 0.351, 0.672, 0.167, 2.288, 0.201, 0.248, 0.134, 0.158, 0.232, 0.191, 0.537, ...
$ age      <dbl> 50, 31, 32, 21, 33, 30, 26, 29, 53, 54, 30, 34, 57, 59, 51, 32, 31, 31, 33, 32, 27, ...
$ diabetes <fct> pos, neg, pos, neg, pos, neg, pos, neg, pos, pos, neg, pos, neg, pos, pos, pos, pos,...

Data=PimaIndiansDiabetes

4.0 StepWise Selection

Stepwise selection is a method that allows for variables to be added or remove in either direction. The model performance is measured in AIC. Akaike information criterion (AIC) estimates the quality of each model relative to each model with the lowest AIC being the best model.

4.1a Backward Selection

The backward procedure begins with a general model that includes all variables and eliminates one variable at a time.

#glm(diabetes ~ ., data=Data, family="binomial")
step(glm(diabetes ~ ., data=Data, family="binomial"),direction="backward")

Start:  AIC=741.45
diabetes ~ pregnant + glucose + pressure + triceps + insulin + 
    mass + pedigree + age

           Df Deviance    AIC
- triceps   1   723.45 739.45
- insulin   1   725.19 741.19
<none>          723.45 741.45
- age       1   725.97 741.97
- pressure  1   729.99 745.99
- pedigree  1   733.78 749.78
- pregnant  1   738.68 754.68
- mass      1   764.22 780.22
- glucose   1   838.37 854.37

Step:  AIC=739.45
diabetes ~ pregnant + glucose + pressure + insulin + mass + pedigree + 
    age

           Df Deviance    AIC
<none>          723.45 739.45
- insulin   1   725.46 739.46
- age       1   725.97 739.97
- pressure  1   730.13 744.13
- pedigree  1   733.92 747.92
- pregnant  1   738.69 752.69
- mass      1   768.77 782.77
- glucose   1   840.87 854.87

Call:  glm(formula = diabetes ~ pregnant + glucose + pressure + insulin + 
    mass + pedigree + age, family = "binomial", data = Data)

Coefficients:
(Intercept)     pregnant      glucose     pressure      insulin         mass     pedigree  
  -8.405136     0.123172     0.035112    -0.013214    -0.001157     0.090089     0.947595  
        age  
   0.014789  

Degrees of Freedom: 767 Total (i.e. Null);  760 Residual
Null Deviance:      993.5 
Residual Deviance: 723.5    AIC: 739.5

4.1b Forward Selection

The forward method begins with a simple model then adds suitable variable one at a time until the best model is obtained.

#glm(diabetes ~ ., data=Data, family="binomial")
step(glm(diabetes ~ ., data=Data, family="binomial"),direction="forward")

Start:  AIC=741.45
diabetes ~ pregnant + glucose + pressure + triceps + insulin + 
    mass + pedigree + age


Call:  glm(formula = diabetes ~ pregnant + glucose + pressure + triceps + 
    insulin + mass + pedigree + age, family = "binomial", data = Data)

Coefficients:
(Intercept)     pregnant      glucose     pressure      triceps      insulin         mass  
  -8.404696     0.123182     0.035164    -0.013296     0.000619    -0.001192     0.089701  
   pedigree          age  
   0.945180     0.014869  

Degrees of Freedom: 767 Total (i.e. Null);  759 Residual
Null Deviance:      993.5 
Residual Deviance: 723.4    AIC: 741.4

4.1c Both Selection

The both method is the combination of backwrad and forward procedures.

#glm(diabetes ~ ., data=Data, family="binomial")
step(glm(diabetes ~ ., data=Data, family="binomial"),direction="both")

Start:  AIC=741.45
diabetes ~ pregnant + glucose + pressure + triceps + insulin + 
    mass + pedigree + age

           Df Deviance    AIC
- triceps   1   723.45 739.45
- insulin   1   725.19 741.19
<none>          723.45 741.45
- age       1   725.97 741.97
- pressure  1   729.99 745.99
- pedigree  1   733.78 749.78
- pregnant  1   738.68 754.68
- mass      1   764.22 780.22
- glucose   1   838.37 854.37

Step:  AIC=739.45
diabetes ~ pregnant + glucose + pressure + insulin + mass + pedigree + 
    age

           Df Deviance    AIC
<none>          723.45 739.45
- insulin   1   725.46 739.46
- age       1   725.97 739.97
+ triceps   1   723.45 741.45
- pressure  1   730.13 744.13
- pedigree  1   733.92 747.92
- pregnant  1   738.69 752.69
- mass      1   768.77 782.77
- glucose   1   840.87 854.87

Call:  glm(formula = diabetes ~ pregnant + glucose + pressure + insulin + 
    mass + pedigree + age, family = "binomial", data = Data)

Coefficients:
(Intercept)     pregnant      glucose     pressure      insulin         mass     pedigree  
  -8.405136     0.123172     0.035112    -0.013214    -0.001157     0.090089     0.947595  
        age  
   0.014789  

Degrees of Freedom: 767 Total (i.e. Null);  760 Residual
Null Deviance:      993.5 
Residual Deviance: 723.5    AIC: 739.5

Seems that the model without Triceps is Best

5.0 Boruta

Boruta is a feature ranking and selection algorithm based on random forest algorithm. The advantages of using this package are the ease of variables selection and the ability to adjust variable selection.

Below I fitted the Boruta function with the dataset for evaluation.

library('Boruta')
set.seed(143)
boruta_output = Boruta(diabetes ~ ., data=na.omit(Data), doTrace=0) 

5.1 Significant Variables

The significant variables can be extracted from the selection. Tentative variables are variables that can be dropped or kept.

Significant_vars = getSelectedAttributes(boruta_output, withTentative = TRUE)
Significant_vars

[1] "pregnant" "glucose"  "pressure" "triceps"  "insulin"  "mass"     "pedigree" "age"     

Seems all are significant

Boruta has a method for making the selecting tentative variable for the user.

roughFixMod = TentativeRoughFix(boruta_output)
boruta_signif = getSelectedAttributes(roughFixMod)
boruta_signif

[1] "pregnant" "glucose"  "pressure" "triceps"  "insulin"  "mass"     "pedigree" "age"     

The importance of variables can be shown by the below method with being the most important variable.

5.2 Importance of Variables

The importance of variables can be shown by the below method with vari being the most important variable.

# Variable Importance Scores
imps = attStats(roughFixMod)
imps2 = imps[imps$decision != 'Rejected', c('meanImp', 'decision')]
head(imps2[order(-imps2$meanImp), ],10)  # descending sort

NA
NA

5.3 Plotting Importance of Variables

plot(boruta_output, cex.axis=.7, las=2, xlab="", main="Variable Importance")  

6.0 Information Value and Weights of Evidence (Categorical Variables)

The Information Value can be used to judge how important a given categorical variable is in explaining the binary Y variable. It goes well with logistic regression and other classification models that can model binary variables.

Let’s try to find out how important the categorical variables are in predicting if an individual will earn >50k from the ‘adult.csv’ dataset. Just run the code below to import the dataset.

Alright, let’s now find the information value for the categorical variables in the inputData.

Here is what the quantum of Information Value means:

Less than 0.02, then the predictor is not useful for modeling (separating the Goods from the Bads) 0.02 to 0.1, then the predictor has only a weak relationship. 0.1 to 0.3, then the predictor has a medium strength relationship. 0.3 or higher, then the predictor has a strong relationship. That was about IV. Then what is Weight of Evidence?

Weights of evidence can be useful to find out how important a given categorical variable is in explaining the ‘events’ (called ‘Goods’ in below table.)

WOE = ln(%good of all good/%bad of all bad)

Here is what the quantum of Information Value means:

Less than 0.02, then the predictor is not useful for modeling (separating the Goods from the Bads) 0.02 to 0.1, then the predictor has only a weak relationship. 0.1 to 0.3, then the predictor has a medium strength relationship. 0.3 or higher, then the predictor has a strong relationship. That was about IV. Then what is Weight of Evidence?

Weights of evidence can be useful to find out how important a given categorical variable is in explaining the ‘events’ (called ‘Goods’ in below table.)

# The ‘Information Value’ of the categorical variable can then be derived from the respective WOE values.
# 
# IV?=?(perc good of all goods?perc bad of all bads)?*?WOE
# 
# The ‘WOETable’ below given the computation in more detail.

WOETable(X=inputData[, 'workclass'], Y=inputData$income)

The total IV of a variable is the sum of IV�s of its categories.

#—————————————————————————————————————- ## TUTORIAL USING OPTIMAL BINNING () # ——————————————————————————————————————-

7.0 WOE FROM Information Pcakage

7.1 Woe Continuous Variables and Factor Variables

summary(mydata)

       X             admit             gre             gpa        rank   
 Min.   :  1.0   Min.   :0.0000   Min.   :220.0   Min.   :2.260   1: 61  
 1st Qu.:100.8   1st Qu.:0.0000   1st Qu.:520.0   1st Qu.:3.130   2:151  
 Median :200.5   Median :0.0000   Median :580.0   Median :3.395   3:121  
 Mean   :200.5   Mean   :0.3175   Mean   :587.7   Mean   :3.390   4: 67  
 3rd Qu.:300.2   3rd Qu.:1.0000   3rd Qu.:660.0   3rd Qu.:3.670          
 Max.   :400.0   Max.   :1.0000   Max.   :800.0   Max.   :4.000          

7.2 Make an independent Variable a Factor

It is important to note here the number of bins for ‘rank’ variable. Since it is a categorical variable, the number of bins would be according to unique values of the factor variable. The parameter bins=10 does not work for a factor variable.

mydata$rank=factor(mydata$rank)

7.3 Compute The Info Value

7.4 Print the IV Tables

To get WOE table for variable gre, you need to call Tables list from IV list.

We can do this for any of the variables

print(IV$Tables$gre, row.names=FALSE)

print(IV$Tables$gpa, row.names=FALSE)

print(IV$Tables$rank, row.names=FALSE)

print(IV$Tables$X, row.names=FALSE)

7.5 Put The Tables if Dataframes

To get WOE table for variable gre, you need to call Tables list from IV list.

We can do this for any of the variables

gre = data.frame(IV$Tables$gre)

gpa = data.frame(IV$Tables$gpa)

X = data.frame(IV$Tables$X)

rank = data.frame(IV$Tables$rank)

7.6 Plot Woe Scores For 1 Variable

We can plot 1 at a time

plot_infotables(IV, "gre")

7.7 Plot Woe Scores For Many Variable

We can plot many at a time

plot_infotables(IV, IV$Summary$Variable[1:4], same_scale=FALSE)

8.0 WOE Bimmnning

This package generates, visualizes, tabulates and deploys a supervised weight of evidence (WOE) binning of variables.

Details

This package generates, visualizes, tabulates and deploys a supervised weight of evidence (WOE) binning of variables.

The package woeBinning automates the process of binning of numeric variables and factors with respect to a dichotomous target variable. Additionally, it visualizes the realized binning solution, tabulates it and deploys it to (new) data. All functions can be used with single variables or an entire data frame.

8.1 Binning

woe.binning generates a supervised fine and coarse classing of numeric variables and factors.

woe.tree.binning generates a supervised tree-like segmentation of numeric variables and factors.

woe.binning.plot visualizes the binning solution generated and saved via woe.binning or woe.tree.binning.

woe.binning.table tabulates the binning solution generated and saved via woe.binning or woe.tree.binning.

woe.binning.deploy deploys the binning solution generated and saved via woe.binning or woe.tree.binning to (new) data.

References

Siddiqi, N. 2006: Credit Risk Scorecards: Developing and Implementing Intelligent Credit Scoring. Hoboken, New Jersey: John Wiley & Sons.

Anderson, R. 2007: The Credit Scoring Toolkit: Theory and Practice for Retail Credit Risk Management and Decision Automation. Oxford / New York: Oxford University Press.

8.1.1 Binning of Numeric Variables

Numeric variables (continuous and ordinal) are binned by merging initial classes with similar frequencies. The number of initial bins results from the min.perc.total parameter: min.perc.total will result in trunc(1/min.perc.total) initial bins, whereby trunc is needed to guarantee bins with similar frequencies. For example min.perc.total=0.07 will cause trunc(14.3)=14 initial classes. Next, if min.perc.class>0, bins with sparse target classes will be merged with the next upper bin, and in case of the last bin with the next lower one. NAs have their own bin and will not be merged with others. Finally nearby bins with most similar weight of evidence (WOE) values are joined step by step until the information value (IV) decreases more than specified by a percentage value (stop.limit parameter) or until two bins are reached.

8.1.2 Binning of Factors

Factors (categorical variables) are binned by merging factor levels. As a start sparse levels (defined via the min.perc.total and min.perc.class parameters) are merged to a ‘miscellaneous’ level: if possible, respective levels (including sparse NAs) are bundled as ‘misc. level pos.’ (associated with positive WOE values), respectively as ‘misc. level neg.’ (associated with negative WOE values). In case a misc. level contains only NAs it will be named ‘Missing’. Afterwards levels with similar WOE values are joined step by step until the information value (IV) decreases more than specified by a percentage value (stop.limit parameter) or until two bins are reached.

8.1.3 Adjustment of 0 Frequencies

In case the crosstab of the bins with the target classes contains frequencies = 0 the column percentages are adjusted to be able to compute the WOE and IV values: the offset 0.0001 (=0.01%) is added to each column percentage cell and the column percentages are recomputed then. This allows considering bins associated with one target class only, but may cause extreme WOE values for these bins. If a correction is not appropriate choose min.perc.class>0; bins with sparse target classes will be merged then before computing any WOE or IV value.

8.1.4 See Also

Other binning functions: woe.tree.binning

8.1.5 Binning Duration on Tree

woe.tree.binning generates a supervised tree-like segmentation of numeric variables and factors with respect to a dichotomous target variable. Its parameters provide flexibility in finding a binning that fits specific data characteristics and practical needs.

https://rdrr.io/cran/woeBinning/man/woe.binning.html # Shift Click to open link

8.2 Examples

9.0 Post Processing on The WoE Above

9.1 Correlation Between Categporical and Numeric Variables

The WoE is now a continuous unit so we can tell the correlation between contuous and categorical Units

We conduct Correlation using the data frame compiled in the WoE Binning Above

9.1a Lets glimpse the DF and get a Sample

glimpse(df.with.binned.vars.added)

Rows: 1,000
Columns: 13
$ creditability                        <fct> good, bad, good, good, bad, good, good, good, good, bad,...
$ credit.amount                        <dbl> 1169, 5951, 2096, 7882, 4870, 9055, 2835, 6948, 3059, 52...
$ duration.in.month                    <dbl> 6, 48, 12, 42, 24, 36, 24, 36, 12, 30, 12, 48, 12, 24, 1...
$ savings.account.and.bonds            <fct> unknown/ no savings account, ... < 100 DM, ... < 100 DM,...
$ purpose                              <fct> radio/television, radio/television, education, furniture...
$ duration.in.month.binned             <fct> "(-Inf,6]", "(30, Inf]", "(6,15]", "(30, Inf]", "(15,30]...
$ woe.duration.in.month.binned         <dbl> -124.59370, 76.63288, -36.53869, 76.63288, 10.83411, 76....
$ savings.account.and.bonds.binned     <fct> unknown/ no savings account + 500 <= ... < 1000 DM, ... ...
$ woe.savings.account.and.bonds.binned <dbl> -70.47080, 25.24534, 25.24534, 25.24534, 25.24534, -70.4...
$ purpose.binned                       <fct> radio/television + car (used), radio/television + car (u...
$ woe.purpose.binned                   <dbl> -50.02820, -50.02820, 27.99201, 27.99201, 27.99201, 27.9...
$ credit.amount.binned                 <fct> "(-Inf,5969.95]", "(-Inf,5969.95]", "(-Inf,5969.95]", "(...
$ woe.credit.amount.binned             <dbl> -13.83897, -13.83897, -13.83897, 44.18328, -13.83897, 44...

library("corrplot")

library(dplyr)



df=df.with.binned.vars.added%>%select(creditability,contains("woe"))

df$creditability=ifelse(df$creditability=="bad",1,0)




cor_mat_df = cor(df  ,use="pairwise.complete.obs", method = "pearson")
corrplot(cor_mat_df, method = "color", 
         type = "upper", order = "original", number.cex = .7,
         addCoef.col = "black", # Add coefficient of correlation
         tl.col = "black", tl.srt = 90, # Text label color and rotation
                  # hide correlation coefficient on the principal diagonal
         diag = TRUE)

NA
NA
NA

Result is by correlation with 1/0 it seems that there is an order of importance , but lets check this by variable importance and 2 ways, Via WOE and By Factors , ie the binning

9.2 Variable importance By Woe

We will use RFE to gauge variable importance of the above frame

library('Boruta')
set.seed(143)
boruta_output_woe = Boruta(creditability ~ ., data=na.omit(df), doTrace=0) 

Significant_vars = getSelectedAttributes(boruta_output_woe, withTentative = TRUE)
Significant_vars

[1] "woe.duration.in.month.binned"         "woe.savings.account.and.bonds.binned"
[3] "woe.purpose.binned"                   "woe.credit.amount.binned"            

# Neede Roughfix conversionj
roughFixMod <- TentativeRoughFix(boruta_output_woe)

There are no Tentative attributes! Returning original object.

boruta_signif <- getSelectedAttributes(roughFixMod)
print(boruta_signif)

[1] "woe.duration.in.month.binned"         "woe.savings.account.and.bonds.binned"
[3] "woe.purpose.binned"                   "woe.credit.amount.binned"            

plot(boruta_output_woe, cex.axis=.7, las=2, xlab="", main="Variable Importance Woe")  

imps = attStats(roughFixMod)
imps2 = imps[imps$decision != 'Rejected', c('meanImp', 'decision')]
head(imps2[order(-imps2$meanImp), ],10)  # descending sort

# Agrees with Correlation

#
# Now we look at the binning bFactors

df_factors=df.with.binned.vars.added%>%select(creditability,contains("binned"))%>%select_if(is.factor)

glimpse(df_factors)

Rows: 1,000
Columns: 5
$ creditability                    <fct> good, bad, good, good, bad, good, good, good, good, bad, bad...
$ duration.in.month.binned         <fct> "(-Inf,6]", "(30, Inf]", "(6,15]", "(30, Inf]", "(15,30]", "...
$ savings.account.and.bonds.binned <fct> unknown/ no savings account + 500 <= ... < 1000 DM, ... < 10...
$ purpose.binned                   <fct> radio/television + car (used), radio/television + car (used)...
$ credit.amount.binned             <fct> "(-Inf,5969.95]", "(-Inf,5969.95]", "(-Inf,5969.95]", "(5969...

set.seed(143)
boruta_output_binned = Boruta(creditability ~ ., data=na.omit(df_factors), doTrace=0) 

Significant_vars = getSelectedAttributes(boruta_output_binned, withTentative = TRUE)
Significant_vars

[1] "duration.in.month.binned"         "savings.account.and.bonds.binned"
[3] "purpose.binned"                   "credit.amount.binned"            

# Neede Roughfix conversionj
roughFixMod <- TentativeRoughFix(boruta_output_binned)

There are no Tentative attributes! Returning original object.

boruta_signif <- getSelectedAttributes(roughFixMod)
print(boruta_signif)

[1] "duration.in.month.binned"         "savings.account.and.bonds.binned"
[3] "purpose.binned"                   "credit.amount.binned"            

imps = attStats(roughFixMod)
imps2 = imps[imps$decision != 'Rejected', c('meanImp', 'decision')]
head(imps2[order(-imps2$meanImp), ],10)  # descending sort

plot(boruta_output_binned, cex.axis=.7, las=2, xlab="", main="Variable Importance Binned")  

NA
NA

9.3 Decision Trees

We Now use the Binns to do the decision Trees

Use rpart

rpart.plot(tree)

Bad 'data' field in model 'call' (expected a data.frame or a matrix).
To silence this warning:
    Call rpart.plot with roundint=FALSE,
    or rebuild the rpart model with model=TRUE.

10.0 Excerpts From Credit Risk Modelling

10.1a Creditability Vs Duration Of Loan

cross_CredVsDuration=CrossTable(df_factors$duration.in.month.binned,df_factors$creditability,, prop.r = TRUE,
prop.c = FALSE, prop.t = FALSE, prop.chisq = TRUE,chisq = TRUE)

 
   Cell Contents
|-------------------------|
|                       N |
| Chi-square contribution |
|           N / Row Total |
|-------------------------|

 
Total Observations in Table:  1000 

 
                                    | df_factors$creditability 
df_factors$duration.in.month.binned |       bad |      good | Row Total | 
------------------------------------|-----------|-----------|-----------|
                           (-Inf,6] |         9 |        73 |        82 | 
                                    |     9.893 |     4.240 |           | 
                                    |     0.110 |     0.890 |     0.082 | 
------------------------------------|-----------|-----------|-----------|
                             (6,15] |        80 |       269 |       349 | 
                                    |     5.827 |     2.497 |           | 
                                    |     0.229 |     0.771 |     0.349 | 
------------------------------------|-----------|-----------|-----------|
                            (15,30] |       128 |       268 |       396 | 
                                    |     0.712 |     0.305 |           | 
                                    |     0.323 |     0.677 |     0.396 | 
------------------------------------|-----------|-----------|-----------|
                          (30, Inf] |        83 |        90 |       173 | 
                                    |    18.636 |     7.987 |           | 
                                    |     0.480 |     0.520 |     0.173 | 
------------------------------------|-----------|-----------|-----------|
                       Column Total |       300 |       700 |      1000 | 
------------------------------------|-----------|-----------|-----------|

 
Statistics for All Table Factors


Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  50.09743     d.f. =  3     p =  7.616443e-11 


 

plot(cross_CredVsDuration$prop.tbl,col=c("red","green"),main="Creditability Vs Duration",xlab = "Duration",ylab = "Creditability")

Conclusion: The longer the duration the worse the creditability

10.1b Creditability Vs Amount Of Loan

CrossTable(df_factors$credit.amount.binned,df_factors$creditability, prop.r = TRUE,
prop.c = FALSE, prop.t = FALSE, prop.chisq = FALSE)

 
   Cell Contents
|-------------------------|
|                       N |
|           N / Row Total |
|-------------------------|

 
Total Observations in Table:  1000 

 
                                | df_factors$creditability 
df_factors$credit.amount.binned |       bad |      good | Row Total | 
--------------------------------|-----------|-----------|-----------|
                 (-Inf,5969.95] |       231 |       619 |       850 | 
                                |     0.272 |     0.728 |     0.850 | 
--------------------------------|-----------|-----------|-----------|
               (5969.95,9162.7] |        40 |        60 |       100 | 
                                |     0.400 |     0.600 |     0.100 | 
--------------------------------|-----------|-----------|-----------|
                  (9162.7, Inf] |        29 |        21 |        50 | 
                                |     0.580 |     0.420 |     0.050 | 
--------------------------------|-----------|-----------|-----------|
                   Column Total |       300 |       700 |      1000 | 
--------------------------------|-----------|-----------|-----------|

 

cross_CredVsAmount=CrossTable(df_factors$credit.amount.binned,df_factors$creditability,, prop.r = TRUE,
prop.c = FALSE, prop.t = FALSE, prop.chisq = TRUE,chisq = TRUE)

 
   Cell Contents
|-------------------------|
|                       N |
| Chi-square contribution |
|           N / Row Total |
|-------------------------|

 
Total Observations in Table:  1000 

 
                                | df_factors$creditability 
df_factors$credit.amount.binned |       bad |      good | Row Total | 
--------------------------------|-----------|-----------|-----------|
                 (-Inf,5969.95] |       231 |       619 |       850 | 
                                |     2.259 |     0.968 |           | 
                                |     0.272 |     0.728 |     0.850 | 
--------------------------------|-----------|-----------|-----------|
               (5969.95,9162.7] |        40 |        60 |       100 | 
                                |     3.333 |     1.429 |           | 
                                |     0.400 |     0.600 |     0.100 | 
--------------------------------|-----------|-----------|-----------|
                  (9162.7, Inf] |        29 |        21 |        50 | 
                                |    13.067 |     5.600 |           | 
                                |     0.580 |     0.420 |     0.050 | 
--------------------------------|-----------|-----------|-----------|
                   Column Total |       300 |       700 |      1000 | 
--------------------------------|-----------|-----------|-----------|

 
Statistics for All Table Factors


Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  26.65546     d.f. =  2     p =  1.628696e-06 


 

plot(cross_CredVsAmount$prop.tbl,col=c("red","green"),main="Creditability Vs Amount",xlab = "Amount_Bin",ylab = "Creditability")

10.1c Creditability Vs Savings

library(gmodels)
library(gridExtra)

## Lets Plot Some CrossTables - Normally categorical Data to Show Relationship in a Table
#  Tpo help Understand the Tree
#
CrossTable(df_factors$savings.account.and.bonds.binned,df_factors$creditability, prop.r = TRUE,
prop.c = FALSE, prop.t = FALSE, prop.chisq = FALSE)

 
   Cell Contents
|-------------------------|
|                       N |
|           N / Row Total |
|-------------------------|

 
Total Observations in Table:  1000 

 
                                                   | df_factors$creditability 
       df_factors$savings.account.and.bonds.binned |       bad |      good | Row Total | 
---------------------------------------------------|-----------|-----------|-----------|
                ... < 100 DM + 100 <= ... < 500 DM |       251 |       455 |       706 | 
                                                   |     0.356 |     0.644 |     0.706 | 
---------------------------------------------------|-----------|-----------|-----------|
                                  misc. level neg. |         6 |        42 |        48 | 
                                                   |     0.125 |     0.875 |     0.048 | 
---------------------------------------------------|-----------|-----------|-----------|
unknown/ no savings account + 500 <= ... < 1000 DM |        43 |       203 |       246 | 
                                                   |     0.175 |     0.825 |     0.246 | 
---------------------------------------------------|-----------|-----------|-----------|
                                      Column Total |       300 |       700 |      1000 | 
---------------------------------------------------|-----------|-----------|-----------|

 

cross_CredVsSavings=CrossTable(df_factors$savings.account.and.bonds.binned,df_factors$creditability,, prop.r = TRUE,
prop.c = FALSE, prop.t = FALSE, prop.chisq = TRUE,chisq = TRUE)

 
   Cell Contents
|-------------------------|
|                       N |
| Chi-square contribution |
|           N / Row Total |
|-------------------------|

 
Total Observations in Table:  1000 

 
                                                   | df_factors$creditability 
       df_factors$savings.account.and.bonds.binned |       bad |      good | Row Total | 
---------------------------------------------------|-----------|-----------|-----------|
                ... < 100 DM + 100 <= ... < 500 DM |       251 |       455 |       706 | 
                                                   |     7.255 |     3.109 |           | 
                                                   |     0.356 |     0.644 |     0.706 | 
---------------------------------------------------|-----------|-----------|-----------|
                                  misc. level neg. |         6 |        42 |        48 | 
                                                   |     4.900 |     2.100 |           | 
                                                   |     0.125 |     0.875 |     0.048 | 
---------------------------------------------------|-----------|-----------|-----------|
unknown/ no savings account + 500 <= ... < 1000 DM |        43 |       203 |       246 | 
                                                   |    12.854 |     5.509 |           | 
                                                   |     0.175 |     0.825 |     0.246 | 
---------------------------------------------------|-----------|-----------|-----------|
                                      Column Total |       300 |       700 |      1000 | 
---------------------------------------------------|-----------|-----------|-----------|

 
Statistics for All Table Factors


Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  35.72764     d.f. =  2     p =  1.745187e-08 


 

plot(cross_CredVsSavings$prop.tbl,col=c("red","green"),main="Creditability Vs Savings",xlab = "Savings_Bin",ylab = "Creditability")

10.1d Creditability Vs Purpose

library(gmodels)
library(gridExtra)

## Lets Plot Some CrossTables - Normally categorical Data to Show Relationship in a Table
#  Tpo help Understand the Tree
#
CrossTable(df_factors$purpose.binned,df_factors$creditability, prop.r = TRUE,
prop.c = FALSE, prop.t = FALSE, prop.chisq = FALSE)

 
   Cell Contents
|-------------------------|
|                       N |
|           N / Row Total |
|-------------------------|

 
Total Observations in Table:  1000 

 
                                                                          | df_factors$creditability 
                                                df_factors$purpose.binned |       bad |      good | Row Total | 
--------------------------------------------------------------------------|-----------|-----------|-----------|
education + car (new) + misc. level pos. + business + furniture/equipment |       220 |       388 |       608 | 
                                                                          |     0.362 |     0.638 |     0.608 | 
--------------------------------------------------------------------------|-----------|-----------|-----------|
                                            radio/television + car (used) |        79 |       304 |       383 | 
                                                                          |     0.206 |     0.794 |     0.383 | 
--------------------------------------------------------------------------|-----------|-----------|-----------|
                                                         misc. level neg. |         1 |         8 |         9 | 
                                                                          |     0.111 |     0.889 |     0.009 | 
--------------------------------------------------------------------------|-----------|-----------|-----------|
                                                             Column Total |       300 |       700 |      1000 | 
--------------------------------------------------------------------------|-----------|-----------|-----------|

 

cross_CredVsPurpose=CrossTable(df_factors$purpose.binned,df_factors$creditability,, prop.r = TRUE,
prop.c = FALSE, prop.t = FALSE, prop.chisq = TRUE,chisq = TRUE)

Chi-squared approximation may be incorrect

 
   Cell Contents
|-------------------------|
|                       N |
| Chi-square contribution |
|           N / Row Total |
|-------------------------|

 
Total Observations in Table:  1000 

 
                                                                          | df_factors$creditability 
                                                df_factors$purpose.binned |       bad |      good | Row Total | 
--------------------------------------------------------------------------|-----------|-----------|-----------|
education + car (new) + misc. level pos. + business + furniture/equipment |       220 |       388 |       608 | 
                                                                          |     7.751 |     3.322 |           | 
                                                                          |     0.362 |     0.638 |     0.608 | 
--------------------------------------------------------------------------|-----------|-----------|-----------|
                                            radio/television + car (used) |        79 |       304 |       383 | 
                                                                          |    11.217 |     4.807 |           | 
                                                                          |     0.206 |     0.794 |     0.383 | 
--------------------------------------------------------------------------|-----------|-----------|-----------|
                                                         misc. level neg. |         1 |         8 |         9 | 
                                                                          |     1.070 |     0.459 |           | 
                                                                          |     0.111 |     0.889 |     0.009 | 
--------------------------------------------------------------------------|-----------|-----------|-----------|
                                                             Column Total |       300 |       700 |      1000 | 
--------------------------------------------------------------------------|-----------|-----------|-----------|

 
Statistics for All Table Factors


Pearson's Chi-squared test 
------------------------------------------------------------
Chi^2 =  28.62578     d.f. =  2     p =  6.081227e-07 


 

plot(cross_CredVsPurpose$prop.tbl,col=c("red","green"),main="Creditability Vs Purpose",xlab = "Savings_Bin",ylab = "Creditability")

10.2a Decision Trees Reference

library(rpart)
library(rpart.plot)

set.seed(123)
tree <- rpart(creditability ~ ., data=na.omit(df_factors),control = rpart.control(cp = 0.0001),model = TRUE)

tree

n= 1000 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

1) root 1000 300 good (0.3000000 0.7000000)  
  2) savings.account.and.bonds.binned=... < 100 DM + 100 <= ... < 500 DM 706 251 good (0.3555241 0.6444759)  
    4) duration.in.month.binned=(30, Inf] 125  51 bad (0.5920000 0.4080000) *
    5) duration.in.month.binned=(-Inf,6],(6,15],(15,30] 581 177 good (0.3046472 0.6953528) *
  3) savings.account.and.bonds.binned=misc. level neg.,unknown/ no savings account + 500 <= ... < 1000 DM 294  49 good (0.1666667 0.8333333) *

rpart.plot(tree)

NA
NA

10.2b Adjust Priors and Pruning

library(rpart)
library(rpart.plot)

set.seed(123)
tree <- rpart(creditability ~ ., data=na.omit(df_factors),control = rpart.control(cp = 0.000001),model = TRUE,parms = list(prior = c(0.7, 0.3)))

tree

n= 1000 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

 1) root 1000 300.000000 bad (0.7000000 0.3000000)  
   2) savings.account.and.bonds.binned=... < 100 DM + 100 <= ... < 500 DM 706 195.000000 bad (0.7502135 0.2497865)  
     4) duration.in.month.binned=(15,30],(30, Inf] 400  94.714290 bad (0.8151487 0.1848513)  
       8) purpose.binned=education + car (new) + misc. level pos. + business + furniture/equipment 252  51.857140 bad (0.8549547 0.1450453) *
       9) purpose.binned=radio/television + car (used),misc. level neg. 148  42.857140 bad (0.7232472 0.2767528)  
        18) duration.in.month.binned=(30, Inf] 47   9.857143 bad (0.8503254 0.1496746) *
        19) duration.in.month.binned=(15,30] 101  33.000000 bad (0.6292135 0.3707865)  
          38) credit.amount.binned=(-Inf,5969.95],(9162.7, Inf] 92  29.142860 bad (0.6577181 0.3422819) *
          39) credit.amount.binned=(5969.95,9162.7] 9   0.000000 good (0.0000000 1.0000000) *
     5) duration.in.month.binned=(-Inf,6],(6,15] 306 100.285700 bad (0.6261981 0.3738019)  
      10) purpose.binned=education + car (new) + misc. level pos. + business + furniture/equipment,misc. level neg. 198  61.285710 bad (0.6767956 0.3232044) *
      11) purpose.binned=radio/television + car (used) 108  39.000000 bad (0.5042373 0.4957627)  
        22) duration.in.month.binned=(6,15] 89  30.857140 bad (0.5624578 0.4375422) *
        23) duration.in.month.binned=(-Inf,6] 19   0.000000 good (0.0000000 1.0000000) *
   3) savings.account.and.bonds.binned=misc. level neg.,unknown/ no savings account + 500 <= ... < 1000 DM 294 105.000000 bad (0.5212766 0.4787234)  
     6) purpose.binned=education + car (new) + misc. level pos. + business + furniture/equipment 161  54.000000 bad (0.6019656 0.3980344)  
      12) duration.in.month.binned=(15,30],(30, Inf] 87  27.428570 bad (0.6617733 0.3382267)  
        24) savings.account.and.bonds.binned=unknown/ no savings account + 500 <= ... < 1000 DM 73  21.857140 bad (0.7013663 0.2986337) *
        25) savings.account.and.bonds.binned=misc. level neg. 14   2.333333 good (0.2951807 0.7048193) *
      13) duration.in.month.binned=(-Inf,6],(6,15] 74  26.571430 bad (0.5130890 0.4869110)  
        26) savings.account.and.bonds.binned=misc. level neg. 18   6.000000 bad (0.6086957 0.3913043) *
        27) savings.account.and.bonds.binned=unknown/ no savings account + 500 <= ... < 1000 DM 56  18.666670 good (0.4757282 0.5242718) *
     7) purpose.binned=radio/television + car (used),misc. level neg. 133  32.666670 good (0.3904382 0.6095618)  
      14) credit.amount.binned=(5969.95,9162.7],(9162.7, Inf] 27   9.000000 bad (0.6086957 0.3913043) *
      15) credit.amount.binned=(-Inf,5969.95] 106  18.666670 good (0.3076923 0.6923077) *

rpart.plot(tree)

prunetree=prune(tree,cp=0.01)


rpart.plot(prunetree)

10.2c Loss Matrix

Change the code provided in the video such that a decision tree is constructed using a loss matrix penalizing 10 times more heavily for misclassified defaults.

set.seed(123)
treelossm <- rpart(creditability ~ ., data=na.omit(df_factors),control = rpart.control(cp = 0.001),model = TRUE,parms = list(loss = matrix(c(0, 3, 1, 0), ncol = 2)))
tree

n= 1000 

node), split, n, loss, yval, (yprob)
      * denotes terminal node

1) root 1000 300 good (0.3000000 0.7000000)  
  2) duration.in.month.binned=(30, Inf] 173  83 good (0.4797688 0.5202312) *
  3) duration.in.month.binned=(-Inf,6],(6,15],(15,30] 827 217 good (0.2623942 0.7376058)  
    6) credit.amount.binned=(9162.7, Inf] 11   6 bad (0.8181818 0.1818182) *
    7) credit.amount.binned=(-Inf,5969.95],(5969.95,9162.7] 816 208 good (0.2549020 0.7450980) *

rpart.plot(treelossm)

10.3 Pruning The Tree

10.3.1 Using Printcp and Plotcp

Using printcp we get the complexity parameter ( We can use this to set the tree depth, ie pruning ), we set the complexity parameter to get to the complexity level

Need to minimise xerror in the second last column of the table below

Need to set a seed to get the results reproducable

# We use printcp on the above tree 
printcp(tree)

Classification tree:
rpart(formula = creditability ~ ., data = na.omit(df_factors), 
    model = TRUE, parms = list(prior = c(0.7, 0.3)), control = rpart.control(cp = 1e-06))

Variables actually used in tree construction:
[1] credit.amount.binned             duration.in.month.binned         purpose.binned                   savings.account.and.bonds.binned

Root node error: 300/1000 = 0.3

n= 1000 

         CP nsplit rel error  xerror     xstd
1 0.0305556      0   1.00000 1.00000 0.020702
2 0.0166667      2   0.93889 1.01397 0.038818
3 0.0090476      3   0.92222 0.97571 0.036479
4 0.0057143      6   0.89508 0.99508 0.041416
5 0.0042857      9   0.87794 0.98365 0.041492
6 0.0000010     12   0.86508 0.94111 0.039309

plotcp(tree)

NA
NA

10.3.2 Using prune() to Prune the tree

ptree_pruned = prune(tree,cp=0.0002
                     )


rpart.plot(ptree_pruned,uniform=TRUE)

# prp function

prp(ptree_pruned,extra=1)

10.3.3 Other Tree Options

Weights: can use this to counter imbalance in tree

Other arguments are the rpart.control arguments

minsplit: min number in a node to split: default is 20 ..may be useful to lower it when data is unbalanced minbuckets: minimum in a leaf node (default is 1/3 of minsplit)

Evaluation : use the predict function

confusionMatrix(pred_tree_class,df_factors$creditability)

Confusion Matrix and Statistics

          Reference
Prediction bad good
      bad  283  513
      good  17  187
                                          
               Accuracy : 0.47            
                 95% CI : (0.4387, 0.5015)
    No Information Rate : 0.7             
    P-Value [Acc > NIR] : 1               
                                          
                  Kappa : 0.1429          
                                          
 Mcnemar's Test P-Value : <2e-16          
                                          
            Sensitivity : 0.9433          
            Specificity : 0.2671          
         Pos Pred Value : 0.3555          
         Neg Pred Value : 0.9167          
             Prevalence : 0.3000          
         Detection Rate : 0.2830          
   Detection Prevalence : 0.7960          
      Balanced Accuracy : 0.6052          
                                          
       'Positive' Class : bad             
                                          

10.4 Evaluating The Credit Risk Model

10.4.1 Strategy Curve

We dont want to exceed a certain % of defaults

10.4.2 ROC Curve

What is the best model overall ?

Accuracy is Maximized when a very high cutoff is chosen; when all testset items are classified as non default

ROC SHOWS CROSSOVER .. Harder to tell which is better

auc(ROC_tree)

Area under the curve: 0.7114

Which Variables do many AUC and compare

AUC Based Pruning..

# Example 1

# # Build four models each time deleting one variable in log_3_remove_ir
# log_4_remove_amnt <- glm(loan_status ~ grade + annual_inc + emp_cat, 
#                          family = binomial, data = training_set) 
# log_4_remove_grade <- glm(loan_status ~ loan_amnt  + annual_inc + emp_cat, 
#                           family = binomial, data = training_set)
# log_4_remove_inc <- glm(loan_status ~ loan_amnt + grade  + emp_cat , 
#                         family = binomial, data = training_set)
# log_4_remove_emp <- glm(loan_status ~ loan_amnt + grade + annual_inc, 
#                         family = binomial, data = training_set)
# 
# # Make PD-predictions for each of the models
# pred_4_remove_amnt <- predict(log_4_remove_amnt, newdata = test_set, type = "response")
# pred_4_remove_grade <- predict(log_4_remove_grade, newdata = test_set, type = "response")
# pred_4_remove_inc <- predict(log_4_remove_inc, newdata = test_set, type = "response")
# pred_4_remove_emp <- predict(log_4_remove_emp, newdata = test_set, type = "response")
# 
# # Compute the AUCs
# auc(test_set$loan_status, pred_4_remove_amnt)
# auc(test_set$loan_status, pred_4_remove_grade)
# auc(test_set$loan_status, pred_4_remove_inc)
# auc(test_set$loan_status, pred_4_remove_emp)


# Example 2

# # Build four models each time deleting one variable in log_3_remove_ir
# log_4_remove_amnt <- glm(loan_status ~ grade + annual_inc + emp_cat, 
#                          family = binomial, data = training_set) 
# log_4_remove_grade <- glm(loan_status ~ loan_amnt  + annual_inc + emp_cat, 
#                           family = binomial, data = training_set)
# log_4_remove_inc <- glm(loan_status ~ loan_amnt + grade  + emp_cat , 
#                         family = binomial, data = training_set)
# log_4_remove_emp <- glm(loan_status ~ loan_amnt + grade + annual_inc, 
#                         family = binomial, data = training_set)
# 
# # Make PD-predictions for each of the models
# pred_4_remove_amnt <- predict(log_4_remove_amnt, newdata = test_set, type = "response")
# pred_4_remove_grade <- predict(log_4_remove_grade, newdata = test_set, type = "response")
# pred_4_remove_inc <- predict(log_4_remove_inc, newdata = test_set, type = "response")
# pred_4_remove_emp <- predict(log_4_remove_emp, newdata = test_set, type = "response")
# 
# # Compute the AUCs
# auc(test_set$loan_status, pred_4_remove_amnt)
# auc(test_set$loan_status, pred_4_remove_grade)
# auc(test_set$loan_status, pred_4_remove_inc)
# auc(test_set$loan_status, pred_4_remove_emp)

11.0 Model Explanation