Capítulo 10 Trabajo Practico

10.1 Reglas del Trabajo practico

  1. Integrantes: máximo 3 por grupo.

  2. Extensión: máximo 8 carillas (hoja A4, 12pts, etc.). La página 9 no se corrige.

  3. Formato: se debe realizar utilizando Markdown y deben entregar un archivo .pdf y el .Rmd

  4. Copia o plagio: trabajo desaprobado, grupo fuera del régimen de promoción.

  5. Redacción: Formal.

  6. Presentación: tablas/cuadros bien descriptas y ordenadas.

  7. Bases de datos: a definir diferenciando por grupos.

  8. Sugerencias bibliográficas:

10.2 Enunciado del Trabajo Practico

En base a lo desarrollado en las clases teóricas se busca que elaboren un modelo de scoring que permita discriminar entre buenos y malos deudores.

  1. Realizar una revisión de la literatura teórica y empírica para elaborar una sección que describa:
  • ¿Por qué es importante el problema a analizar?
  • Ventajas y desventajas de enfoques tradicionales vs. machine learning.
  • Los principales resultados encontrados.

  1. Presentar, describir y analizar los datos utilizados.

  2. Presentar e interpretar los principales resultados de un modelo logit en comparación con técnicas de árboles.

  3. Elaborar conclusiones de política.

10.3 Aplicacion practica

La irrupción de las firmas BigTech en la provisión de crédito está modificando la estructura del sistema financiero. Si bien la actividad principal de estas compañías es la provisión de servicios digitales como el e-commerce y servicios de pago, paulatinamente han ido incorporando otros productos como la provisión de crédito, seguros, inversiones y ahorro.

El modelo de negocios de las BigTech difiere del modelo de las entidades financieras tradicionales principalmente por dos factores distintivos: efectos de red (generados por las plataformas de e-commerce, aplicaciones de mensajería y redes sociales); el uso de la tecnología (inteligencia artificial utilizando big data).

La utilización de nuevas técnicas de análisis y fuentes de datos alternativos brindan a las empresas tecnológicas una ventaja informativa para la evaluación de deudores respecto de las entidades financieras, que utilizan métodos econométricos convencionales (ej. estimaciones logit) menos flexibles para capturar la información contenida en grandes volúmenes de datos.

En esta sección se utiliza una base de datos bancaria para predecir la probabilidad de default con distintas metodologías, un modelo logit, un árbol simple y un random forest, para comparar las capacidades predictivas.

Primero se limpia la memoria y se cargan las librerías que vamos a utilizar.

# Limpiar memoria
rm(list=ls())
gc()
##           used  (Mb) gc trigger  (Mb) max used  (Mb)
## Ncells 3138965 167.7    6226840 332.6  4049954 216.3
## Vcells 5511369  42.1   12255594  93.6 10146179  77.5
# Librerias
library(tidyverse)
library(rsample)
library(yardstick)
library(rpart)
library(rpart.plot)
library(ranger)
library(caret)

Se cargan los datos desde un archivo separado por comas y se analizan los valores missing. En la Figura de abajo las variables están en el eje y y las observaciones en el eje x.

# Cargar datos
score_data_raw <- read_csv('./data/CreditScore1.csv', show_col_types = FALSE)

score_data_raw %>%
  is.na() %>%
  reshape2::melt() %>%
  ggplot(aes(Var2, Var1, fill=value)) + 
  geom_raster() + coord_flip() +
  scale_y_continuous(NULL, expand = c(0, 0)) +
  scale_fill_discrete(name = "", 
                  labels = c("Dato", 
                             "Missing")) +
  labs( x= 'Variable', title= 'Valores missing') +
  theme_minimal() +
  theme(axis.text.y  = element_text(size = 7)) +
  NULL

La variable ingreso es la que mayormente presenta observaciones con valores missing, dado que no existe una forma obvia de imputación y que tenemos 150.000 observaciones en total (y no se pierden tantas…), se eliminan las filas correspondientes.

