Customer Churn Prediction Exercise

Dataset: Credit card customers

Load packages

library(tidyverse)
library(skimr)
library(Hmisc)
library(ggsci)
library(caret)
library(ggpubr)
library(gridExtra)
library(corrplot)
library(rpart)
library(rpart.plot)
library(rattle)
library(randomForest)
library(caret)
library(xgboost)
library(pscl)
library(pROC)

Import data

df = read.csv("BankChurners.csv", header=TRUE)
head(df)

#delete naive bayes columns
df = df[,2:21]

dim(df)

[1] 10127    20

colnames(df)

 [1] "Attrition_Flag"           "Customer_Age"             "Gender"                   "Dependent_count"         
 [5] "Education_Level"          "Marital_Status"           "Income_Category"          "Card_Category"           
 [9] "Months_on_book"           "Total_Relationship_Count" "Months_Inactive_12_mon"   "Contacts_Count_12_mon"   
[13] "Credit_Limit"             "Total_Revolving_Bal"      "Avg_Open_To_Buy"          "Total_Amt_Chng_Q4_Q1"    
[17] "Total_Trans_Amt"          "Total_Trans_Ct"           "Total_Ct_Chng_Q4_Q1"      "Avg_Utilization_Ratio"   

Target variable

Hmisc::describe(df$Attrition_Flag)

df$Attrition_Flag 
       n  missing distinct 
   10127        0        2 
                                              
Value      Attrited Customer Existing Customer
Frequency               1627              8500
Proportion             0.161             0.839

#convert target variable to numeric 
df$label= ifelse(df$Attrition_Flag =="Attrited Customer","1","0")
df$label = as.factor(df$label)
Hmisc::describe(df$label)

df$label 
       n  missing distinct 
   10127        0        2 
                      
Value          0     1
Frequency   8500  1627
Proportion 0.839 0.161

#drop attrition flag 
df = subset(df, select = -c(Attrition_Flag))

• For the purpose of this exercise, Attrited Customer and Existing Customer classes are releveled to 1 and 0 respectively.
• Out of 10127 obs, 1627 (16.1%) are attrited customers suggesting an imbalanced dataset.

#summary
df = df %>% mutate_at(vars(label, Gender,Education_Level,Income_Category,Marital_Status,Card_Category),list(factor))
skim(df)

── Data Summary ────────────────────────
                           Values
Name                       df    
Number of rows             10127 
Number of columns          20    
_______________________          
Column type frequency:           
  factor                   6     
  numeric                  14    
________________________         
Group variables            None  

── Variable type: factor ──────────────────────────────────────────────────────────────────────────────────────────
  skim_variable   n_missing complete_rate ordered n_unique top_counts                                
1 Gender                  0             1 FALSE          2 F: 5358, M: 4769                          
2 Education_Level         0             1 FALSE          7 Gra: 3128, Hig: 2013, Unk: 1519, Une: 1487
3 Marital_Status          0             1 FALSE          4 Mar: 4687, Sin: 3943, Unk: 749, Div: 748  
4 Income_Category         0             1 FALSE          6 Les: 3561, $40: 1790, $80: 1535, $60: 1402
5 Card_Category           0             1 FALSE          4 Blu: 9436, Sil: 555, Gol: 116, Pla: 20    
6 label                   0             1 FALSE          2 0: 8500, 1: 1627                          

── Variable type: numeric ─────────────────────────────────────────────────────────────────────────────────────────
   skim_variable            n_missing complete_rate     mean       sd    p0      p25      p50       p75      p100
 1 Customer_Age                     0             1   46.3      8.02    26    41       46        52        73    
 2 Dependent_count                  0             1    2.35     1.30     0     1        2         3         5    
 3 Months_on_book                   0             1   35.9      7.99    13    31       36        40        56    
 4 Total_Relationship_Count         0             1    3.81     1.55     1     3        4         5         6    
 5 Months_Inactive_12_mon           0             1    2.34     1.01     0     2        2         3         6    
 6 Contacts_Count_12_mon            0             1    2.46     1.11     0     2        2         3         6    
 7 Credit_Limit                     0             1 8632.    9089.    1438. 2555     4549     11068.    34516    
 8 Total_Revolving_Bal              0             1 1163.     815.       0   359     1276      1784      2517    
 9 Avg_Open_To_Buy                  0             1 7469.    9091.       3  1324.    3474      9859     34516    
