0.1 Introducción

En este caso práctico se se realizan predicciónes de aprobación de credito bancario del set de datos Credit Aprroval Data Set del UC Irvine Machine Learning Repository. La información sobre el set de datos se encuentra en este enlace, Se emplean los modelos de selección automática de parámetros Ridge regression y LASSO.

0.2 Configuración inicial

0.2.1 Opciones globales

knitr::opts_chunk$set(message = FALSE, warning = FALSE, fig.align = "center", table.align = "center", echo = TRUE)

0.2.2 Carga de librerías

library(tidyverse)
library(here)
library(xray)
library(gridExtra)
library(missForest)
library(forcats)
library(FSelector)
library(glmnet)
library(caret)
library(knitr)
library(kableExtra)
library(broom)
library(influence.ME)

0.3 Carga de datos

datos <-read.csv(here::here("datos", "CASO_FINAL_crx.data"))

0.4 Formateo y limpieza inicial

# renombramos columnas
encabezados <- paste0("A", 1:16)
names(datos) <- encabezados


kable(summary(datos)) %>% kable_styling(font_size = 11, full_width = T) %>% 
  row_spec(0, color = "white", background = "steelblue")
A1 
   A2 </th>
 
   A3 </th>
 A4 A5 
   A6 </th>
 
   A7 </th>
 
   A8 </th>
 A9 A10 
  A11 </th>
 A12 A13 
  A14 </th>
 
  A15 </th>
 A16
?: 12 ? : 12 Min. : 0.000 ?: 6 ? : 6 c :137 v :398 Min. : 0.000 f:329 f:395 Min. : 0.000 f:373 g:624 00000 :132 Min. : 0 -:383
a:210 22.67 : 9 1st Qu.: 1.000 l: 2 g :518 q : 78 h :138 1st Qu.: 0.165 t:360 t:294 1st Qu.: 0.000 t:316 p: 8 00120 : 35 1st Qu.: 0 +:306
b:467 20.42 : 7 Median : 2.750 u:518 gg: 2 w : 63 bb : 59 Median : 1.000 NA NA Median : 0.000 NA s: 57 00200 : 35 Median : 5 NA
NA 18.83 : 6 Mean : 4.766 y:163 p :163 i : 59 ff : 57 Mean : 2.225 NA NA Mean : 2.402 NA NA 00160 : 34 Mean : 1019 NA
NA 19.17 : 6 3rd Qu.: 7.250 NA NA aa : 54 ? : 9 3rd Qu.: 2.625 NA NA 3rd Qu.: 3.000 NA NA 00080 : 30 3rd Qu.: 396 NA
NA 20.67 : 6 Max. :28.000 NA NA ff : 53 j : 8 Max. :28.500 NA NA Max. :67.000 NA NA 00100 : 30 Max. :100000 NA
NA (Other):643 NA NA NA (Other):245 (Other): 20 NA NA NA NA NA NA (Other):393 NA NA

Observamos en primer lugar, que hay algunos datos etiquetados con “?”, que corresponden a valores perdidos, los convertiremos en la nomenclatura de valores perdidos para R

# valores perdidos están identificados con "?". Cambiamos a valor perdido por defecto de R
datos <- datos %>% 
  na_if("?")

# eliminaicon de nivel "?" de factors
datos[c("A1", "A4", "A5", "A6", "A7", "A14")] <- map(datos[c("A1", "A4", "A5", "A6", "A7", "A14")], ~fct_drop(.x, "?"))

Hay algunas variables que al contener el caracter “?” se han detectado como factores a pesar de ser numéricas. Procedemos al cambio de clase.

map_chr(datos, class)
##        A1        A2        A3        A4        A5        A6        A7 
##  "factor"  "factor" "numeric"  "factor"  "factor"  "factor"  "factor" 
##        A8        A9       A10       A11       A12       A13       A14 
## "numeric"  "factor"  "factor" "integer"  "factor"  "factor"  "factor" 
##       A15       A16 
## "integer"  "factor"
datos[c("A2", "A14")] <- map(datos[c("A2", "A14")], ~as.numeric(.x))
map_chr(datos, class)
##        A1        A2        A3        A4        A5        A6        A7 
##  "factor" "numeric" "numeric"  "factor"  "factor"  "factor"  "factor" 
##        A8        A9       A10       A11       A12       A13       A14 
## "numeric"  "factor"  "factor" "integer"  "factor"  "factor" "numeric" 
##       A15       A16 
## "integer"  "factor"

Para facilitar la obtencion de resultados tras utilizar los algoritmos de prediccion, cambiamos la nomenclatura de las categorías de la variable predictora, A16, a 0 para no aprobación de crédito, y 1 para aprobación de crédito.

datos <- datos %>% 
  mutate(A16 = as.factor(if_else(A16 == "-", 0, 1)))