score_data_tbl <- score_data_raw %>%
  dplyr::select(-id) %>% drop_na()

Se realiza una inspección inicial de la base de datos.

head(score_data_tbl)
## # A tibble: 6 × 10
##   default personal_total  edad nro_atraso3059 gastos_ingreso ingreso
##     <dbl>          <dbl> <dbl>          <dbl>          <dbl>   <dbl>
## 1       1          0.766    45              2         0.803     9120
## 2       0          0.957    40              0         0.122     2600
## 3       0          0.658    38              1         0.0851    3042
## 4       0          0.234    30              0         0.0360    3300
## 5       0          0.907    49              1         0.0249   63588
## 6       0          0.213    74              0         0.376     3500
## # ℹ 4 more variables: lineas_credito <dbl>, nro_atraso90 <dbl>,
## #   nro_hipoteca <dbl>, familia <dbl>
glimpse(score_data_tbl)
## Rows: 120,269
## Columns: 10
## $ default        <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1…
## $ personal_total <dbl> 0.76612663, 0.95715100, 0.65818012, 0.23380977, 0.90723…
## $ edad           <dbl> 45, 40, 38, 30, 49, 74, 39, 57, 30, 51, 46, 40, 76, 64,…
## $ nro_atraso3059 <dbl> 2, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0…
## $ gastos_ingreso <dbl> 0.80298215, 0.12187620, 0.08511338, 0.03604968, 0.02492…
## $ ingreso        <dbl> 9120, 2600, 3042, 3300, 63588, 3500, 3500, 23684, 2500,…
## $ lineas_credito <dbl> 13, 4, 2, 5, 7, 3, 8, 9, 5, 7, 13, 9, 6, 7, 7, 7, 2, 10…
## $ nro_atraso90   <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0…
## $ nro_hipoteca   <dbl> 6, 0, 0, 0, 1, 1, 0, 4, 0, 2, 2, 1, 1, 1, 0, 1, 0, 2, 1…
## $ familia        <dbl> 2, 1, 0, 0, 0, 1, 0, 2, 0, 2, 2, 2, 0, 2, 0, 2, 0, 0, 2…
summary(score_data_tbl)
##     default        personal_total          edad        nro_atraso3059   
##  Min.   :0.00000   Min.   :    0.00   Min.   :  0.00   Min.   : 0.0000  
##  1st Qu.:0.00000   1st Qu.:    0.04   1st Qu.: 40.00   1st Qu.: 0.0000  
##  Median :0.00000   Median :    0.18   Median : 51.00   Median : 0.0000  
##  Mean   :0.06949   Mean   :    5.90   Mean   : 51.29   Mean   : 0.3818  
##  3rd Qu.:0.00000   3rd Qu.:    0.58   3rd Qu.: 61.00   3rd Qu.: 0.0000  
##  Max.   :1.00000   Max.   :50708.00   Max.   :103.00   Max.   :98.0000  
##  gastos_ingreso        ingreso        lineas_credito    nro_atraso90    
##  Min.   :    0.00   Min.   :      0   Min.   : 0.000   Min.   : 0.0000  
##  1st Qu.:    0.14   1st Qu.:   3400   1st Qu.: 5.000   1st Qu.: 0.0000  
##  Median :    0.30   Median :   5400   Median : 8.000   Median : 0.0000  
##  Mean   :   26.60   Mean   :   6670   Mean   : 8.758   Mean   : 0.2119  
##  3rd Qu.:    0.48   3rd Qu.:   8249   3rd Qu.:11.000   3rd Qu.: 0.0000  
##  Max.   :61106.50   Max.   :3008750   Max.   :58.000   Max.   :98.0000  
##   nro_hipoteca       familia       
##  Min.   : 0.000   Min.   : 0.0000  
##  1st Qu.: 0.000   1st Qu.: 0.0000  
##  Median : 1.000   Median : 0.0000  
##  Mean   : 1.055   Mean   : 0.8518  
##  3rd Qu.: 2.000   3rd Qu.: 2.0000  
##  Max.   :54.000   Max.   :20.0000
sort(unique(score_data_tbl$familia))
##  [1]  0  1  2  3  4  5  6  7  8  9 10 13 20
table(score_data_tbl$default)
## 
##      0      1 
## 111912   8357
# Porcentaje de positivos
8357 / (8357 + 111912)
## [1] 0.0694859