10 Total_Amt_Chng_Q4_Q1             0             1    0.760    0.219    0     0.631    0.736     0.859     3.40 
11 Total_Trans_Amt                  0             1 4404.    3397.     510  2156.    3899      4741     18484    
12 Total_Trans_Ct                   0             1   64.9     23.5     10    45       67        81       139    
13 Total_Ct_Chng_Q4_Q1              0             1    0.712    0.238    0     0.582    0.702     0.818     3.71 
14 Avg_Utilization_Ratio            0             1    0.275    0.276    0     0.023    0.176     0.503     0.999
   hist 
 1 ▂▆▇▃▁
 2 ▇▇▇▅▁
 3 ▁▃▇▃▂
 4 ▇▇▆▆▆
 5 ▅▇▇▁▁
 6 ▅▇▇▃▁
 7 ▇▂▁▁▁
 8 ▇▅▇▇▅
 9 ▇▂▁▁▁
10 ▅▇▁▁▁
11 ▇▅▁▁▁
12 ▂▅▇▂▁
13 ▇▆▁▁▁
14 ▇▂▂▂▁

Exploratory data analysis

p1 = df %>% group_by(label,Contacts_Count_12_mon) %>% tally() %>% mutate(prop=n/sum(n)) %>% ggplot(aes(x=Contacts_Count_12_mon, y=prop,fill=label)) + geom_col(position="dodge") + scale_fill_jama() + labs(y="proportion") + theme_light() + theme(legend.position="bottom")
p2 = df %>% group_by(label,Months_Inactive_12_mon) %>% tally() %>% mutate(prop=n/sum(n)) %>% ggplot(aes(x=Months_Inactive_12_mon, y=prop,fill=label)) + geom_col(position="dodge") + scale_fill_jama() + labs(y="proportion") + theme_light() + theme(legend.position="bottom")
grid.arrange(p1,p2,ncol=2,nrow=1)

• Contacts_Count_12_mon: customers with 3 or more contacts count in the past 12 months have a higher proportion of attrition
• Months_Inactive_12_mon: customers with 3 or more inactive months in the past 12 months have a higher proportion of attrition

p3 = df %>% group_by(label,Total_Relationship_Count) %>% tally() %>% mutate(prop=n/sum(n)) %>% ggplot(aes(x=Total_Relationship_Count, y=prop,fill=label)) + geom_col(position="dodge") + scale_fill_jama() + labs(y="proportion") + theme_light() + theme(legend.position="bottom") 
p4 = df %>% group_by(label,Dependent_count) %>% tally() %>% mutate(prop=n/sum(n)) %>% ggplot(aes(x=Dependent_count, y=prop,fill=label)) + geom_col(position="dodge") + scale_fill_jama() + labs(y="proportion") + theme_light() + theme(legend.position="bottom")
grid.arrange(p3,p4,ncol=2,nrow=1)

• Total_Relationship_Count: customers with 3 or less total relationship count has a higher proportion of attrition
• Dependent_count: customers with 3 or more dependents have a higher proportion of attrition

p5 = df %>% group_by(Gender,label) %>% tally() %>% mutate(prop=n/sum(n)) %>% ggplot(aes(x=label, y=prop,fill=Gender)) + geom_col(position="dodge") + scale_fill_jama() + labs(y="proportion") + theme_light() + theme(legend.position="bottom")
p6 = df %>% group_by(label,Education_Level) %>% tally() %>% mutate(prop=n/sum(n)) %>% ggplot(aes(x=Education_Level, y=prop,fill=label)) + geom_col(position="dodge") + scale_fill_jama() + labs(y="proportion") + theme_light() + theme(legend.position="bottom") + coord_flip()
grid.arrange(p5,p6,ncol=2,nrow=1)

• Gender: Female has a higher proportion of attrition compared to males.
• Education_Level: The attrition class has a higher proportion of doctorate, post-graduate and unknown education level compared to the existing customer class.

Density plots (inspired by Who’s gonna churn? by Carmine Minichini)