0.5 Preexploración

0.5.1 Detección de anomalías

A11 y A15 tienen un alto porcentaje de ceros. En principio no es algo que deba resultar preocupante.

kable(anomalies(datos)$variables) %>% kable_styling(font_size = 11, 
                                                    bootstrap_options = c("striped", "hover", "condensed", "responsive"), 
                                                    position = "right") %>% 
  row_spec(0, color = "white", background = "steelblue") %>% 
  row_spec(c(1,3),background = "#FBF3A9") %>% 
  column_spec(1, bold = T)
Variable q qNA pNA qZero pZero qBlank pBlank qInf pInf qDistinct type anomalous_percent
A11 689 0
395 57.33% 0
0
23 Integer 57.33%
A16 689 0
383 55.59% 0
0
2 Factor 55.59%
A15 689 0
294 42.67% 0
0
240 Integer 42.67%
A8 689 0
70 10.16% 0
0
132 Numeric 10.16%
A3 689 0
18 2.61% 0
0
215 Numeric 2.61%
A14 689 13 1.89% 0
0
0
170 Numeric 1.89%
A1 689 12 1.74% 0
0
0
3 Factor 1.74%
A2 689 12 1.74% 0
0
0
349 Numeric 1.74%
A7 689 9 1.31% 0
0
0
10 Factor 1.31%
A6 689 9 1.31% 0
0
0
15 Factor 1.31%
A4 689 6 0.87% 0
0
0
4 Factor 0.87%
A5 689 6 0.87% 0
0
0
4 Factor 0.87%
A9 689 0
0
0
0
2 Factor
A10 689 0
0
0
0
2 Factor
A12 689 0
0
0
0
2 Factor
A13 689 0
0
0
0
3 Factor

0.6 EDA

0.6.1 Distribución de variables explicativas en función de variable dependiente

Usaremos una función auxiliar para graficar en varios paneles la distribución de las variables explicativas contra la variable dependiente

grafica_por_tipo <- function(datos, var) {
  
  # character -> ggplot
  
  # Dibuja un grafico dada una variable "A16" presente en un data.frame, contra otra de las variables del midmo data.frame.
  # El tipo de grafico depende del tipo de variable
  # ASUME: el data.frame contiene la variable "A16"

  # aislamiento de la variable elegida para detectar suclase
  vector_var <- datos[ ,var]

  # guardamos grafico de barras si la variable explicativa es factor
  if (is.factor(vector_var)) {
    g <- datos %>% 
      ggplot(aes_string(x = var)) + 
      geom_bar(aes(fill = A16), position = "dodge") # dividimos A16 por color

  }
  
  else {
    
    # guardamos diagrama de caja si la variable explicativa es numerica
   g <-  datos %>% 
      ggplot(aes(x = A16)) + 
     geom_boxplot(aes_string(y = var, fill = "A16")) # una caja para cada uno de los valores de A16
  }
  
  # dibujo y formateo del grafico
  g + labs(x = "", y = "", fill = "", title = var) + # etiquetas de ejes y titulo
    scale_fill_manual(values = c("orange", "steelblue")) + # colores para los valores de a16
    theme_bw() + # tema generico del grafico
    
    # establecimiento de tamaño de fuentes
    theme(plot.title = element_text(size = 9), 
          axis.text.y = element_text(size = 6))
}

# creacion de lista de graficos de A16 contra el resto de variables
lista_graficos <- map(paste0("A", 1:15), ~grafica_por_tipo(datos, .x))

# dibujo de los gráficos en un panel de 5 columnas
grid.arrange(grobs = lista_graficos, ncol = 5, 
             top = "Distribución de las variables condicionada a A16")

Las variables con más influencia sobre la dependiente parecen ser la A9 y la A10. El calculo de la ganancia de información confirma la influencia de estas dos variables. Pero aparecen también en lugar destacado las variables A11 y A15, cuya influencia pasa desapercibida en el gráfico debido a que hay presencia de valores extremos. Se da el caso de que son las mismas variables con una grán presencia de ceros.