Luego, se calculan algunas estadísticas descriptivas.

stat = score_data_tbl %>% 
  dplyr::select_if(is.numeric) %>% 
  pivot_longer(everything(), names_to = 'Variable', values_to = 'Value') %>%
  group_by(Variable) %>% 
  summarise(
    Obs = n(),
    Media = mean(Value, na.rm = T),
    Mediana = median(Value, na.rm = T),
    SD = sd(Value, na.rm = T),
    Min = min(Value, na.rm = T),
    Max = max(Value, na.rm = T)) %>% 
  ungroup()
stat
## # A tibble: 10 × 7
##    Variable          Obs     Media  Mediana        SD   Min      Max
##    <chr>           <int>     <dbl>    <dbl>     <dbl> <dbl>    <dbl>
##  1 default        120269    0.0695    0         0.254     0       1 
##  2 edad           120269   51.3      51        14.4       0     103 
##  3 familia        120269    0.852     0         1.15      0      20 
##  4 gastos_ingreso 120269   26.6       0.296   424.        0   61106.
##  5 ingreso        120269 6670.     5400     14385.        0 3008750 
##  6 lineas_credito 120269    8.76      8         5.17      0      58 
##  7 nro_atraso3059 120269    0.382     0         3.50      0      98 
##  8 nro_atraso90   120269    0.212     0         3.47      0      98 
##  9 nro_hipoteca   120269    1.05      1         1.15      0      54 
## 10 personal_total 120269    5.90      0.177   257.        0   50708

Se procede a realizar el feature engineering o creación de variables…notar que la capacidad de clasificación depende de los atributos disponibles y los valores de los hiperpárametros.23

# Transformacion
score_data_tbl = score_data_tbl %>%
                 mutate(
                    default = factor(default),
                    ingreso = log(1+ingreso), 
                    familia_bin = case_when(familia == 0 ~ 1,
                                            familia == 1 ~ 2,
                                            familia >= 2&familia < 5 ~ 3,
                                            familia >= 5 ~ 4))
score_data_tbl = score_data_tbl %>% dplyr::select(-familia)

Se divide la muestra en \(80\%\) para entrenamiento y \(20\%\) para test.

# Train / Test split
set.seed(1234)
train_test_split <- initial_split(score_data_tbl, prop = 0.8)
train_test_split
## <Training/Testing/Total>
## <96215/24054/120269>
train_tbl <- training(train_test_split)
test_tbl  <- testing(train_test_split)

Se definen dos objetos para utilizar más abajo.

# Formula
formula  <-  formula(default ~ .)

# Y observado a 0/1 para confusionMatrix
obs =  factor(test_tbl$default)

Se estima el modelo lineal (default vs. resto de variables).

lm.mod = lm(as.numeric(default) ~ ., data = train_tbl)
summary(lm.mod)
## 
## Call:
## lm(formula = as.numeric(default) ~ ., data = train_tbl)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.90421 -0.08651 -0.06467 -0.03990  1.18753 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1.157e+00  6.296e-03 183.764  < 2e-16 ***
## personal_total -3.760e-06  3.623e-06  -1.038 0.299402    
## edad           -1.585e-03  5.894e-05 -26.887  < 2e-16 ***
## nro_atraso3059  2.144e-02  1.059e-03  20.242  < 2e-16 ***
## gastos_ingreso -5.645e-06  1.962e-06  -2.877 0.004010 ** 
## ingreso        -2.315e-03  6.831e-04  -3.389 0.000701 ***
## lineas_credito -6.870e-04  1.776e-04  -3.868 0.000110 ***
## nro_atraso90   -1.354e-02  1.071e-03 -12.651  < 2e-16 ***
## nro_hipoteca    1.685e-03  7.908e-04   2.131 0.033066 *  
## familia_bin     7.086e-03  9.798e-04   7.232 4.78e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2504 on 96205 degrees of freedom
## Multiple R-squared:  0.02673,    Adjusted R-squared:  0.02664 
## F-statistic: 293.6 on 9 and 96205 DF,  p-value: < 2.2e-16