p7 = df %>% ggplot(aes(x=Total_Trans_Ct, fill=label)) + geom_density(alpha=0.6) + scale_fill_jama() + theme_light() 
p8 = df %>% ggplot(aes(x=Total_Ct_Chng_Q4_Q1, fill=label)) + geom_density(alpha=0.6) + scale_fill_jama() + theme_light()
p9 = df %>% ggplot(aes(x=Total_Revolving_Bal, fill=label)) + geom_density(alpha=0.6) + scale_fill_jama() + theme_light()
p10 = df %>% ggplot(aes(x=Avg_Utilization_Ratio, fill=label)) + geom_density(alpha=0.6) + scale_fill_jama() + theme_light()
p11 = df %>% ggplot(aes(x=Total_Trans_Amt, fill=label)) + geom_density(alpha=0.6) + scale_fill_jama() + theme_light()
p12 = df %>% ggplot(aes(x=Credit_Limit, fill=label)) + geom_density(alpha=0.6) + scale_fill_jama() + theme_light()
grid.arrange(p7,p8,ncol=1,nrow=2)

grid.arrange(p9,p10,ncol=1,nrow=2)

grid.arrange(p11,p12,ncol=1,nrow=2)

• The attrition class has a lower total transaction count, total count change, total revolving balance, average utilization ratio and total transactions amount compared to the existing customers class, as expected.

Feature selection

#check correlation of all numeric variables
df_num = select_if(df,is.numeric)
df_num = data.frame(lapply(df_num, function(x) as.numeric(as.character(x))))
res=cor(df_num)
corrplot(res, type="upper", tl.col="#636363",tl.cex=0.5 )

• Avg_Open_To_Buy is highly correlated to Credit_Limit
• Total_Trans amount is highly correlated to Total_Trans_Amt
• Total_Amt_Chng_Q4_Q1 is correlated to Total_Ct_Cng_Q4_Q1
• Avg_Utilization_Ratio is correlated to Avg_Open_To_BUy, Total_Revolving_Bal and Credit_Limit

#drop Months_on_book,Total_Trans_Amt, Total_Amt_Chng_Q4_Q1, Avg_Utilization_Ratio
df1 = df %>% select(-c(Months_on_book,Total_Trans_Amt, Total_Amt_Chng_Q4_Q1, Avg_Utilization_Ratio, Avg_Open_To_Buy))
dim(df1)

#check correlation after dropping variables
df1_num = select_if(df1,is.numeric)
df1_num = data.frame(lapply(df1_num, function(x) as.numeric(as.character(x))))
res2=cor(df1_num)
corrplot(res2, type="lower", tl.col="#636363",tl.cex=0.5 )

Test and train set

trainIndex <- createDataPartition(df1$label, p = .75,list=FALSE)
training <- df1[trainIndex,]
testing <- df1[-trainIndex,]

Hmisc::describe(training$label)

training$label 
       n  missing distinct 
    7596        0        2 
                      
Value          0     1
Frequency   6375  1221
Proportion 0.839 0.161

Hmisc::describe(testing$label)

testing$label 
       n  missing distinct 
    2531        0        2 
                    
Value         0    1
Frequency  2125  406
Proportion 0.84 0.16

Logistic regression

model1= glm(label ~., data=training, family = "binomial")
summary(model1) 

Call:
glm(formula = label ~ ., family = "binomial", data = training)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8443  -0.4139  -0.2076  -0.0890   3.5346  

Coefficients:
                                Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    5.250e+00  4.985e-01  10.533  < 2e-16 ***
