Default Credit Cards

1 Introducción

La base de datos consiste en clientes de una institución bancaria donde el objetivo es entender el default de pago en tarjetas de crédito en el mercado de Taiwan con el fin de dar claridad acerca de la base de clientes al equipo de riesgo de crédito.

La variable dependiente (target) es ‘default payment next month’ consiste en una variable binaria definida como

\[ y = \begin{cases} 1 & \text{El cliente ha tenido default en pago TC} \\ 0 & \text{El cliente no ha tenido default en pago TC} \end{cases} \]

Source: https://archive.ics.uci.edu/dataset/350/default+of+credit+card+clients

2 Descripción de Variables

Las variables que provienen en el dataset se componen en la sig. estructura

Variable	Descripción
ID	Identificador único por cliente
LIMIT_BAL	Balance aprobado en TC
SEX	Género cliente (1 = masculino; 2 = femenino)
EDUCATION	Nivel Educativo (0=sin escolaridad; 1 =escuela; 2 = universidad; 3 = bachillerato)
MARRIAGE	Estado Civil (1 = casado; 2 = soltero; 3 = otros)
AGE	Edad
PAY_0,PAY_2 - PAY_6	Histórico de pagos desde Abril a Sep 2005. PAY_6: pago mes Sep, PAY_5: pago mes Ago, … , PAY_1: pago Abr. Los valores de status de pago: -1 = pago al día, 1 = retraso pago 1 mes, 2 = retraso pago 2 meses, … , 9= retraso de pago por ≥ 9 meses
BILL_AMT1 - BILL_AMT6	Monto total adeudado por mes desde Abril a Sep 2005
PAY_AMT1 - PAY_AMT6	Monto total pagado por mes desde Abril a Sep 2005

3 Tratamiento de los Datos

Se realizarán correcciones y ajustes a los datos que permitan generar el EDA (Exploratory Data Analysis) de las variables. El primer paso consistirá en leer los datos, clasificar las variables cuantitativas, cualitativas, target, limpieza de algunas variables, determinar si los ID de los clientes son únicos

library(tidyverse)
library(dplyr)
library(readxl)
library(ggplot2)
library(tidymodels)
library(skedastic)
library(tseries)
library(regclass)
library(lmtest)
library(vip)
library(ranger)
library(xgboost)

setwd("D:/Documents/Big Data Master/Visualización de Datos/Default Credit Card")

data <- read_excel("default of credit card clients.xls",
                                             sheet = "Data", skip = 1)

str(data)

tibble [30,000 × 25] (S3: tbl_df/tbl/data.frame)
 $ ID                        : num [1:30000] 1 2 3 4 5 6 7 8 9 10 ...
 $ LIMIT_BAL                 : num [1:30000] 20000 120000 90000 50000 50000 50000 500000 100000 140000 20000 ...
 $ SEX                       : num [1:30000] 2 2 2 2 1 1 1 2 2 1 ...
 $ EDUCATION                 : num [1:30000] 2 2 2 2 2 1 1 2 3 3 ...
 $ MARRIAGE                  : num [1:30000] 1 2 2 1 1 2 2 2 1 2 ...
 $ AGE                       : num [1:30000] 24 26 34 37 57 37 29 23 28 35 ...
 $ PAY_0                     : num [1:30000] 2 -1 0 0 -1 0 0 0 0 -2 ...
 $ PAY_2                     : num [1:30000] 2 2 0 0 0 0 0 -1 0 -2 ...
 $ PAY_3                     : num [1:30000] -1 0 0 0 -1 0 0 -1 2 -2 ...
 $ PAY_4                     : num [1:30000] -1 0 0 0 0 0 0 0 0 -2 ...
 $ PAY_5                     : num [1:30000] -2 0 0 0 0 0 0 0 0 -1 ...
 $ PAY_6                     : num [1:30000] -2 2 0 0 0 0 0 -1 0 -1 ...
 $ BILL_AMT1                 : num [1:30000] 3913 2682 29239 46990 8617 ...
 $ BILL_AMT2                 : num [1:30000] 3102 1725 14027 48233 5670 ...
 $ BILL_AMT3                 : num [1:30000] 689 2682 13559 49291 35835 ...
 $ BILL_AMT4                 : num [1:30000] 0 3272 14331 28314 20940 ...
 $ BILL_AMT5                 : num [1:30000] 0 3455 14948 28959 19146 ...
 $ BILL_AMT6                 : num [1:30000] 0 3261 15549 29547 19131 ...
 $ PAY_AMT1                  : num [1:30000] 0 0 1518 2000 2000 ...
 $ PAY_AMT2                  : num [1:30000] 689 1000 1500 2019 36681 ...
 $ PAY_AMT3                  : num [1:30000] 0 1000 1000 1200 10000 657 38000 0 432 0 ...
 $ PAY_AMT4                  : num [1:30000] 0 1000 1000 1100 9000 ...
 $ PAY_AMT5                  : num [1:30000] 0 0 1000 1069 689 ...
 $ PAY_AMT6                  : num [1:30000] 0 2000 5000 1000 679 ...
 $ default payment next month: num [1:30000] 1 1 0 0 0 0 0 0 0 0 ...

Las variables Sex, Marriage, default payment next month son categóricas nominales, Education categórica ordinal, se procederán a transformar en factor

# Esto garantiza cambiar el tipo de dato para múltiples variables al mismo tiempo

data<- data %>% 
        mutate_at(.vars = vars(SEX,EDUCATION,MARRIAGE,"default payment next month"),.fun = factor)

summary(data)

       ID          LIMIT_BAL       SEX       EDUCATION MARRIAGE 
 Min.   :    1   Min.   :  10000   1:11888   0:   14   0:   54  
 1st Qu.: 7501   1st Qu.:  50000   2:18112   1:10585   1:13659  
 Median :15000   Median : 140000             2:14030   2:15964  
 Mean   :15000   Mean   : 167484             3: 4917   3:  323  
 3rd Qu.:22500   3rd Qu.: 240000             4:  123            
 Max.   :30000   Max.   :1000000             5:  280            
                                             6:   51            
      AGE            PAY_0             PAY_2             PAY_3        
 Min.   :21.00   Min.   :-2.0000   Min.   :-2.0000   Min.   :-2.0000  
 1st Qu.:28.00   1st Qu.:-1.0000   1st Qu.:-1.0000   1st Qu.:-1.0000  
 Median :34.00   Median : 0.0000   Median : 0.0000   Median : 0.0000  
 Mean   :35.49   Mean   :-0.0167   Mean   :-0.1338   Mean   :-0.1662  
 3rd Qu.:41.00   3rd Qu.: 0.0000   3rd Qu.: 0.0000   3rd Qu.: 0.0000  
 Max.   :79.00   Max.   : 8.0000   Max.   : 8.0000   Max.   : 8.0000  
                                                                      
     PAY_4             PAY_5             PAY_6           BILL_AMT1      
 Min.   :-2.0000   Min.   :-2.0000   Min.   :-2.0000   Min.   :-165580  
 1st Qu.:-1.0000   1st Qu.:-1.0000   1st Qu.:-1.0000   1st Qu.:   3559  
 Median : 0.0000   Median : 0.0000   Median : 0.0000   Median :  22382  
 Mean   :-0.2207   Mean   :-0.2662   Mean   :-0.2911   Mean   :  51223  
 3rd Qu.: 0.0000   3rd Qu.: 0.0000   3rd Qu.: 0.0000   3rd Qu.:  67091  
 Max.   : 8.0000   Max.   : 8.0000   Max.   : 8.0000   Max.   : 964511  
                                                                        
   BILL_AMT2        BILL_AMT3         BILL_AMT4         BILL_AMT5     
 Min.   :-69777   Min.   :-157264   Min.   :-170000   Min.   :-81334  
 1st Qu.:  2985   1st Qu.:   2666   1st Qu.:   2327   1st Qu.:  1763  
 Median : 21200   Median :  20089   Median :  19052   Median : 18105  
 Mean   : 49179   Mean   :  47013   Mean   :  43263   Mean   : 40311  
 3rd Qu.: 64006   3rd Qu.:  60165   3rd Qu.:  54506   3rd Qu.: 50191  
 Max.   :983931   Max.   :1664089   Max.   : 891586   Max.   :927171  
                                                                      
   BILL_AMT6          PAY_AMT1         PAY_AMT2          PAY_AMT3     
 Min.   :-339603   Min.   :     0   Min.   :      0   Min.   :     0  
 1st Qu.:   1256   1st Qu.:  1000   1st Qu.:    833   1st Qu.:   390  
 Median :  17071   Median :  2100   Median :   2009   Median :  1800  
 Mean   :  38872   Mean   :  5664   Mean   :   5921   Mean   :  5226  
 3rd Qu.:  49198   3rd Qu.:  5006   3rd Qu.:   5000   3rd Qu.:  4505  
 Max.   : 961664   Max.   :873552   Max.   :1684259   Max.   :896040  
                                                                      
    PAY_AMT4         PAY_AMT5           PAY_AMT6       
 Min.   :     0   Min.   :     0.0   Min.   :     0.0  
 1st Qu.:   296   1st Qu.:   252.5   1st Qu.:   117.8  
 Median :  1500   Median :  1500.0   Median :  1500.0  
 Mean   :  4826   Mean   :  4799.4   Mean   :  5215.5  
 3rd Qu.:  4013   3rd Qu.:  4031.5   3rd Qu.:  4000.0  
 Max.   :621000   Max.   :426529.0   Max.   :528666.0  
                                                       
 default payment next month
 0:23364                   
 1: 6636                   
                           
                           
                           
                           
                           