Alternativamente:

broom::tidy(lm.mod)
## # A tibble: 10 × 5
##    term              estimate  std.error statistic   p.value
##    <chr>                <dbl>      <dbl>     <dbl>     <dbl>
##  1 (Intercept)     1.16       0.00630       184.   0        
##  2 personal_total -0.00000376 0.00000362     -1.04 2.99e-  1
##  3 edad           -0.00158    0.0000589     -26.9  1.20e-158
##  4 nro_atraso3059  0.0214     0.00106        20.2  6.41e- 91
##  5 gastos_ingreso -0.00000564 0.00000196     -2.88 4.01e-  3
##  6 ingreso        -0.00232    0.000683       -3.39 7.01e-  4
##  7 lineas_credito -0.000687   0.000178       -3.87 1.10e-  4
##  8 nro_atraso90   -0.0135     0.00107       -12.7  1.19e- 36
##  9 nro_hipoteca    0.00169    0.000791        2.13 3.31e-  2
## 10 familia_bin     0.00709    0.000980        7.23 4.78e- 13
t(broom::glance(lm.mod))
##                        [,1]
## r.squared      2.673354e-02
## adj.r.squared  2.664249e-02
## sigma          2.504424e-01
## statistic      2.936161e+02
## p.value        0.000000e+00
## df             9.000000e+00
## logLik        -3.305975e+03
## AIC            6.633949e+03
## BIC            6.738167e+03
## deviance       6.034112e+03
## df.residual    9.620500e+04
## nobs           9.621500e+04

Se estima el modelo logit.24

glm.mod <- glm(formula, data = train_tbl, family = binomial)
summary(glm.mod)
## 
## Call:
## glm(formula = formula, family = binomial, data = train_tbl)
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -1.142e+00  1.007e-01 -11.341  < 2e-16 ***
## personal_total -8.713e-05  9.148e-05  -0.953 0.340843    
## edad           -2.770e-02  1.027e-03 -26.979  < 2e-16 ***
## nro_atraso3059  2.548e-01  1.393e-02  18.288  < 2e-16 ***
## gastos_ingreso -2.669e-04  7.694e-05  -3.469 0.000523 ***
## ingreso        -4.414e-02  1.181e-02  -3.738 0.000186 ***
## lineas_credito -9.410e-03  3.021e-03  -3.114 0.001843 ** 
## nro_atraso90   -2.120e-01  1.415e-02 -14.982  < 2e-16 ***
## nro_hipoteca    4.254e-02  1.204e-02   3.533 0.000410 ***
## familia_bin     1.265e-01  1.478e-02   8.555  < 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: 48424  on 96214  degrees of freedom
## Residual deviance: 46551  on 96205  degrees of freedom
## AIC: 46571
## 
## Number of Fisher Scoring iterations: 6
# Efectos marginales ver:
# library(mfx)
# logitmfx(formula, data)

glm.probs <- predict(glm.mod, test_tbl,  type = 'response')
glm.class <- factor(ifelse(glm.probs > 0.5, 1, 0))

cm_logit = confusionMatrix(glm.class, obs, positive = '1')

cm_logit
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction     0     1
##          0 22347  1679
##          1    11    17
##                                           
##                Accuracy : 0.9297          
##                  95% CI : (0.9264, 0.9329)
##     No Information Rate : 0.9295          
##     P-Value [Acc > NIR] : 0.4463          
##                                           
##                   Kappa : 0.0175          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.0100236       
##             Specificity : 0.9995080       
##          Pos Pred Value : 0.6071429       
##          Neg Pred Value : 0.9301174       
##              Prevalence : 0.0705080       
##          Detection Rate : 0.0007067       
##    Detection Prevalence : 0.0011640       
##       Balanced Accuracy : 0.5047658       
##                                           
##        'Positive' Class : 1               
##