Customer_Age                  -9.093e-03  5.134e-03  -1.771 0.076564 .  
GenderM                       -7.165e-01  1.600e-01  -4.478 7.55e-06 ***
Dependent_count                1.162e-01  3.308e-02   3.512 0.000444 ***
Education_LevelDoctorate       2.583e-01  2.324e-01   1.111 0.266391    
Education_LevelGraduate       -6.983e-02  1.539e-01  -0.454 0.649965    
Education_LevelHigh School     7.434e-02  1.630e-01   0.456 0.648452    
Education_LevelPost-Graduate   4.239e-01  2.226e-01   1.905 0.056843 .  
Education_LevelUneducated      1.176e-01  1.724e-01   0.682 0.495204    
Education_LevelUnknown         9.088e-02  1.708e-01   0.532 0.594734    
Marital_StatusMarried         -3.225e-01  1.682e-01  -1.917 0.055228 .  
Marital_StatusSingle           1.483e-01  1.690e-01   0.878 0.380112    
Marital_StatusUnknown          9.930e-02  2.172e-01   0.457 0.647562    
Income_Category$40K - $60K    -4.664e-01  2.304e-01  -2.024 0.042935 *  
Income_Category$60K - $80K    -3.511e-01  2.073e-01  -1.693 0.090437 .  
Income_Category$80K - $120K    1.470e-01  1.903e-01   0.773 0.439754    
Income_CategoryLess than $40K -3.622e-01  2.493e-01  -1.453 0.146253    
Income_CategoryUnknown        -4.259e-01  2.597e-01  -1.640 0.101059    
Card_CategoryGold              1.389e+00  4.092e-01   3.394 0.000690 ***
Card_CategoryPlatinum          1.483e+00  8.517e-01   1.741 0.081654 .  
Card_CategorySilver            4.279e-01  2.312e-01   1.851 0.064126 .  
Total_Relationship_Count      -5.489e-01  3.058e-02 -17.947  < 2e-16 ***
Months_Inactive_12_mon         4.717e-01  4.106e-02  11.488  < 2e-16 ***
Contacts_Count_12_mon          4.393e-01  3.942e-02  11.145  < 2e-16 ***
Credit_Limit                  -9.382e-06  6.983e-06  -1.344 0.179063    
Total_Revolving_Bal           -9.144e-04  5.151e-05 -17.751  < 2e-16 ***
Total_Trans_Ct                -6.611e-02  2.517e-03 -26.264  < 2e-16 ***
Total_Ct_Chng_Q4_Q1           -2.767e+00  2.004e-01 -13.809  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 6698.1  on 7595  degrees of freedom
Residual deviance: 3871.1  on 7568  degrees of freedom
AIC: 3927.1

Number of Fisher Scoring iterations: 6

pR2(model1)  

fitting null model for pseudo-r2
          llh       llhNull            G2      McFadden          r2ML          r2CU 
-1935.5629265 -3349.0692506  2827.0126481     0.4220594     0.3107638     0.5303478 

anova(model1, test= "Chisq")

Analysis of Deviance Table

Model: binomial, link: logit

Response: label

Terms added sequentially (first to last)

                         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                                      7595     6698.1              
Customer_Age              1     1.71      7594     6696.4 0.1911947    
Gender                    1    12.06      7593     6684.4 0.0005162 ***
Dependent_count           1     1.76      7592     6682.6 0.1849018    
Education_Level           6     8.61      7586     6674.0 0.1969515    
Marital_Status            3     3.54      7583     6670.5 0.3154854    
Income_Category           5     7.35      7578     6663.1 0.1961685    
Card_Category             3     2.66      7575     6660.5 0.4477147    
Total_Relationship_Count  1   192.81      7574     6467.7 < 2.2e-16 ***
Months_Inactive_12_mon    1   170.79      7573     6296.9 < 2.2e-16 ***
Contacts_Count_12_mon     1   362.38      7572     5934.5 < 2.2e-16 ***
Credit_Limit              1     6.66      7571     5927.8 0.0098542 ** 
Total_Revolving_Bal       1   457.22      7570     5470.6 < 2.2e-16 ***
Total_Trans_Ct            1  1345.97      7569     4124.6 < 2.2e-16 ***
Total_Ct_Chng_Q4_Q1       1   253.50      7568     3871.1 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

• 8 variables (Gender, Total_Relationship_Count, Months_Inactive_12_mon, Contacts_Count_12_mon, Credit_Limit, Total_Revolving_Bal, Total_Trans_Ct and Total_Ct_Chng_Q4_Q1) are significant variables in predicting customer churn for credit card services.
  • Analysis of deviance table shows that Total_Trans_Ct has most significant variable.