Se puede observar que las variables han sido transformadas en factor, se puede detallar valores de variables que se encuentran fuera del contexto de descripción observado en la tabla de referencia anterior, se procederá al tratamiento de las mismas. Se identificará unicidad de ID lo cual indica que 1 cliente aparece sólo 1 vez y no hay duplicidad de clientes

# Verificando unicidad de ID
length(data$ID)

[1] 30000

n_distinct(data$ID)

[1] 30000

En efecto, los clientes son únicos. Acorde a las descripciones de la fuente de datos, los valores de la variable EDUCATION poseen los valores 0=sin escolaridad; 1 =escuela; 2 = universidad; 3 = bachillerato; por tanto,los valores 5,6 que aparecen en la BD serán reemplazados por la categoría 2 - universidad. Para esta variable, se codificará con los valores ordinales definiendo, así, el orden de los valores.

Los valores MARRIAGE sólo poseen las descripciones Estado Civil (1 = casado; 2 = soltero; 3 = otros) el valor 0 no pertenece a ninguna categoría, por tanto, se colocará dentro de la categoría otros

La variable sexo será transformada a los valores Género del cliente (1 = masculino; 2 = femenino)

data<- data %>% 
  mutate(EDUCATION = recode(EDUCATION,'0'='Sin Escolaridad','4'='Universidad','5'='Universidad','6'='Universidad','1'='Escuela','2'='Universidad','3'='Bachillerato')) %>% 
  mutate(SEX = recode(SEX,'1'='Masculino','2'='Femenino')) %>% 
  mutate(MARRIAGE = recode(MARRIAGE,'0'='Otros','1'='Casado','2'='Soltero','3'='Otros'))

# Establecer EDUCATION categórica ordinal

data$EDUCATION <- factor(data$EDUCATION, 
                         levels = c("Sin Escolaridad","Escuela","Bachillerato","Universidad"),
                         ordered = TRUE)

levels(data$EDUCATION)

[1] "Sin Escolaridad" "Escuela"         "Bachillerato"    "Universidad"    

La variable PAY_i posee rangos de valores defasados por 1 unidad, se incrementarán en 1 para ajustar acorde a la descripción de la base de datos. PAY_0 será renombrado como PAY_1 como inicio de pago

for(i in 7:12){
  data[i]<-data[i]+1
}

data <- data %>% rename(PAY_1 = PAY_0)

min(data$PAY_2)

[1] -1

max(data$PAY_2)

[1] 9

Se revisarán si existen valores nulos en los datos

sum(apply(data,2, function(x) is.na(x)))

[1] 0

El total de registros con NA es 0, por tanto, no es necesario hacer limpieza de valores NA

Se chequeará revisión de correlación entre las variables, determinando aquellas que posean un índice mayor de 80%. Primero, se definen los tipos de variables

continuas <- c("LIMIT_BAL","AGE","BILL_AMT1",   "BILL_AMT2",    "BILL_AMT3","BILL_AMT4",    "BILL_AMT5",    "BILL_AMT6","PAY_AMT1", "PAY_AMT2", "PAY_AMT3","PAY_AMT4",  "PAY_AMT5", "PAY_AMT6")
discretas <- c("PAY_1", "PAY_2",    "PAY_3",
               "PAY_4", "PAY_5",    "PAY_6")
categoricas <- c("SEX","EDUCATION","MARRIAGE")

target <- "default payment next month"

res <- cor(data[c(continuas,discretas)])
round(res, 2)

          LIMIT_BAL   AGE BILL_AMT1 BILL_AMT2 BILL_AMT3 BILL_AMT4 BILL_AMT5
LIMIT_BAL      1.00  0.14      0.29      0.28      0.28      0.29      0.30
AGE            0.14  1.00      0.06      0.05      0.05      0.05      0.05
BILL_AMT1      0.29  0.06      1.00      0.95      0.89      0.86      0.83
BILL_AMT2      0.28  0.05      0.95      1.00      0.93      0.89      0.86
BILL_AMT3      0.28  0.05      0.89      0.93      1.00      0.92      0.88
BILL_AMT4      0.29  0.05      0.86      0.89      0.92      1.00      0.94
BILL_AMT5      0.30  0.05      0.83      0.86      0.88      0.94      1.00
BILL_AMT6      0.29  0.05      0.80      0.83      0.85      0.90      0.95
PAY_AMT1       0.20  0.03      0.14      0.28      0.24      0.23      0.22
PAY_AMT2       0.18  0.02      0.10      0.10      0.32      0.21      0.18
PAY_AMT3       0.21  0.03      0.16      0.15      0.13      0.30      0.25
PAY_AMT4       0.20  0.02      0.16      0.15      0.14      0.13      0.29
PAY_AMT5       0.22  0.02      0.17      0.16      0.18      0.16      0.14
PAY_AMT6       0.22  0.02      0.18      0.17      0.18      0.18      0.16
PAY_1         -0.27 -0.04      0.19      0.19      0.18      0.18      0.18
PAY_2         -0.30 -0.05      0.23      0.24      0.22      0.22      0.22
PAY_3         -0.29 -0.05      0.21      0.24      0.23      0.23      0.23
PAY_4         -0.27 -0.05      0.20      0.23      0.24      0.25      0.24
PAY_5         -0.25 -0.05      0.21      0.23      0.24      0.27      0.27
PAY_6         -0.24 -0.05      0.21      0.23      0.24      0.27      0.29
          BILL_AMT6 PAY_AMT1 PAY_AMT2 PAY_AMT3 PAY_AMT4 PAY_AMT5 PAY_AMT6 PAY_1
LIMIT_BAL      0.29     0.20     0.18     0.21     0.20     0.22     0.22 -0.27
AGE            0.05     0.03     0.02     0.03     0.02     0.02     0.02 -0.04
BILL_AMT1      0.80     0.14     0.10     0.16     0.16     0.17     0.18  0.19
BILL_AMT2      0.83     0.28     0.10     0.15     0.15     0.16     0.17  0.19
BILL_AMT3      0.85     0.24     0.32     0.13     0.14     0.18     0.18  0.18
BILL_AMT4      0.90     0.23     0.21     0.30     0.13     0.16     0.18  0.18
BILL_AMT5      0.95     0.22     0.18     0.25     0.29     0.14     0.16  0.18
BILL_AMT6      1.00     0.20     0.17     0.23     0.25     0.31     0.12  0.18
PAY_AMT1       0.20     1.00     0.29     0.25     0.20     0.15     0.19 -0.08
PAY_AMT2       0.17     0.29     1.00     0.24     0.18     0.18     0.16 -0.07
PAY_AMT3       0.23     0.25     0.24     1.00     0.22     0.16     0.16 -0.07
PAY_AMT4       0.25     0.20     0.18     0.22     1.00     0.15     0.16 -0.06
PAY_AMT5       0.31     0.15     0.18     0.16     0.15     1.00     0.15 -0.06
PAY_AMT6       0.12     0.19     0.16     0.16     0.16     0.15     1.00 -0.06
PAY_1          0.18    -0.08    -0.07    -0.07    -0.06    -0.06    -0.06  1.00
PAY_2          0.22    -0.08    -0.06    -0.06    -0.05    -0.04    -0.04  0.67
PAY_3          0.22     0.00    -0.07    -0.05    -0.05    -0.04    -0.04  0.57
PAY_4          0.24    -0.01     0.00    -0.07    -0.04    -0.03    -0.03  0.54
PAY_5          0.26    -0.01     0.00     0.01    -0.06    -0.03    -0.02  0.51
PAY_6          0.29     0.00    -0.01     0.01     0.02    -0.05    -0.03  0.47
          PAY_2 PAY_3 PAY_4 PAY_5 PAY_6
LIMIT_BAL -0.30 -0.29 -0.27 -0.25 -0.24
AGE       -0.05 -0.05 -0.05 -0.05 -0.05
BILL_AMT1  0.23  0.21  0.20  0.21  0.21
BILL_AMT2  0.24  0.24  0.23  0.23  0.23
BILL_AMT3  0.22  0.23  0.24  0.24  0.24
BILL_AMT4  0.22  0.23  0.25  0.27  0.27
BILL_AMT5  0.22  0.23  0.24  0.27  0.29
BILL_AMT6  0.22  0.22  0.24  0.26  0.29
PAY_AMT1  -0.08  0.00 -0.01 -0.01  0.00
PAY_AMT2  -0.06 -0.07  0.00  0.00 -0.01
PAY_AMT3  -0.06 -0.05 -0.07  0.01  0.01
PAY_AMT4  -0.05 -0.05 -0.04 -0.06  0.02
PAY_AMT5  -0.04 -0.04 -0.03 -0.03 -0.05
PAY_AMT6  -0.04 -0.04 -0.03 -0.02 -0.03
PAY_1      0.67  0.57  0.54  0.51  0.47
PAY_2      1.00  0.77  0.66  0.62  0.58
PAY_3      0.77  1.00  0.78  0.69  0.63
PAY_4      0.66  0.78  1.00  0.82  0.72
PAY_5      0.62  0.69  0.82  1.00  0.82
PAY_6      0.58  0.63  0.72  0.82  1.00