Se estima un árbol simple.

set.seed(4321)
rpart.mod = rpart(formula,
                  data = train_tbl,  
                  control = rpart.control(minsplit = 20,  # obs minimas por nodo para partir
                                          minbucket = 6,  # obs minimas en nodo terminal
                                          cp = 0,         # complexidad
                                          xval = 0,       # numero de cross-validations.
                                          maxdepth = 16)) # profundidad del arbol
names(rpart.mod)
##  [1] "frame"               "where"               "call"               
##  [4] "terms"               "cptable"             "method"             
##  [7] "parms"               "control"             "functions"          
## [10] "numresp"             "splits"              "variable.importance"
## [13] "y"                   "ordered"
rpart.prob = predict(rpart.mod, test_tbl)
rpart.class = factor(ifelse(rpart.prob[, '1']>0.5, 1, 0))

cm_rpart = confusionMatrix(rpart.class, obs, positive = '1')

cm_rpart
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction     0     1
##          0 21956  1372
##          1   402   324
##                                           
##                Accuracy : 0.9262          
##                  95% CI : (0.9229, 0.9295)
##     No Information Rate : 0.9295          
##     P-Value [Acc > NIR] : 0.9754          
##                                           
##                   Kappa : 0.2352          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.19104         
##             Specificity : 0.98202         
##          Pos Pred Value : 0.44628         
##          Neg Pred Value : 0.94119         
##              Prevalence : 0.07051         
##          Detection Rate : 0.01347         
##    Detection Prevalence : 0.03018         
##       Balanced Accuracy : 0.58653         
##                                           
##        'Positive' Class : 1               
## 
prp(rpart.mod, extra=101, digits=2, branch=1, type=4, varlen=0, faclen=0)

rpartVarImp = as_tibble_row(rpart.mod$variable.importance) %>%
 pivot_longer(everything(), names_to = 'Variable', values_to = 'Value') %>%
 arrange(desc(Value)) 
rpartVarImp
## # A tibble: 9 × 2
##   Variable        Value
##   <chr>           <dbl>
## 1 nro_atraso90   1451. 
## 2 personal_total  748. 
## 3 nro_atraso3059  436. 
## 4 gastos_ingreso  428. 
## 5 ingreso         395. 
## 6 edad            265. 
## 7 lineas_credito  226. 
## 8 nro_hipoteca    112. 
## 9 familia_bin      69.5

Se estima un random forest.

set.seed(1234)
ranger.mod  = ranger(formula,
                     data = train_tbl,  
                     probability = TRUE,      # estima una probabilidad
                     num.trees = 300,         # numero de arboles
                     min.node.size = 15,      # tamaño mínimo de nodo a partir 
                     mtry = 3,                # nro de vars para separar por cada nodo
                     splitrule = 'gini',       
                     importance = 'impurity')   
names(ranger.mod)
##  [1] "predictions"               "num.trees"                
##  [3] "num.independent.variables" "mtry"                     
##  [5] "min.node.size"             "variable.importance"      
##  [7] "prediction.error"          "forest"                   
##  [9] "splitrule"                 "treetype"                 
## [11] "call"                      "importance.mode"          
## [13] "num.samples"               "replace"
ranger.prob = predict(ranger.mod, test_tbl)
ranger.class = factor(ifelse(ranger.prob$predictions[, '1']>0.5, 1, 0))

cm_ranger = confusionMatrix(ranger.class, obs, positive = '1')

# Error del RF

1 - cm_ranger$overall[['Accuracy']]
## [1] 0.066517
rangerVarImp = as_tibble_row(ranger.mod$variable.importance) %>%
               pivot_longer(everything(), names_to = 'Variable', values_to = 'Value') %>%
               arrange(desc(Value)) 