• Logistic regression model suggests that:
  • The log odds of a male customer churning is 0.7 lower than female.
  • For one unit increase in dependent count, the log odds of the customer churning increases by 0.1.
  • Versus income category 120K+, the 40K-60K category decreases the log odds of the customer churning by 0.5.
  • Versus the blue card category, gold category increases the log odds of the customer churning by 1.4
  • For one unit increase in total relationship count, the log odds of the customer churning decreases by 0.5.
  • For one unit increase in months inactive in the last 12 months, the log odds of the customer churning increases by 0.5.
  • For one unit increase in total revolving balance, the log odds of the customer churning decreases by 0.0009.
  • For one unit increase in total transaction count, the log odds of the customer churning decreases by 0.07.
  • For one unit increase Change in Transaction Count (Q4 over Q1) , the log odds of the customer churning decreases by 2.8.

prob=predict(model1,testing,type="response")
prob1=rep(0,2531)
prob1[prob>0.2]=1
cmlr = confusionMatrix(as.factor(prob1), testing$label, positive="1")
cmlr

Confusion Matrix and Statistics

          Reference
Prediction    0    1
         0 1870   87
         1  255  319
                                         
               Accuracy : 0.8649         
                 95% CI : (0.8509, 0.878)
    No Information Rate : 0.8396         
    P-Value [Acc > NIR] : 0.0002231      
                                         
                  Kappa : 0.5703         
                                         
 Mcnemar's Test P-Value : < 2.2e-16      
                                         
            Sensitivity : 0.7857         
            Specificity : 0.8800         
         Pos Pred Value : 0.5557         
         Neg Pred Value : 0.9555         
             Prevalence : 0.1604         
         Detection Rate : 0.1260         
   Detection Prevalence : 0.2268         
      Balanced Accuracy : 0.8329         
                                         
       'Positive' Class : 1              
                                         

round(cmlr$byClass["F1"], 4)

   F1 
0.651 

roc_lr2 = roc(testing$label, prob1, plot=TRUE, print.auc=TRUE)

Setting levels: control = 0, case = 1
Setting direction: controls < cases

• 319 out of 406 positive instances were predicted correctly (recall of 0.786) with logistic regression (probability value 0.2)

Decision tree

mt = rpart(label ~., data = training, method = "class")
plotcp(mt)

mt_prune = prune(mt,cp=0.036)
fancyRpartPlot(mt_prune)

printcp(mt_prune)

Classification tree:
rpart(formula = label ~ ., data = training, method = "class")

Variables actually used in tree construction:
[1] Total_Relationship_Count Total_Revolving_Bal      Total_Trans_Ct          

Root node error: 1221/7596 = 0.16074

n= 7596 

        CP nsplit rel error  xerror     xstd
1 0.165029      0   1.00000 1.00000 0.026217
2 0.074529      2   0.66994 0.68223 0.022304
3 0.036000      3   0.59541 0.60852 0.021204

mt_prune$variable.importance

          Total_Trans_Ct      Total_Revolving_Bal Total_Relationship_Count      Total_Ct_Chng_Q4_Q1 
              356.304278               335.977874               141.636180               108.826027 
            Credit_Limit    Contacts_Count_12_mon   Months_Inactive_12_mon             Customer_Age 
               53.946861                 8.022617                 3.047734                 2.168455 

tree.p = predict(mt_prune, testing, type = "class")
cmt = confusionMatrix(tree.p, testing$label, positive ="1")
cmt

Confusion Matrix and Statistics

          Reference
Prediction    0    1
         0 2039  157
         1   86  249
                                          
               Accuracy : 0.904           
                 95% CI : (0.8918, 0.9152)
    No Information Rate : 0.8396          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.6164          
                                          
 Mcnemar's Test P-Value : 7.106e-06       
                                          
            Sensitivity : 0.61330         
            Specificity : 0.95953         
         Pos Pred Value : 0.74328         
         Neg Pred Value : 0.92851         
             Prevalence : 0.16041         
         Detection Rate : 0.09838         
   Detection Prevalence : 0.13236         
      Balanced Accuracy : 0.78641         
                                          
       'Positive' Class : 1               
                                          

round(cmt$byClass["F1"], 4)

    F1 
0.6721 

testing$tp1= tree.p
roc_t= roc(response= testing$label, predictor = factor(testing$tp1, ordered=TRUE), plot=TRUE, print.auc=TRUE)