Para acceder visualmente a las correlaciones más altas, se realizará un mapa de calor

col<- colorRampPalette(c("blue", "white", "red"))(20)
heatmap(x = res, col = col, symm = TRUE)

Se puede observar tanto en la matriz de correlación, como en el mapa de calor, que las variables BILL_AMT1 a BILL_AMT6 se encuentran altamente correlacionadas, índice de correlación >0.8 se quitarán del modelo previamente en el tratamiento de recipes

summary(data)

       ID          LIMIT_BAL              SEX                  EDUCATION    
 Min.   :    1   Min.   :  10000   Masculino:11888   Sin Escolaridad:   14  
 1st Qu.: 7501   1st Qu.:  50000   Femenino :18112   Escuela        :10585  
 Median :15000   Median : 140000                     Bachillerato   : 4917  
 Mean   :15000   Mean   : 167484                     Universidad    :14484  
 3rd Qu.:22500   3rd Qu.: 240000                                            
 Max.   :30000   Max.   :1000000                                            
    MARRIAGE          AGE            PAY_1             PAY_2        
 Otros  :  377   Min.   :21.00   Min.   :-1.0000   Min.   :-1.0000  
 Casado :13659   1st Qu.:28.00   1st Qu.: 0.0000   1st Qu.: 0.0000  
 Soltero:15964   Median :34.00   Median : 1.0000   Median : 1.0000  
                 Mean   :35.49   Mean   : 0.9833   Mean   : 0.8662  
                 3rd Qu.:41.00   3rd Qu.: 1.0000   3rd Qu.: 1.0000  
                 Max.   :79.00   Max.   : 9.0000   Max.   : 9.0000  
     PAY_3             PAY_4             PAY_5             PAY_6        
 Min.   :-1.0000   Min.   :-1.0000   Min.   :-1.0000   Min.   :-1.0000  
 1st Qu.: 0.0000   1st Qu.: 0.0000   1st Qu.: 0.0000   1st Qu.: 0.0000  
 Median : 1.0000   Median : 1.0000   Median : 1.0000   Median : 1.0000  
 Mean   : 0.8338   Mean   : 0.7793   Mean   : 0.7338   Mean   : 0.7089  
 3rd Qu.: 1.0000   3rd Qu.: 1.0000   3rd Qu.: 1.0000   3rd Qu.: 1.0000  
 Max.   : 9.0000   Max.   : 9.0000   Max.   : 9.0000   Max.   : 9.0000  
   BILL_AMT1         BILL_AMT2        BILL_AMT3         BILL_AMT4      
 Min.   :-165580   Min.   :-69777   Min.   :-157264   Min.   :-170000  
 1st Qu.:   3559   1st Qu.:  2985   1st Qu.:   2666   1st Qu.:   2327  
 Median :  22382   Median : 21200   Median :  20089   Median :  19052  
 Mean   :  51223   Mean   : 49179   Mean   :  47013   Mean   :  43263  
 3rd Qu.:  67091   3rd Qu.: 64006   3rd Qu.:  60165   3rd Qu.:  54506  
 Max.   : 964511   Max.   :983931   Max.   :1664089   Max.   : 891586  
   BILL_AMT5        BILL_AMT6          PAY_AMT1         PAY_AMT2      
 Min.   :-81334   Min.   :-339603   Min.   :     0   Min.   :      0  
 1st Qu.:  1763   1st Qu.:   1256   1st Qu.:  1000   1st Qu.:    833  
 Median : 18105   Median :  17071   Median :  2100   Median :   2009  
 Mean   : 40311   Mean   :  38872   Mean   :  5664   Mean   :   5921  
 3rd Qu.: 50191   3rd Qu.:  49198   3rd Qu.:  5006   3rd Qu.:   5000  
 Max.   :927171   Max.   : 961664   Max.   :873552   Max.   :1684259  
    PAY_AMT3         PAY_AMT4         PAY_AMT5           PAY_AMT6       
 Min.   :     0   Min.   :     0   Min.   :     0.0   Min.   :     0.0  
 1st Qu.:   390   1st Qu.:   296   1st Qu.:   252.5   1st Qu.:   117.8  
 Median :  1800   Median :  1500   Median :  1500.0   Median :  1500.0  
 Mean   :  5226   Mean   :  4826   Mean   :  4799.4   Mean   :  5215.5  
 3rd Qu.:  4505   3rd Qu.:  4013   3rd Qu.:  4031.5   3rd Qu.:  4000.0  
 Max.   :896040   Max.   :621000   Max.   :426529.0   Max.   :528666.0  
 default payment next month
 0:23364                   
 1: 6636                   
                           
                           
                           
                           

Se realizarán estudios gráficos de las variables para determinar posibles outliers, distribuciones, comportamientos en los datos

#---------------------------------------
# Histogramas

plot_hist <- function(x)
{
  ggplot(data,aes(.data[[x]])) +
    geom_histogram(aes(y=..density..),colour="darkcyan", fill="white") +
    geom_density(alpha=0.08,fill="darkslategray1",colour="aquamarine2",linewidth=0.6)
}

names(continuas) <- continuas
names(discretas) <- discretas

lapply(c(continuas,discretas), plot_hist)

$LIMIT_BAL

$AGE

$BILL_AMT1

$BILL_AMT2

$BILL_AMT3

$BILL_AMT4

$BILL_AMT5

$BILL_AMT6

$PAY_AMT1

$PAY_AMT2

$PAY_AMT3

$PAY_AMT4

$PAY_AMT5

$PAY_AMT6

$PAY_1

$PAY_2

$PAY_3

$PAY_4

$PAY_5

$PAY_6

Las variables en general, son sesgadas a derecha para las variables continuas dado que se confirma en los summary de datos que la mediana < media por lo cual confirma el sesgo a derecha. Para las variables discretas PAY_i i=1,6, se determina el comportamiento de conteo, con mayores datos representados en 1 mes de atraso para los meses de Abril - Sep 2005.

Se testearán con otros bines del histograma para observar en detalle.

Histogramas a división de 50 bines

plot_hist_50 <- function(x)
{
  ggplot(data,aes(.data[[x]])) +
    geom_histogram(bins=50,aes(y=..density..),colour="darkcyan", fill="white") +
    geom_density(alpha=0.08,fill="darkslategray1",colour="aquamarine2",linewidth=0.6)
}


lapply(c(continuas,discretas), plot_hist_50)

$LIMIT_BAL

$AGE

$BILL_AMT1

$BILL_AMT2

$BILL_AMT3

$BILL_AMT4

$BILL_AMT5

$BILL_AMT6

$PAY_AMT1

$PAY_AMT2

$PAY_AMT3

$PAY_AMT4

$PAY_AMT5

$PAY_AMT6

$PAY_1

$PAY_2

$PAY_3

$PAY_4

$PAY_5

$PAY_6

Histogramas a división de 15 bines

plot_hist_15 <- function(x)
{
  ggplot(data,aes(.data[[x]])) +
    geom_histogram(bins=15,aes(y=..density..),colour="darkcyan", fill="white") +
    geom_density(alpha=0.08,fill="darkslategray1",colour="aquamarine2",linewidth=0.6)
}

lapply(c(continuas,discretas), plot_hist_15)

$LIMIT_BAL

$AGE

$BILL_AMT1

$BILL_AMT2

$BILL_AMT3

$BILL_AMT4

$BILL_AMT5

$BILL_AMT6

$PAY_AMT1

$PAY_AMT2

$PAY_AMT3

$PAY_AMT4

$PAY_AMT5

$PAY_AMT6

$PAY_1

$PAY_2

$PAY_3

$PAY_4

$PAY_5

$PAY_6

Se observan los mismos patrones de sesgo a derecha para las variables a diferentes bines de división. Para el caso de las variables discretas, se puede observar un mejor agrupamiento de los datos numéricos en los bines = 15, dado que la representación muestra un sesgo a derecha de los datos.

Para detección de outliers, se graficarán los boxplots

plot_box <- function(x)
{
  ggplot(data,aes(`default payment next month`,.data[[x]])) +
    geom_boxplot(fill = "white", colour = "cadetblue3",outlier.colour = "deeppink", outlier.shape = 1)}


lapply(c(continuas,discretas), plot_box)

$LIMIT_BAL

$AGE

$BILL_AMT1

$BILL_AMT2

$BILL_AMT3

$BILL_AMT4

$BILL_AMT5

$BILL_AMT6

$PAY_AMT1