grafica_ganancia_informacion <- function(datos) {
  
# data.frame -> ggplot
  
# grafica en un diagrama de barras la ganancia de informacion de las variables de un data.frame al desagruparlas
# por la variable A!6
# ASUME: el data.frame contiene la variable A16
  
# calculo de ganancia de informacion y formato del data.frame resultante 
information.gain(A16~., data = datos) %>% 
  mutate(variables = rownames(.)) %>% 
  select(variables, attr_importance) %>% 
  
  # generacion del grafico
  ggplot(aes(x = fct_reorder(variables, attr_importance), y = attr_importance)) + 
    
  # columnas
  geom_col(fill = "steelblue", alpha = 0.8) +
    
  # particion de columnas con linea blanca
  geom_hline(yintercept = seq(0, 0.3, by = 0.05), color = "white") + 
    
  # inversion de coordenadas para convertir columnas en barras
  coord_flip() + 
    
  # adicion de etiquetas numericas con ganancia de información
  geom_text(aes(label = round(attr_importance, 3)), nudge_y = 0.015, color = "grey45") +
    
    # Formato del grafico
    labs(x = "Variables", y = "Ganancia de informacion", 
         title = "Ganancia de información al desagrupar las variables independientes por A16") + 
  theme_bw() + 
  theme(axis.line = element_line(color = "black"), 
        axis.text = element_text(),
        panel.border = element_blank(), 
        legend.background = element_rect(size = 1), 
        legend.key = element_rect(size = 0.2), 
        panel.grid = element_blank(), 
        plot.subtitle = element_text(color = "darkgoldenrod4"))
}

# ejecucion de la funcion anterior para el data.frame "datos"
grafica_ganancia_informacion(datos)

Si graficamos A11 y A15 tomando en una escala logarítmica se observa mejor la influencia de la variable A16.

lista_graf_extremos <- map(c("A11", "A15"), 
                           ~grafica_por_tipo(datos, .x) + 
                             scale_y_log10() + 
                             theme(plot.subtitle = element_text(size = 6)))
      

grid.arrange(grobs = lista_graf_extremos, ncol = 2, 
             top = "Distribución de variables con valores extremos condicionadas a A16 en escala logarítmica")

0.7 Imputación de valores perdidos

Usamos un bosque aleatorio para la imputación de los valores perdidos, ya que las funciones para los algoritmos de clasificacion Ridge y LASSO no admiten valores perdidos

# imputacion de valores
datos_completos <- missForest(datos)$ximp 
##   missForest iteration 1 in progress...done!
##   missForest iteration 2 in progress...done!
##   missForest iteration 3 in progress...done!
##   missForest iteration 4 in progress...done!
##   missForest iteration 5 in progress...done!
##   missForest iteration 6 in progress...done!
##   missForest iteration 7 in progress...done!
##   missForest iteration 8 in progress...done!

0.8 Entreno de modelos con Rigde y Lasso

En primer lugar dividimos el data set en un training set y un test set

train <- datos_completos[1:590,]
test <- datos_completos[591:690, ]

A continuación preparamos los datos para la estimación

x_train <- data.matrix(train[, 1:15 ])
y_train <- data.matrix(train$A16)

x_test <- data.matrix(test[, 1:15])

# preparando y test en formato 0 y 1
y_test <- data.matrix(test$A16)

Construimos los modelos Ridge y lasso. Hacemos tres particiones para cada valor de lambda en el cross validation.

# establecimiento de 3 folds para el cross validation
n_folds <- 3

# fijando seed para elatoriedad
set.seed(123)


cv_ridge <- cv.glmnet(x_train, y_train, family='binomial', alpha=0, 
                      parallel=TRUE, standardize=TRUE, type.measure='auc',
                      nfolds = n_folds)

cv_lasso <- cv.glmnet(x_train, y_train, family='binomial', alpha=1, 
                      parallel=TRUE, standardize=TRUE, type.measure='auc',
                      nfolds = n_folds)

Calculamos el mejor AUC (Area Under de Curve). Tanto la estimación puntual como el error estandar del AUC estimado

# estimacion del mejor AUC en cada modelo
auc_ridge <- max(cv_ridge$cvm)
auc_lasso <- max(cv_lasso$cvm)

# Estimacion del error estandar del mejor AUC
se_auc_ridge <- cv_ridge$cvsd[which.max(cv_ridge$cvm)]/sqrt(n_folds) 
se_auc_lasso <- cv_lasso$cvsd[which.max(cv_lasso$cvm)]/sqrt(n_folds)

resultados_auc <- data.frame("AUC" = c("V. esperado", "Err. estandar"),
                             "Ridge" = round(c(auc_ridge, se_auc_ridge), 3), 
                            "LASSO" = round(c(auc_lasso, se_auc_lasso), 3)
                            )

# row.names(resultados_auc) <- c("AUC", "SE")

kable(resultados_auc) %>% kable_styling(full_width = FALSE, bootstrap_options = "condensed") %>% 
  row_spec(0, color = "white", background = "steelblue")
AUC Ridge LASSO
V. esperado 0.910 0.907
Err. estandar 0.008 0.010

Ambos modelos dan resultados muy similares. Aunque el modelo con Ridge tiene mejor AUC estimado, los modelos están a menos de un error estandar de distancia en cuanto a la estimación de esta medida. Por principio de parsimonia escogemos el modelo LASSO, que es que incluye menos variables.