Setting levels: control = 0, case = 1
Setting direction: controls < cases

*249 out of 486 positive instances were predicted correctly (recall of 0.613) with classification tree.

Random forest

trControl <- trainControl(method = "cv",
    number = 10,
    search = "grid")

set.seed(1234)
rf1 = train(label ~ .,data = training,method="rf",metric ="Accuracy",trControl = trControl)
print(rf1)

plot(rf1)

varImp(rf1)

rf variable importance

  only 20 most important variables shown (out of 27)

rfpred = predict(rf1, testing)
cmrf = confusionMatrix(rfpred, testing$label,positive="1")
cmrf

Confusion Matrix and Statistics

          Reference
Prediction    0    1
         0 2076  129
         1   49  277
                                         
               Accuracy : 0.9297         
                 95% CI : (0.919, 0.9393)
    No Information Rate : 0.8396         
    P-Value [Acc > NIR] : < 2.2e-16      
                                         
                  Kappa : 0.7163         
                                         
 Mcnemar's Test P-Value : 3.194e-09      
                                         
            Sensitivity : 0.6823         
            Specificity : 0.9769         
         Pos Pred Value : 0.8497         
         Neg Pred Value : 0.9415         
             Prevalence : 0.1604         
         Detection Rate : 0.1094         
   Detection Prevalence : 0.1288         
      Balanced Accuracy : 0.8296         
                                         
       'Positive' Class : 1              
                                         

round(cmrf$byClass["F1"], 4)

    F1 
0.7568 

testing$rfp= rfpred
roc_rf= roc(response= testing$label, predictor = factor(testing$rfp, ordered=TRUE), plot=TRUE, print.auc=TRUE)

Setting levels: control = 0, case = 1
Setting direction: controls < cases

• 277 out of 486 positive instances were predicted correctly (recall of 0.682) with random forest.

XGBoost

XGBoost code reference: Who’s gonna churn? by Carmine Minichini

target_column = df1$label
data =  df1 %>% select(-label)
dmy = dummyVars(" ~ .", data = data)
train_data = data.frame(predict(dmy, newdata = data))
data <- cbind(train_data,target_column)
names(data)[33] <- 'label'

trainIndex <- createDataPartition(data$label,p=0.75,list=FALSE)
data_train <- data[trainIndex,]
data_test <-  data[-trainIndex,]

grid_train = data_train
levels(grid_train$label) <- c("X0","X1")
#grid parameters
xgb_grid_1 = expand.grid(
    nrounds = 10,
    eta = seq(2,10,by=1)/10,
    max_depth = c(6, 8, 10),
    gamma = 0,
    subsample = c(0.5, 0.75, 1),
    min_child_weight = c(1,2) ,
    colsample_bytree = c(0.3,0.5)
  )
# pack the training control parameters
  xgb_trcontrol_1 = trainControl(
    method = "cv",
    number = 2,
    search='grid',
    verboseIter = FALSE,
    returnData = TRUE,
    returnResamp = "all", # save losses across all models
    classProbs = TRUE, # set to TRUE for AUC to be computed
    summaryFunction = prSummary, # probability summary(AUC)
    allowParallel = TRUE,
  )

xgb_train_1 = train(
    x = as.matrix(grid_train %>% select(-label)),
    y = factor(grid_train$label),
    trControl = xgb_trcontrol_1,
    tuneGrid = xgb_grid_1,
    method = "xgbTree",
    metric= 'Recall'
  )

best_tune <- xgb_train_1$bestTune
results <- xgb_train_1$results
trained_model <- xgb_train_1
cat(paste("",
            paste('With a recall of:',results[rownames(best_tune),"Recall"]),
            'Best GRIDSEARCH Hyperparameters:',
            '',
            sep='\n\n'))

With a recall of: 0.997489999169304

Best GRIDSEARCH Hyperparameters:

  
rownames(best_tune) <- 'Value'
print(t(best_tune))

                 Value
nrounds          10.00
max_depth         6.00
eta               0.20
gamma             0.00
colsample_bytree  0.30
min_child_weight  1.00
subsample         0.75

#out dataframe
gridresults = results
best_tune = best_tune
train_data = data_train
test_data = data_test