$PAY_AMT2

$PAY_AMT3

$PAY_AMT4

$PAY_AMT5

$PAY_AMT6

$PAY_1

$PAY_2

$PAY_3

$PAY_4

$PAY_5

$PAY_6

Con el fin de estudiar relaciones existentes con la variable dependiente, se puede observar que no hay diferencia estadística representativa entre las variables númericas dado que las cajas se solapan entre ellas afirman que no hay diferencia significativa, excepto por la variable límite, que presenta un leve distanciamiento, lo cual parece indicar que a mayores límites de TC aprobados, menor es la propensión a default, mientras que límites inferiores de aprobación de TC pueden tender a casos de default. Los outliers son en su gran mayoría presentes en los datos, lo cual puede indicar un mercado volátil en Taiwán o clientes con comportamientos atípicos en la entidad financiera que no se estandariza hacia un comportamiento con agrupamientos naturales. Para el caso de las variables PAY_AMT_i, se realizarán zooms para las variables donde los boxplots han quedado más pequeños y no se pueden visualizar correctamente dada la alta concentración de valores outliers

data_box <- data %>% 
  select(`default payment next month`,PAY_AMT1,PAY_AMT2,PAY_AMT3,PAY_AMT4,PAY_AMT5,PAY_AMT6) %>% 
  filter(PAY_AMT1 <= 5000 & 
           PAY_AMT2 <= 5000 & 
           PAY_AMT3 <= 5000 & 
           PAY_AMT4 <= 5000 & 
           PAY_AMT5 <= 5000 & 
           PAY_AMT6 <= 5000)

min(data_box$PAY_AMT6)

[1] 0

max(data_box$PAY_AMT6)

[1] 5000

plot_box_pay <- function(x)
{
  ggplot(data_box,aes(`default payment next month`,.data[[x]])) +
    geom_boxplot(fill = "white", colour = "cadetblue3",outlier.colour = "deeppink", outlier.shape = 1)
  }


pay_amount_cols<-colnames(data_box[2:7])

names(pay_amount_cols) <- pay_amount_cols

lapply(pay_amount_cols, plot_box_pay)

$PAY_AMT1

$PAY_AMT2

$PAY_AMT3

$PAY_AMT4

$PAY_AMT5

$PAY_AMT6

Hacemos facets para visualizar los histogramas con respecto a la target, en estas variables que se harán enfoque

plot_hist_target <- function(x)
{
  ggplot(data_box,aes(.data[[x]],fill=`default payment next month`)) +
    geom_histogram()+
    facet_wrap(~ `default payment next month`,ncol=1)
  }

lapply(pay_amount_cols, plot_hist_target)

$PAY_AMT1

$PAY_AMT2

$PAY_AMT3

$PAY_AMT4

$PAY_AMT5

$PAY_AMT6

Se pueden observar valores de sesgo a derecha para las variables de PAY_AMT_i los valores más altos de montos pagados por los clientes se observan en promedio para los meses de Abril - Mayo, Julio - Septiembre, el mes de Junio posee una reducción sobre los montos pagados respecto a los valores promedio de los otros meses.

Para el caso de las variables PAY_AMT1 y PAY_AMT2 son estadísticamente significativas respecto a la variable de default de TC dado que no hay solapamiento completo de las cajas, en PAY_AMT3 hay un pequeño distanciamiento pero no es tan evidente como los distanciamientos de PAY_AMT1 y PAY_AMT2, esto quiere decir que, los meses de Abril y Mayo 2005 son representativos para influir sobre los valores de la target.

Estudio de pairplots para ver relaciones entre las variables respecto a la target.

Se seleccionarán para realizar las gráficas de scatter plot las variables contínuas AGE, LIMIT_BAL con BILL_AMOUNT_1,2 y BILL_AMOUNT_5,BILL_AMOUNT_6. En un 2° gráfico, se visualizarán las variables contínuas AGE, LIMIT_BAL con PAY_AMOUNT_1,2 y PAY_AMOUNT_5,PAY_AMOUNT_6 para mejor visualización de los datos

names(target) <- target


continuas2<-c("LIMIT_BAL","AGE","BILL_AMT1","BILL_AMT2","BILL_AMT5",    "BILL_AMT6")
continuas3<-c("LIMIT_BAL","AGE","PAY_AMT1", "PAY_AMT2","PAY_AMT5",  "PAY_AMT6")

names(continuas2) <- continuas2
names(continuas3) <- continuas3


pairs(data[continuas2],
      col = c( "cornflowerblue", "purple")[data$`default payment next month`])

Se pueden observar relaciones lineales entre las variables BILL_AMT_i, hecho que se puede determinar desde la matriz de correlación donde las correlaciones son mayores al 80%

pairs(data[continuas3],
      col = c( "cornflowerblue", "purple")[data$`default payment next month`])

Los valores de PAY_AMT_i con LIMIT_BAL y AGE no poseen relación lineal. No existen separaciones de grupos representativos por target entre las variables cuantitativas.

data_scatter1 <- data[c(discretas,target)]

pairs(data_scatter1,
      col = c( "cornflowerblue", "purple")[data_scatter1$`default payment next month`])

Se visualizan sobre las variables PAY_3, PAY_4, PAY_5, PAY_6 una ligera relación lineal, efectos que pueden visualizarse en la correlación entre dichas variables con correlaciones promedio de 55%.

#---------------------------------------
# Bar plots


plot_bar <- function(x)
{
    ggplot(data = data,aes(x = .data[[x]], fill = `default payment next month`)) +
    geom_bar(alpha=0.6) +
    scale_fill_brewer(palette="Pastel2") +
    geom_text(stat = 'count', aes(label=..count..), position=position_dodge(width = 1)
              ,vjust=-0.5, size=3,color="midnightblue")
}

names(categoricas) <- categoricas


lapply(categoricas, plot_bar)

$SEX

$EDUCATION

$MARRIAGE

La población se distribuye en mayor volumen a mujeres que a hombres, a nivel porcentual, la tasa de default en hombres es del 24% comparado con la tasa de default en mujeres que representa un 21%.

A nivel educativo, la población posee en su mayoría grados universitarios seguidos de escuela, a nivel relativo, las tasas de default son mayores para las personas con bachillerato representando un 25%, mientras que los universitarios poseen una tasa del 23% y las personas con escolaridad representan la menor tasa del 19%, esto podría inferir que existe un mercado informal que maneja negocios propios y posee un manejo adecuado del dinero apesar de poseer sólo grado escolar.

En cuanto al estado civil, se puede observar que más de 15.000 clientes se distribuye en estado civil soltero, porcentualmente, las tasas de default mayor ocurren para el estado civil de casado con un 23%, mientras que para los solteros, la tasa es de 21% se podría inferir que al adquirir mayores responsabilidades como hijos, hipotecas podría implicar un default mayor para pagos de tarjetas de crédito.

#---------------------------------------
# Polar coord


plot_polar <- function(x)
{
  bar <- ggplot(data = data) +
    geom_bar(
      mapping = aes(x = .data[[x]], fill = .data[[x]]),
      show.legend = FALSE,
      width = 1
    ) +
    scale_fill_brewer(palette="Pastel2")+
    theme(aspect.ratio = 1) +
    labs(x = NULL, y = NULL) +
    ggtitle(x) +
    theme(plot.title = element_text(size=10, face="bold"))
  bar + coord_flip()
  bar + coord_polar()
  
}


lapply(c(categoricas,target), plot_polar)

$SEX

$EDUCATION

$MARRIAGE

$`default payment next month`

table(data$`default payment next month`)

    0     1 
23364  6636 

4 Fase de Modelado - Modelo GLM

Se procede al ajuste del modelo, se puede observar que la variable respuesta tiene una distribución en menor proporción de los clientes con default vs. los clientes sin default, la relación de datos es de un 22% para los casos de interés y 78% para los casos sin default, la relación puede ser considerada para la aplicación del modelo sin realizar reajuste de balanceo de la muestra, la variable target es una variable dicotómica, por tal razón, se deben aplicar modelos de clasificación supervisados lo que indica que son modelos predictivos en base a los valores de default o no en TC. Se aplicará un modelo lineal generalizado GLM con función vínculo logit para evaluar los ajustes a la fuente de datos. Se utilizará un método denominado ‘recipes’ de la librería tidymodels que ajusta los pasos a través de secuencias ordenadas

Previo al ajuste, se seleccionarán las muestras de entrenamiento y test

# Fijar semilla para reproducción de resultados
set.seed(123)

# Corte a un 80% training, 20$ test
data_split <- rsample::initial_split(data,prop=0.8)
data_train <- rsample::training(data_split)
data_test <- rsample::testing(data_split)

# Ajuste con cross validation

set.seed(123)
k<-5
folds<-rsample::vfold_cv(data_train,v=k)

# Ajuste de recipe