Para el modelo de lasso, los resultados no cambian significativamente hasta que la penalización, \(\lambda\), alcanza valores relativamente elevados. Consecuentemente el modelo final contempla solamente cuatro varibles con un efecto significativo sobre A16.

plot(cv_lasso)

0.8.1 Predicciones sobre el test set

y_pred_lasso <- as.numeric(predict.glmnet(cv_lasso$glmnet.fit, newx=x_test, s=cv_lasso$lambda.min)>.5)
conf_matrix <- confusionMatrix(as.factor(y_pred_lasso), as.factor(y_test), mode = "everything", positive = "1")
conf_matrix
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction  0  1
##          0 83  6
##          1  3  7
##                                           
##                Accuracy : 0.9091          
##                  95% CI : (0.8344, 0.9576)
##     No Information Rate : 0.8687          
##     P-Value [Acc > NIR] : 0.1475          
##                                           
##                   Kappa : 0.5583          
##  Mcnemar's Test P-Value : 0.5050          
##                                           
##             Sensitivity : 0.53846         
##             Specificity : 0.96512         
##          Pos Pred Value : 0.70000         
##          Neg Pred Value : 0.93258         
##               Precision : 0.70000         
##                  Recall : 0.53846         
##                      F1 : 0.60870         
##              Prevalence : 0.13131         
##          Detection Rate : 0.07071         
##    Detection Prevalence : 0.10101         
##       Balanced Accuracy : 0.75179         
##                                           
##        'Positive' Class : 1               
##

0.9 Log ODDS

Los log odds son los coeficientes del modelo, porque los modelos logísticos son lineales en el logaritmo del odds

# Extraccion de los coeficientes del mejor modelo en formato data.frame
logodds <- tidy(coef(cv_lasso, s = cv_lasso$lambda.min))[c(1, 3)]
names(logodds) <- c("Variable", "Beta")

logodds$Beta <- round(logodds$Beta, 4)

# generacion de tabla
kable(logodds) %>% kable_styling(full_width = F, bootstrap_options = "condensed") %>% 
  row_spec(0, color = "white", background = "steelblue")
Variable Beta
(Intercept) -5.0663
A9 2.5684
A10 0.5743
A11 0.0130
A15 0.0000

Vemos como hay variables que no aparecen. Son las variables a las que el modelo ha asignado un coeficiente de 0.

La interpretación del los coeficientes entonces, se lee en términos del odds ratio. Por ejemplo, el coeficiente de A9 nos dice que al incrementar el valor de esta variable en una unidad, el odds ratio se multiplica por \(e^{2.5684} = 13.04\)

0.10 Rentabilidad esperada

Por cada verdadero positivo, hay una ganacia de 100 €. Por cada falso positivo hay una pérdida de 20 €.

Un verdadero positivo es la sensibilidad. 1 - falso positivo es la especificidad. Calculamos a partir de aquí, y los resultados del modelo, la rentabilidad del modelo.

sensibilidad <- round(conf_matrix$byClass["Sensitivity"], 3)
especificidad <- round(conf_matrix$byClass["Specificity"], 3)
rent_esp <- sensibilidad*100 - especificidad*20
names(rent_esp) <- NULL

\[R = FPR \cdot 100 - FPR \cdot 20 = 34.5\]

La rentabilidad esperada es de 34.5 € por caso

0.11 Conclusiones

El modelo final tiene una precisión elevada pero esto se debe en parte a la asimetría de los datos en el test set. Un 0.86 de los valores para A16 en el test set son 0. Y es ahí precisamente donde el modelo se comporta bien. El modelo tiene una especifidad de 0.965: cuando un caso es negativo, lo predice correctamente un 96.5% de las veces. En cambio, cuando un caso es positivo, solo lo precide correctamente un 53.8% de las veces.

Obviamente, la selección no aleatoria del test set es un problema. Es un problema manifiesto porque la distribución de la variable A16 en el test set no parece ser representativa de la población.Se debería reformular el problema con una selección aleatoria del training y el test set.

bind_rows(
  tibble(A16 = as.factor(y_train), Set = "Train"),
  tibble(A16 = as.factor(y_test), Set = "Test")
  ) %>%
  filter(!is.na(A16)) %>% 
  mutate(
    Set = fct_drop(Set),
    Set = fct_relevel(Set, c("Train", "Test"))
    ) %>% 
  ggplot(aes(x = Set)) + 
  geom_bar(aes(fill = A16), position = "dodge", width = 0.6) + 
  scale_fill_manual(values = c("orange", "steelblue")) + 
  labs(x = "Set", 
       y = "Frecuencia",
       title = "Disparidad en la distribución de A16 en train y test set") + 
  theme_classic() + 
  theme()

La presencia de valores perdidos tambien es un problema potencial que debería examinarse. Una resolución alternativa del problema omitiendo algunos de los valores extremos podría dar resultados diferentes.