#best hyperparameters from gridsearch
best_tune <- best_tune
#train
data_train <- train_data %>% select(-label)
label_train <- train_data$label
#test
data_test <- test_data %>% select(-label)
label_test <- test_data$label
# as matrix
data_train <- as.matrix(data_train)
data_test <- as.matrix(data_test)
# as numeric
label_train <- as.numeric(label_train)
label_test <- as.numeric(label_test)
#relevel
label_train= ifelse(label_train>1,1,0)
label_test= ifelse(label_test>1,1,0)

#XGB matrix
dtrain <-  xgb.DMatrix(data_train,label=label_train)
dtest <- xgb.DMatrix(data_test,label=label_test)

#XGB model
 model <- xgboost(data= dtrain, 
                   objective = "binary:logistic",
                   max_depth = best_tune$max_depth,
                   nrounds=100,
                   colsample_bytree = best_tune$colsample_bytree,
                   gamma = best_tune$gamma,
                   min_child_weight = best_tune$min_child_weight,
                   eta = best_tune$eta, 
                   subsample = best_tune$subsample,
                   print_every_n = 20,
                   scale_pos_weight=5.22,
                   max_delta_step=1,
                   eval_metric='aucpr',
                   verbose=1,
                   nthread = 4)

[1] train-aucpr:0.452852 
[21]    train-aucpr:0.905627 
[41]    train-aucpr:0.946657 
[61]    train-aucpr:0.963228 
[81]    train-aucpr:0.974350 
[100]   train-aucpr:0.982677 

cv  <-  xgb.cv(data = dtrain, 
                 nround = 50, 
                 print_every_n= 10,
                 verbose = TRUE,
                 metrics = list("aucpr"),
                 nfold = 5, 
                 nthread = 4,
                 objective = "binary:logistic",
                 prediction=F)

[1] train-aucpr:0.820787+0.011257   test-aucpr:0.754658+0.031137 
[11]    train-aucpr:0.938609+0.003525   test-aucpr:0.848662+0.024433 
[21]    train-aucpr:0.963758+0.004984   test-aucpr:0.862755+0.019337 
[31]    train-aucpr:0.978737+0.003634   test-aucpr:0.867039+0.020116 
[41]    train-aucpr:0.988214+0.001604   test-aucpr:0.866213+0.021019 
[50]    train-aucpr:0.993660+0.001419   test-aucpr:0.867678+0.021359 

out <- list(data_train = data_train,
              dtest = dtest,
              label_test = label_test,
              model = model)

data_train <- data_train
pred <- predict(model,dtest)
prediction <- as.numeric(pred > 0.5)
cm <- confusionMatrix(factor(prediction),factor(label_test),positive="1")
cm

Confusion Matrix and Statistics

          Reference
Prediction    0    1
         0 1989   72
         1  136  334
                                          
               Accuracy : 0.9178          
                 95% CI : (0.9064, 0.9282)
    No Information Rate : 0.8396          
    P-Value [Acc > NIR] : < 2.2e-16       
                                          
                  Kappa : 0.7132          
                                          
 Mcnemar's Test P-Value : 1.252e-05       
                                          
            Sensitivity : 0.8227          
            Specificity : 0.9360          
         Pos Pred Value : 0.7106          
         Neg Pred Value : 0.9651          
             Prevalence : 0.1604          
         Detection Rate : 0.1320          
   Detection Prevalence : 0.1857          
      Balanced Accuracy : 0.8793          
                                          
       'Positive' Class : 1               
                                          

round(cm$byClass["F1"], 4)

    F1 
0.7626 

roc.curve = roc(response = label_test,
                  predictor = prediction,
                  levels=c(0, 1),quiet = T) 
plot(roc.curve,print.auc=TRUE)

*334 out of 486 positive instances are predicted correctly with XGB (0.823)

importance_matrix <- xgb.importance(colnames(data_train), model = model)
  xgb.plot.importance(importance_matrix,
                      top_n=10,
                      main='Features Importance',
                      measure = 'Frequency')

Summary of exercise

	Recall	AUC	F1 Score
Logistic regression	0.786	0.833	0.651
Decision tree	0.613	0.786	0.672
Random Forest	0.682	0.83	0.757
XGBoost	0.823	0.879	0.763