recip1 <- recipe(`default payment next month`~., data=data_train) %>%  # Se consideran todas las variables del modelo con la data de train
  update_role(ID,new_role = "ID") %>%  # Se identifica que la variable ID es un ID y por tanto no se toma en cuenta como variable predictora
  step_corr(all_numeric_predictors(), threshold = .8) %>%  # Elimina variables con más del 80% de correlación
  step_normalize(all_numeric_predictors()) %>%  # Estandarización de valores para comparación de escalas entre variables numéricas
  step_dummy(all_nominal_predictors(), one_hot=FALSE) %>% # Transformación de las variables categóricas nominales, usa internamente la decodificación de las clases en binarias quitando una de las variables y ordinales que al ser declarada ordinal previamente se asigna un contrast_poly un ajuste polinómico para dar orden 
  step_zv(all_predictors()) # Remueve variables que sólo contengan 1 valor en los datos, por ello se aplica luego del paso step_dummy

print(summary(recip1),n=25)

# A tibble: 25 × 4
   variable                   type      role      source  
   <chr>                      <list>    <chr>     <chr>   
 1 ID                         <chr [2]> ID        original
 2 LIMIT_BAL                  <chr [2]> predictor original
 3 SEX                        <chr [3]> predictor original
 4 EDUCATION                  <chr [2]> predictor original
 5 MARRIAGE                   <chr [3]> predictor original
 6 AGE                        <chr [2]> predictor original
 7 PAY_1                      <chr [2]> predictor original
 8 PAY_2                      <chr [2]> predictor original
 9 PAY_3                      <chr [2]> predictor original
10 PAY_4                      <chr [2]> predictor original
11 PAY_5                      <chr [2]> predictor original
12 PAY_6                      <chr [2]> predictor original
13 BILL_AMT1                  <chr [2]> predictor original
14 BILL_AMT2                  <chr [2]> predictor original
15 BILL_AMT3                  <chr [2]> predictor original
16 BILL_AMT4                  <chr [2]> predictor original
17 BILL_AMT5                  <chr [2]> predictor original
18 BILL_AMT6                  <chr [2]> predictor original
19 PAY_AMT1                   <chr [2]> predictor original
20 PAY_AMT2                   <chr [2]> predictor original
21 PAY_AMT3                   <chr [2]> predictor original
22 PAY_AMT4                   <chr [2]> predictor original
23 PAY_AMT5                   <chr [2]> predictor original
24 PAY_AMT6                   <chr [2]> predictor original
25 default payment next month <chr [3]> outcome   original

Las variables predictoras y target son separadas durante el proceso de steps de recipe, lo que garantiza que no se está incluyendo como variable independiente, al igual que la variable ID

mod_glm <- logistic_reg() %>% 
  set_engine("glm")

workflow_glm <- workflow() %>% 
  add_model(mod_glm) %>% 
  add_recipe(recip1)

keep_pred <- control_resamples(save_pred = TRUE)
perf_metr<-metric_set(specificity, recall, accuracy,precision)

cv_glm <- fit_resamples(workflow_glm,
                resamples=folds,
                metrics = perf_metr,
                control = keep_pred)

collect_metrics(cv_glm)

# A tibble: 4 × 6
  .metric     .estimator  mean     n  std_err .config             
  <chr>       <chr>      <dbl> <int>    <dbl> <chr>               
1 accuracy    binary     0.810     5 0.00388  Preprocessor1_Model1
2 precision   binary     0.818     5 0.00386  Preprocessor1_Model1
3 recall      binary     0.972     5 0.000732 Preprocessor1_Model1
4 specificity binary     0.244     5 0.00790  Preprocessor1_Model1

fit_glm<-fit(workflow_glm, data_train)
fit_glm

══ Workflow [trained] ══════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: logistic_reg()

── Preprocessor ────────────────────────────────────────────────────────────────
4 Recipe Steps

• step_corr()
• step_normalize()
• step_dummy()
• step_zv()

── Model ───────────────────────────────────────────────────────────────────────

Call:  stats::glm(formula = ..y ~ ., family = stats::binomial, data = data)

Coefficients:
     (Intercept)         LIMIT_BAL               AGE             PAY_1  
        -4.21897          -0.10913           0.05682           0.65391  
           PAY_2             PAY_3             PAY_6         BILL_AMT1  
         0.09819           0.12796           0.05043          -0.12684  
        PAY_AMT1          PAY_AMT2          PAY_AMT3          PAY_AMT4  
        -0.19650          -0.13445          -0.03803          -0.09156  
        PAY_AMT5          PAY_AMT6      SEX_Femenino       EDUCATION_1  
        -0.05751          -0.04130          -0.12245           7.20050  
     EDUCATION_2       EDUCATION_3   MARRIAGE_Casado  MARRIAGE_Soltero  
        -5.41583           2.52395           0.22545           0.04276  

Degrees of Freedom: 23999 Total (i.e. Null);  23980 Residual
Null Deviance:      25460 
Residual Deviance: 22350    AIC: 22390

print(tidy(fit_glm), n=27)

# A tibble: 20 × 5
   term             estimate std.error statistic   p.value
   <chr>               <dbl>     <dbl>     <dbl>     <dbl>
 1 (Intercept)       -4.22     21.5       -0.196 8.45e-  1
 2 LIMIT_BAL         -0.109     0.0229    -4.77  1.80e-  6
 3 AGE                0.0568    0.0191     2.97  2.96e-  3
 4 PAY_1              0.654     0.0222    29.5   2.79e-191
 5 PAY_2              0.0982    0.0268     3.66  2.49e-  4
 6 PAY_3              0.128     0.0264     4.85  1.25e-  6
 7 PAY_6              0.0504    0.0209     2.41  1.59e-  2
 8 BILL_AMT1         -0.127     0.0222    -5.72  1.09e-  8
 9 PAY_AMT1          -0.196     0.0396    -4.96  6.88e-  7
10 PAY_AMT2          -0.134     0.0424    -3.17  1.53e-  3
11 PAY_AMT3          -0.0380    0.0285    -1.34  1.82e-  1
12 PAY_AMT4          -0.0916    0.0303    -3.02  2.54e-  3
13 PAY_AMT5          -0.0575    0.0262    -2.19  2.84e-  2
14 PAY_AMT6          -0.0413    0.0250    -1.65  9.84e-  2
15 SEX_Femenino      -0.122     0.0343    -3.57  3.50e-  4
16 EDUCATION_1        7.20     57.8        0.125 9.01e-  1
17 EDUCATION_2       -5.42     43.1       -0.126 9.00e-  1
18 EDUCATION_3        2.52     19.3        0.131 8.96e-  1
19 MARRIAGE_Casado    0.225     0.151      1.50  1.34e-  1
20 MARRIAGE_Soltero   0.0428    0.152      0.281 7.79e-  1

drop1(fit_glm |> extract_fit_engine(), .~., test = "Chisq")

Single term deletions

Model:
..y ~ LIMIT_BAL + AGE + PAY_1 + PAY_2 + PAY_3 + PAY_6 + BILL_AMT1 + 
    PAY_AMT1 + PAY_AMT2 + PAY_AMT3 + PAY_AMT4 + PAY_AMT5 + PAY_AMT6 + 
    SEX_Femenino + EDUCATION_1 + EDUCATION_2 + EDUCATION_3 + 
    MARRIAGE_Casado + MARRIAGE_Soltero
                 Df Deviance   AIC    LRT  Pr(>Chi)    
<none>                 22347 22387                     
LIMIT_BAL         1    22371 22409  23.26 1.418e-06 ***
AGE               1    22356 22394   8.79 0.0030232 ** 
PAY_1             1    23214 23252 866.67 < 2.2e-16 ***
PAY_2             1    22361 22399  13.34 0.0002593 ***
PAY_3             1    22371 22409  23.33 1.367e-06 ***
PAY_6             1    22353 22391   5.82 0.0158879 *  
BILL_AMT1         1    22381 22419  33.54 6.985e-09 ***
PAY_AMT1          1    22382 22420  34.15 5.108e-09 ***
PAY_AMT2          1    22360 22398  12.60 0.0003850 ***
PAY_AMT3          1    22349 22387   1.94 0.1634766    
PAY_AMT4          1    22359 22397  11.22 0.0008109 ***
PAY_AMT5          1    22353 22391   5.46 0.0194274 *  
PAY_AMT6          1    22350 22388   2.93 0.0872075 .  
SEX_Femenino      1    22360 22398  12.74 0.0003579 ***
EDUCATION_1       1    22352 22390   4.09 0.0431580 *  
EDUCATION_2       1    22352 22390   4.46 0.0347038 *  
EDUCATION_3       1    22354 22392   6.48 0.0109315 *  
MARRIAGE_Casado   1    22350 22388   2.31 0.1283877    
MARRIAGE_Soltero  1    22348 22386   0.08 0.7779573    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El modelo en base a la devianza D^2=((25460-22350)/25460)*100=12,45% el modelo explica el 12,45% de la variabilidad. Se observan que las variables PAY_AMT3, PAY_AMT6, las variables asociadas a MARRIAGE no son significativas p_valor_variable > 0.05, esto implica que no se puede rechazar H0: beta_i =0, por tanto, dichas variables no aportan al modelo, se eliminarán estas variables del modelo para determinar si se observan mejoría en las métricas de performance del modelo, sin embargo, se observa que el coeficiente del criterio de información de Akaike AIC no disminuye si se elimina algunas de las variables del modelo, al no presentar de reducción de AIC, parece no ser recomendada la eliminación, a modo de testing se realiza la exclusión de variables no significativas