rangerVarImp
## # A tibble: 9 × 2
##   Variable       Value
##   <chr>          <dbl>
## 1 personal_total 1741.
## 2 gastos_ingreso 1330.
## 3 ingreso        1186.
## 4 nro_atraso90   1149.
## 5 edad            839.
## 6 nro_atraso3059  611.
## 7 lineas_credito  607.
## 8 nro_hipoteca    237.
## 9 familia_bin     204.

Se pueden realizar gráficos de dependencia parcial para distintas variables. Estos gráficos ilustran el efecto marginal de las variables seleccionadas en la respuesta después de integrar las otras variables. En este caso la probabilidad de default al principio cae con el ingreso aunque luego de un determinado umbral comienza a aumentar.

library(iml)
# Se elimina la variable target
train_tbl1 = train_tbl %>% dplyr::select(-default)

# Se crea una funcion para usar como argumento mas abajo
pfun = function(object, newdata) {predict(object, data = newdata)$predictions[, 2]}

# Sentencia para el grafico de dependencia parcial
predictor = Predictor$new(ranger.mod, data = train_tbl1, y = factor(train_tbl$default), predict.fun = pfun)

pdp = FeatureEffect$new(predictor, feature = 'ingreso', method = 'pdp') |> plot()
pdp

Se presentan los resultados (no se analizan…) en tabla resumen.

tab_acc = tibble(logit  = cm_logit$overall[['Accuracy']],
                 rpart  = cm_rpart$overall[['Accuracy']],
                 ranger = cm_ranger$overall[['Accuracy']])
                
tab_acc = tab_acc |>
          pivot_longer(everything(), names_to='Modelo', values_to='Accuracy') |>
          arrange(desc(Accuracy))
tab_acc 
## # A tibble: 3 × 2
##   Modelo Accuracy
##   <chr>     <dbl>
## 1 ranger    0.933
## 2 logit     0.930
## 3 rpart     0.926

Finalmente se calculan las AUC:

tab = tibble(obs = factor(test_tbl$default),
             logit = glm.probs,
             rpart = rpart.prob[, '1'],
             ranger = ranger.prob$predictions[,'1'])

tab_auc = tibble(logit = roc_auc(tab, truth=obs, logit, event_level = 'second')$.estimate,
                 rpart = roc_auc(tab, truth=obs, rpart, event_level = 'second')$.estimate,
                 ranger = roc_auc(tab, truth=obs, ranger, event_level = 'second')$.estimate)

tab_auc = tab_auc |> 
          pivot_longer(everything(), names_to="Modelo", values_to="AUC") |>
          arrange(desc(AUC))
tab_auc
## # A tibble: 3 × 2
##   Modelo   AUC
##   <chr>  <dbl>
## 1 ranger 0.842
## 2 rpart  0.814
## 3 logit  0.648

Bibliografia

Bazarbash, Majid. 2019. «FinTech in Financial Inclusion Machine Learning Applications in Assessing Credit Risk». IMF Working Paper, n.º 109. https://www.imf.org/~/media/Files/Publications/WP/2019/WPIEA2019109.ashx.
Frost, Jon, Leonardo Gambacorta, Yi Huang, Hyun Song Shin, y Pablo Zbinden. 2019. «BigTech and the changing structure of financial intermediation». BIS Working Papers, n.º 779. https://www.bis.org/publ/work779.htm.
Petropoulos, Anastasios, Vasilis Siakoulis, Evaggelos Stavroulakis, y Aristotelis Klamargias. 2018. «A robust machine learning approach for credit risk analysis of large loan-level datasets using deep learning and extreme gradient boosting». Irving Fisher Committee. https://www.bis.org/ifc/publ/ifcb49_49.pdf.

  1. Es importante señalar que la estadística descriptiva sugiere realizar más modificaciones.↩︎

  2. Qué sucede si para definir la clase se modifica el umbral \(p > 0.5\)?↩︎