recip2 <- recipe(`default payment next month`~ LIMIT_BAL+AGE+PAY_1+ PAY_2+  PAY_3+PAY_6+
                   BILL_AMT1+PAY_AMT1+  PAY_AMT2+PAY_AMT4+PAY_AMT5+SEX+
                   EDUCATION, data=data_train) %>% 
  step_normalize(all_numeric_predictors()) %>% 
  step_dummy(all_nominal_predictors(), one_hot=FALSE) %>% 
  step_zv(all_predictors())

mod_glm2 <- logistic_reg() %>% 
  set_engine("glm")

workflow_glm2 <- workflow() %>% 
  add_model(mod_glm2) %>% 
  add_recipe(recip2)


cv_glm <- fit_resamples(workflow_glm2,
                resamples=folds,
                metrics = perf_metr,
                control = keep_pred)

collect_metrics(cv_glm)

# A tibble: 4 × 6
  .metric     .estimator  mean     n  std_err .config             
  <chr>       <chr>      <dbl> <int>    <dbl> <chr>               
1 accuracy    binary     0.808     5 0.00348  Preprocessor1_Model1
2 precision   binary     0.817     5 0.00355  Preprocessor1_Model1
3 recall      binary     0.972     5 0.000446 Preprocessor1_Model1
4 specificity binary     0.240     5 0.00705  Preprocessor1_Model1

fit_glm2<-fit(workflow_glm2, data_train)

fit_glm2

══ Workflow [trained] ══════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: logistic_reg()

── Preprocessor ────────────────────────────────────────────────────────────────
3 Recipe Steps

• step_normalize()
• step_dummy()
• step_zv()

── Model ───────────────────────────────────────────────────────────────────────

Call:  stats::glm(formula = ..y ~ ., family = stats::binomial, data = data)

Coefficients:
 (Intercept)     LIMIT_BAL           AGE         PAY_1         PAY_2  
    -4.11315      -0.10646       0.09567       0.65653       0.09786  
       PAY_3         PAY_6     BILL_AMT1      PAY_AMT1      PAY_AMT2  
     0.12864       0.04908      -0.13401      -0.20534      -0.14322  
    PAY_AMT4      PAY_AMT5  SEX_Femenino   EDUCATION_1   EDUCATION_2  
    -0.09695      -0.06465      -0.10813       7.25043      -5.43472  
 EDUCATION_3  
     2.52289  

Degrees of Freedom: 23999 Total (i.e. Null);  23984 Residual
Null Deviance:      25460 
Residual Deviance: 22380    AIC: 22410

print(tidy(fit_glm2), n=27)

# A tibble: 16 × 5
   term         estimate std.error statistic   p.value
   <chr>           <dbl>     <dbl>     <dbl>     <dbl>
 1 (Intercept)   -4.11     21.5       -0.191 8.49e-  1
 2 LIMIT_BAL     -0.106     0.0225    -4.73  2.19e-  6
 3 AGE            0.0957    0.0171     5.61  2.07e-  8
 4 PAY_1          0.657     0.0222    29.6   7.54e-193
 5 PAY_2          0.0979    0.0268     3.65  2.59e-  4
 6 PAY_3          0.129     0.0264     4.88  1.07e-  6
 7 PAY_6          0.0491    0.0209     2.35  1.87e-  2
 8 BILL_AMT1     -0.134     0.0220    -6.09  1.09e-  9
 9 PAY_AMT1      -0.205     0.0398    -5.17  2.40e-  7
10 PAY_AMT2      -0.143     0.0425    -3.37  7.58e-  4
11 PAY_AMT4      -0.0970    0.0305    -3.18  1.49e-  3
12 PAY_AMT5      -0.0646    0.0264    -2.45  1.43e-  2
13 SEX_Femenino  -0.108     0.0341    -3.17  1.51e-  3
14 EDUCATION_1    7.25     57.8        0.125 9.00e-  1
15 EDUCATION_2   -5.43     43.1       -0.126 9.00e-  1
16 EDUCATION_3    2.52     19.3        0.131 8.96e-  1

drop1(fit_glm2 |> extract_fit_engine(), .~., test = "Chisq")

Single term deletions

Model:
..y ~ LIMIT_BAL + AGE + PAY_1 + PAY_2 + PAY_3 + PAY_6 + BILL_AMT1 + 
    PAY_AMT1 + PAY_AMT2 + PAY_AMT4 + PAY_AMT5 + SEX_Femenino + 
    EDUCATION_1 + EDUCATION_2 + EDUCATION_3
             Df Deviance   AIC    LRT  Pr(>Chi)    
<none>             22376 22408                     
LIMIT_BAL     1    22399 22429  22.87 1.733e-06 ***
AGE           1    22407 22437  31.21 2.314e-08 ***
PAY_1         1    23250 23280 873.58 < 2.2e-16 ***
PAY_2         1    22390 22420  13.27 0.0002698 ***
PAY_3         1    22400 22430  23.63 1.169e-06 ***
PAY_6         1    22382 22412   5.53 0.0186821 *  
BILL_AMT1     1    22415 22445  38.24 6.259e-10 ***
PAY_AMT1      1    22414 22444  37.47 9.306e-10 ***
PAY_AMT2      1    22391 22421  14.46 0.0001434 ***
PAY_AMT4      1    22389 22419  12.53 0.0004015 ***
PAY_AMT5      1    22383 22413   6.92 0.0085203 ** 
SEX_Femenino  1    22386 22416  10.04 0.0015352 ** 
EDUCATION_1   1    22381 22411   4.38 0.0364592 *  
EDUCATION_2   1    22381 22411   4.61 0.0316990 *  
EDUCATION_3   1    22383 22413   6.45 0.0110779 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El modelo en base a la devianza D^2=((25460-22380)/25460)*100=12,10% el modelo explica el 12,10% de la variabilidad. El coeficiente de AIC aumenta lo cual indica que el modelo saturado (que incluye todas las variables) es más significativo dado que presenta menor AIC y explica un poco más en comparación con el modelo reducido

Al analizar los valores de AIC del primer modelo, se puede determinar que el AIC no posee un descenso al retirar cualquier variable, esto indica que quitar alguna variable no significativa no implica la mejora de AIC, se verifica con un test anova para comparar ambos modelos 1 y 2

anova(fit_glm2 |> extract_fit_engine(),fit_glm |> extract_fit_engine(), test = "Chisq")

Analysis of Deviance Table

Model 1: ..y ~ LIMIT_BAL + AGE + PAY_1 + PAY_2 + PAY_3 + PAY_6 + BILL_AMT1 + 
    PAY_AMT1 + PAY_AMT2 + PAY_AMT4 + PAY_AMT5 + SEX_Femenino + 
    EDUCATION_1 + EDUCATION_2 + EDUCATION_3
Model 2: ..y ~ LIMIT_BAL + AGE + PAY_1 + PAY_2 + PAY_3 + PAY_6 + BILL_AMT1 + 
    PAY_AMT1 + PAY_AMT2 + PAY_AMT3 + PAY_AMT4 + PAY_AMT5 + PAY_AMT6 + 
    SEX_Femenino + EDUCATION_1 + EDUCATION_2 + EDUCATION_3 + 
    MARRIAGE_Casado + MARRIAGE_Soltero
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1     23984      22376                          
2     23980      22347  4   28.817 8.517e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El test de ANOVA indica con un p_valor de 1.283e-05 se puede rechazar H0 lo que indica que el modelo con menor cantidad de variables no es adecuado para explicar la variable dependiente

Se procederá a realizar el análisis de residuos

#Chequeo de normalidad de los residuos

hist(rstandard(fit_glm |> extract_fit_engine()), freq=FALSE, col="red2",xlab="Residuos estandarizados", ylab="Frecuencia", main="Histograma de los errores estandarizados")

boxplot(rstandard(fit_glm |> extract_fit_engine()),main="Residuos estandarizados",ylab="Residuos estandarizados",col="lightgoldenrodyellow")

qqnorm(rstandard(fit_glm |> extract_fit_engine()),main="Gr?fico Cuantil-Cuantil")
qqline(rstandard(fit_glm |> extract_fit_engine()),lwd=2,col="red")

#H0: los residuos son normales
jarque.bera.test(residuals(fit_glm |> extract_fit_engine()))

    Jarque Bera Test

data:  residuals(extract_fit_engine(fit_glm))
X-squared = 7177.3, df = 2, p-value < 2.2e-16

En el test de Jarque Bera podemos rechazar H0 lo cual indica que los residuos no son normales lo cual es aceptable en un modelo lineal generalizado GLM donde no es exigible el cumplimiento de normalidad de los residuos

#Chequeo de Homocedasticidad
plot(fitted.values(fit_glm |> extract_fit_engine()),rstandard(fit_glm |> extract_fit_engine()),xlab="Valores Ajustados",ylab="Residuos Estandarizados",main="Residuos Estandarizados vs valores ajustados")
abline(0,0,col="red",lwd=2)

#H0: los residuos son homocedásticos
breusch_pagan(fit_glm |> extract_fit_engine())

# A tibble: 1 × 5
  statistic   p.value parameter method                alternative
      <dbl>     <dbl>     <dbl> <chr>                 <chr>      
1      698. 9.00e-136        19 Koenker (studentised) greater    

De igual forma, los residuos no son homocedásticos, al rechazar H0: los residuos son homocedásticos en el test de Breusch Pagan, condición aceptada en un modelo de regresión logística

#Chequeo de Independencia
plot(rstandard(fit_glm |> extract_fit_engine()),main="Residuos estandarizados")
abline(0,0,col="red",lwd=2)

Los residuos se encuentran en el rango entre -3 y 3, se encuentran aglomerados en la parte inferior lo cual indica que no son independientes, en el caso de los modelos GLM esta condición debe garantizarse.

#Chequeo de Autocorrelación de los residuos
#H0: los residuos no son autocorrelacionados

dwtest(fit_glm |> extract_fit_engine(),alternative = "two.sided",data=data_train)

    Durbin-Watson test

data:  extract_fit_engine(fit_glm)
DW = 1.9955, p-value = 0.7295
alternative hypothesis: true autocorrelation is not 0

El p_valor = 0.8128 es demasiado alto como para rechazar H0, por tanto, los residuos no son autocorrelacionados

5 Evaluación - Modelo GLM

Se valuarán los resultados del modelo tanto para la data train como el data test para evaluar el performance y descartar overfitting o underfitting en los resultados

# Resultados data train - GLM
pred_train_glm <- predict(fit_glm, data_train)
pred_test_glm <- predict(fit_glm, data_test)

glm_aug <- augment(fit_glm,data_train)

# Visual de las clases predichas con su probabilidad

glm_aug %>% 
  select(`default payment next month`,.pred_class,.pred_1)

# A tibble: 24,000 × 3
   `default payment next month` .pred_class .pred_1
   <fct>                        <fct>         <dbl>
 1 0                            0            0.180 
 2 0                            0            0.187 
 3 0                            0            0.230 
 4 0                            0            0.213 
 5 0                            0            0.102 
 6 0                            0            0.215 
 7 1                            0            0.0573
 8 0                            0            0.219 
 9 1                            0            0.233 
10 0                            0            0.269 
# ℹ 23,990 more rows

Se obtienen tanto las clases predichas 0, 1 como las probabilidades de default, en este caso, para la data train.

Se procede a realizar el cálculo de ROC AUC sobre data train

# ROC data train
glm_aug %>% 
  roc_auc(truth = `default payment next month`,.pred_0)

# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 roc_auc binary         0.726

glm_aug %>% 
  roc_curve(truth = `default payment next month`,.pred_0) %>% 
  autoplot()

result_train <- glm_aug %>% 
          roc_curve(truth = `default payment next month`,.pred_0)

Se procede a realizar el cálculo de ROC AUC sobre data test

# Resultados data test - GLM

glm_aug2 <- augment(fit_glm,data_test)

glm_aug2 %>% 
  roc_auc(truth = `default payment next month`,.pred_0)

# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 roc_auc binary         0.713

glm_aug2 %>% 
  roc_curve(truth = `default payment next month`,.pred_0) %>% 
  autoplot()

result_test <- glm_aug2 %>% 
          roc_curve(truth = `default payment next month`,.pred_0)

El valor del área bajo la curva ROC-AUC es de 73% sobre los datos de train, el ROC para los datos de test son 71% el porcentaje sobre los datos de test hacen que no se evidencie un comportamiento de bajo performance sobre los datos de test, obteniéndose resultados muy similares de performance, esto implica, que no hay indicios de overfitting ( o de underfitting) en el modelo acorde a los valores de ROC -AUC. Se evaluarán, análogamente, los valores de métricas como especificidad, sensitividad, acuracidad y precisión

# Confusion Matrix - data train

glm_aug %>% 
  conf_mat(truth = `default payment next month`, estimate = .pred_class) %>% 
  autoplot(type = 'heatmap')

# Confusion Matrix - data test

glm_aug2 %>% 
  conf_mat(truth = `default payment next month`, estimate = .pred_class) %>% 
  autoplot(type = 'heatmap')

df_tr_glm<- data.frame(default_payment=data_train$`default payment next month`,preds=pred_train_glm$.pred_class)
df_test_glm<- data.frame(default_payment=data_test$`default payment next month`,preds=pred_test_glm$.pred_class)

dim(df_tr_glm)

[1] 24000     2

dim(df_test_glm)

[1] 6000    2

eval_metrics <- metric_set(specificity, recall, accuracy,precision)

# Métricas para data train

eval_metrics(data = df_tr_glm , truth = default_payment , estimate =preds)

# A tibble: 4 × 3
  .metric     .estimator .estimate
  <chr>       <chr>          <dbl>
1 specificity binary         0.246
2 recall      binary         0.972
3 accuracy    binary         0.811
4 precision   binary         0.818

# Métricas para data test

eval_metrics(data = df_test_glm , truth = default_payment , estimate =preds)

# A tibble: 4 × 3
  .metric     .estimator .estimate
  <chr>       <chr>          <dbl>
1 specificity binary         0.232
2 recall      binary         0.972
3 accuracy    binary         0.813
4 precision   binary         0.822

Para el caso de default en TC, la clase positiva son los clientes que realizan default target =1, en la matriz de confusión son los casos con True Positives = 1.314 y 299 para data train y test, respectivamente, la clase negativa son los clientes que no hicieron default en TC, en la matriz de confusión son los casos con True Negatives = 18.141 y 4.579 target =0 para data train y test, respectivamente, para el caso de considerar más relevante en el estudio un falso negativo indica que el cliente realiza default pero el modelo lo predice como un cliente que no entra en default, esto hace que los FN sean más importantes que los FP en el contexto del significado de la definición de sensitividad = TP / TP+FN entre más pequeño sea el valor de FN, más grande será el valor de la sensitividad, por tanto, se coloca mayor enfásis en que el modelo tenga un valor más alto de recall que de precisión. Los valores de acuracidad son del 81%, sensitividad 97%, precisión 82% y especificidad de 25% para data train y 81% acuracidad, sensitividad 97%, precisión 82% y especificidad de 23% para data test, lo cual garantiza valores de performance similares para el data train y data test descartando la posibilidad de overfitting o de underfitting en los resultados del modelo.

Se calculan 2 coeficientes adicionales de ajuste el coeficiente de Kolmogorov Smirnov que mide la distancia máxima entre las clases de la variable dependiente, entre mayor sea el coeficiente, mejor el ajuste, viene definido por

\[ KS = \max_{a \in [L, H]} \left| F_{n.buenos}(a) - F_{m.malos}(a) \right| \]

El coeficiente de Gini que mide la separación entre las clases de la target viene calculada con la siguiente fórmula

\[ GINI = 1 - \sum_{k=2}^{n+m} [F_{m.malos_k} - F_{m.malos_{k-1}}][F_{n.buenos_k} + F_{n.buenos_{k-1}}] \]

gini_train = 2*(0.726) -1
gini_train

[1] 0.452

gini_test = 2*(0.713) -1
gini_test

[1] 0.426

diff_TPRFPR = result_train[['sensitivity']]-(1-result_train[['specificity']])
KS_train = max(diff_TPRFPR)
KS_train

[1] 0.3755179

diff_TPRFPR_test = result_test[['sensitivity']]-(1-result_test[['specificity']])
KS_test = max(diff_TPRFPR_test)
KS_test

[1] 0.3784625

Los valores de KS y GINI podrían representar mayores valores de discriminación entre clases. Para ello, se buscará implementar la opción a través XGBoost. Para los valores entre train y test no se observan diferencias significativas, lo cual descarta presencia de overfitting o de underfitting.

6 Fase de Modelado - XGBoost

Se ajustará un XGBoost que consiste en una serie de árboles de decisión con boosting donde cada árbol depende de los resultados del árbol previo. Todos los árboles en el ensamble son combinados para producir un resultado final, se detectarán los cores posibles de uso de la máquina para paralelizar la búsqueda cuando se realice el tunning de parámetros, en este caso, con hyperparameter tuning que consiste en buscar los parámatros que mejor optimizan el valor de una métrica obtenida, lo cual forma una grilla de resultados y se elije el mejor

cores <- parallel::detectCores()

cores

[1] 16

set.seed(123)
k<-5
folds<-rsample::vfold_cv(data_train,v=k)

recip_rf <- recipe(`default payment next month`~., data=data_train) %>% 
  update_role(ID,new_role = "ID") %>%  # Se identifica que la variable ID es un ID y por tanto no se toma en cuenta como variable predictora
  step_corr(all_numeric_predictors(), threshold = .8) %>%  # Elimina variables con más del 80% de correlación
  step_normalize(all_numeric_predictors()) %>%  # Estandarización de valores para comparación de escalas entre variables numéricas
  step_dummy(all_nominal_predictors(), one_hot=FALSE) %>% # Transformación de las variables categóricas nominales, usa internamente la decodificación de las clases en binarias quitando una de las variables y ordinales que al ser declarada ordinal previamente se asigna un contrast_poly un ajuste polinómico para dar orden 
  step_zv(all_predictors()) # Remueve variables que sólo contengan 1 valor en los datos, por ello se aplica luego del paso step_dummy


rf_mod <- 
  boost_tree(tree_depth = tune(),trees = tune()) %>%  # tree_depth y trees hiperparámetros a tunear
  set_engine("xgboost", num.threads = cores) %>% 
  set_mode("classification")

rf_workflow <- 
  workflow() %>% 
  add_model(rf_mod) %>% 
  add_recipe(recip_rf)


set.seed(123)
rf_res <- 
  rf_workflow %>% 
  tune_grid(folds,
            grid = 25,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))

rf_res

# Tuning results
# 5-fold cross-validation 
# A tibble: 5 × 5
  splits               id    .metrics          .notes           .predictions
  <list>               <chr> <list>            <list>           <list>      
1 <split [19200/4800]> Fold1 <tibble [25 × 6]> <tibble [0 × 3]> <tibble>    
2 <split [19200/4800]> Fold2 <tibble [25 × 6]> <tibble [0 × 3]> <tibble>    
3 <split [19200/4800]> Fold3 <tibble [25 × 6]> <tibble [0 × 3]> <tibble>    
4 <split [19200/4800]> Fold4 <tibble [25 × 6]> <tibble [0 × 3]> <tibble>    
5 <split [19200/4800]> Fold5 <tibble [25 × 6]> <tibble [0 × 3]> <tibble>    

autoplot(rf_res) + 
  scale_color_viridis_d(direction = -1) + 
  theme(legend.position = "top")

rf_best <- 
  rf_res %>% 
  select_best(metric = "roc_auc")
rf_best

# A tibble: 1 × 3
  trees tree_depth .config              
  <int>      <int> <chr>                
1  1748          1 Preprocessor1_Model21

Los valores del grid search con máx entriopía se encuentran separados, en efecto, los mejores parámteros son trees = 1748, tree_depth = 1

rf_res %>% 
  collect_predictions(parameters = rf_best)

# A tibble: 24,000 × 8
   id    .pred_0 .pred_1  .row trees tree_depth default payment next m…¹ .config
   <chr>   <dbl>   <dbl> <int> <int>      <int> <fct>                    <chr>  
 1 Fold1   0.872  0.128      3  1748          1 0                        Prepro…
 2 Fold1   0.927  0.0729     5  1748          1 0                        Prepro…
 3 Fold1   0.834  0.166      7  1748          1 1                        Prepro…
 4 Fold1   0.856  0.144      8  1748          1 0                        Prepro…
 5 Fold1   0.879  0.121     11  1748          1 0                        Prepro…
 6 Fold1   0.717  0.283     12  1748          1 1                        Prepro…
 7 Fold1   0.391  0.609     18  1748          1 1                        Prepro…
 8 Fold1   0.685  0.315     28  1748          1 1                        Prepro…
 9 Fold1   0.800  0.200     31  1748          1 0                        Prepro…
10 Fold1   0.887  0.113     32  1748          1 0                        Prepro…
# ℹ 23,990 more rows
# ℹ abbreviated name: ¹​`default payment next month`

7 Evaluación Modelo - XGBoost

rf_mod <- 
  boost_tree(tree_depth = 1,trees = 1748) %>% 
  set_engine("xgboost", num.threads = cores) %>% 
  set_mode("classification")


rf_workflow <- 
  workflow() %>% 
  add_model(rf_mod) %>% 
  add_recipe(recip_rf)

# Curva ROC-AUC Train
fit_rf_train<- fit(rf_workflow, data_train)

[01:27:36] WARNING: src/learner.cc:767: 
Parameters: { "num_threads" } are not used.

pred_train_rf <- predict(fit_rf_train, data_train)
pred_test_rf <- predict(fit_rf_train, data_test)

rf_aug <- augment(fit_rf_train,data_train)

rf_aug %>% 
  select(`default payment next month`,.pred_class,.pred_1)

# A tibble: 24,000 × 3
   `default payment next month` .pred_class .pred_1
   <fct>                        <fct>         <dbl>
 1 0                            0            0.103 
 2 0                            0            0.0823
 3 0                            0            0.135 
 4 0                            0            0.0987
 5 0                            0            0.0671
 6 0                            0            0.110 
 7 1                            0            0.166 
 8 0                            0            0.135 
 9 1                            0            0.142 
10 0                            0            0.156 
# ℹ 23,990 more rows

rf_aug %>% 
  roc_auc(truth = `default payment next month`,.pred_0)

# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 roc_auc binary         0.797

rf_aug %>% 
  roc_curve(truth = `default payment next month`,.pred_0) %>% 
  autoplot()

result_train_rf <- rf_aug %>% 
  roc_curve(truth = `default payment next month`,.pred_0)

# Curva ROC-AUC Test
rf_aug_test <- augment(fit_rf_train,data_test)

rf_aug_test %>% 
  select(`default payment next month`,.pred_class,.pred_1)

# A tibble: 6,000 × 3
   `default payment next month` .pred_class .pred_1
   <fct>                        <fct>         <dbl>
 1 1                            1            0.788 
 2 0                            0            0.174 
 3 0                            0            0.101 
 4 0                            0            0.148 
 5 0                            0            0.124 
 6 0                            0            0.189 
 7 0                            0            0.0410
 8 0                            0            0.0615
 9 0                            0            0.148 
10 0                            0            0.0760
# ℹ 5,990 more rows

rf_aug_test %>% 
  roc_auc(truth = `default payment next month`,.pred_0)

# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 roc_auc binary         0.767

rf_aug_test %>% 
  roc_curve(truth = `default payment next month`,.pred_0) %>% 
  autoplot()

result_test_rf <- rf_aug_test %>% 
  roc_curve(truth = `default payment next month`,.pred_0)

Los valores de ROC AUC para train están en un valor del 80% respecto a un 77% de la data de testing, los valores son muy similares respecto a la data de testing, se evaluarán las métricas de acuracidad, recall, precisión

df_tr<- data.frame(default_payment=data_train$`default payment next month`,preds=pred_train_rf$.pred_class)
df_test<- data.frame(default_payment=data_test$`default payment next month`,preds=pred_test_rf$.pred_class)

head(df_tr)

  default_payment preds
1               0     0
2               0     0
3               0     0
4               0     0
5               0     0
6               0     0

head(df_test)

  default_payment preds
1               1     1
2               0     0
3               0     0
4               0     0
5               0     0
6               0     0

dim(df_tr)

[1] 24000     2

dim(df_test)

[1] 6000    2

eval_metrics <- metric_set(specificity, recall, accuracy,precision)
eval_metrics(data = df_tr , truth = default_payment , estimate =preds)

# A tibble: 4 × 3
  .metric     .estimator .estimate
  <chr>       <chr>          <dbl>
1 specificity binary         0.356
2 recall      binary         0.954
3 accuracy    binary         0.820
4 precision   binary         0.838

eval_metrics(data = df_test , truth = default_payment , estimate =preds)

# A tibble: 4 × 3
  .metric     .estimator .estimate
  <chr>       <chr>          <dbl>
1 specificity binary         0.356
2 recall      binary         0.951
3 accuracy    binary         0.823
4 precision   binary         0.843

Los valores de acuracidad son 82% en data train vs. 82.3% data test, recall 95.4% train vs. 95.1% test, precisión 84.8% train vs. 84.3% test, especificidad 35.6% train y test, las métricas no son distantes entre sí lo cual indica un buen performance del modelo, se evaluarán los coeficientes GINI y KS

gini_train = 2*(0.797) -1
gini_train

[1] 0.594

gini_test = 2*(0.767) -1
gini_test

[1] 0.534

diff_TPRFPR = result_train_rf[['sensitivity']]-(1-result_train_rf[['specificity']])
KS_train = max(diff_TPRFPR)
KS_train

[1] 0.4415802

diff_TPRFPR_test = result_test_rf[['sensitivity']]-(1-result_test_rf[['specificity']])
KS_test = max(diff_TPRFPR_test)
KS_test

[1] 0.4138844

Finalmente, el modelo de XGBoost presenta mejores métricas de ajuste para los datos el GINI y el KS fueron 59% train, 53% test y 44% train, 41% test, respectivamente. En general, el modelo de XGBoost posee mejor ajuste en comparación con el GLM, no se presentan a nivel general de los indicadores indicios de overfitting o de underfitting sobre el performance del modelo. Se representan las variables que poseen mayor representatividad de importancia en el XGBoost con el top 20

fit_rf_train %>% 
  extract_fit_parsnip() %>% 
  vip(num_features = 20)

Las variables que mayor impacto generan en la predicción del default en tarjetas de crédito son asociados a las cantidades, montos, días de pago, así como, estatus civil, edad y límites aprobados.