Regresion Logistica

Primera parte

library(dplyr)
library(ggplot2)
library(tidyr)
library(broom)
library(ISLR)
library(GGally)
library(modelr)
library(pROC)
library(cowplot)
library(OneR)
library(rlang)
library(purrr)
library(caret)
set.seed(1992)

La regresión logística es útil para problemas de predicción de clases. El problema que vamos a tratar de resolver es predecir si una persona va defaultear su deuda de tarjeta de crédito en base a ciertos predictores.

Conjunto de datos

Este conjunto de datos proviene de la librería ISLR (Introduction to Statistical Learning Using R) de James, Witten, Hastie y Tibshirani.

default <- Default
glimpse(default)

Observations: 10,000
Variables: 4
$ default <fct> No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No, No,...
$ student <fct> No, Yes, No, No, No, Yes, No, Yes, No, No, Yes, Yes, No, No, No, No, No, Yes, No, No, No, No, No, No, No, No, No, No, No, Yes, ...
$ balance <dbl> 729.5265, 817.1804, 1073.5492, 529.2506, 785.6559, 919.5885, 825.5133, 808.6675, 1161.0579, 0.0000, 0.0000, 1220.5838, 237.0451...
$ income  <dbl> 44361.625, 12106.135, 31767.139, 35704.494, 38463.496, 7491.559, 24905.227, 17600.451, 37468.529, 29275.268, 21871.073, 13268.5...

Tiene 4 variables:

• default: La clase que queremos predecir
• student: Binaria que indica si la persona es estudiante
• balance: Balance promedio que le queda a la persona luego de sus pagos mensuales
• income: Ingreso de la persona

Exploratorias

Analicemos la distribución de la clase

default %>% group_by(default) %>% summarise(numero_casos=n())

Vemos que estamos trabajando con un problema de clasificación con un claro desbalance de clase.

Realizamos un gráfico exploratorio completo para ver el comportamiento y las relaciones entre las variables. El color rojo designa a quienes no defaultean y el azul a los que sí.

ggpairs(default,mapping = aes(colour= default)) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + theme_bw()

¿Qué pueden decir de la relación entre balance y default?

¿Y entre income y balance?

¿Cuáles parecen ser buenas variables para predecir la probabilidad de default de una persona?

Limpieza

Para modelizar va a ser más necesario tener la variable default como numérica. Definimos la variable como {0,1} para los valores {“No”,“Yes”}

default <- default%>% mutate(default= case_when(default=="No"~0,
                                                default=="Yes"~1))

Problema

Queremos estimar \(P(Default=Yes|X)=P(X)\) para cada individuo y partir de ello poder definir un punto de corte para predecir quienes son los que van a entrar en default.

Regresión lineal

En este caso estamos modelando la probabilidad de la siguiente manera:

\(P(X)= \beta_0 + \sum\limits_{j=1}^p \beta_j X\)

Veamos que tan bueno es el modelo lineal para esto, usando balance como predictor.

test_mco <- default %>% 
              lm(formula = default~balance, data = .) 

tdy <- test_mco %>% tidy()
tdy

test_mco %>% glance()

Ambos estimadores son significativos y el test de significatividad global del modelo también es significativo. Veamos un gráfico de nuestro modelo

Parece tener bastantes problemas para estimar la probabilidad de default de los individuos. Por ejemplo, vemos que hay varios individuos a los cuales les asigna una probabilidad negativa.

Regresión logística

Para evitar estos problemas, usamos la funcion logistica

\(P(X)= \frac{e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X}}{1+e^{\beta_0 + \sum\limits_{j=1}^p \beta_j X}}\)

El lado derecho se llama expit

Esta funcion acota el resultado entre 0 y 1, lo cual es mucho mas adecuado para modelar una probabilidad.

Luego de hacer algunas operaciones, podemos llegar a la expresion:

\(\log {\frac{P(x)}{1-P(x)}}= \beta_0 + \sum\limits_{j=1}^p \beta_j X\)

El lado izquierdo es el logaritmo de los odds y se llama logit

Partición Train y Testing

Realizamos una partición entre dataset de entrenamiento (70%) y testeo (30%) usando la función resample_partition del paquete modelr

train_test <- default %>% resample_partition(c(train=0.7,test=0.3))
default <- train_test$train %>% as_tibble()
test <- train_test$test %>% as_tibble()

Creación de fórmulas

Para aplicar la regresion logistica primero usamos la funcion formulas del paquete modelr para crear un objeto que contiene todas las formulas que vamos a utilizar.

logit_formulas <- formulas(.response = ~default, # único lado derecho de las formulas.
                         bal= ~balance, 
                         stud= ~student,  
                         inc= ~income,  
                         bal_stud=~balance+student, 
                         bal_inc=~balance+income, 
                         stud_inc=~student+income,  
                         full= ~balance + income + student  
                         )

Creación de modelos

Procedemos a crear los modelos a partir de estas fórmulas

models <- data_frame(logit_formulas) %>% # dataframe a partir del objeto formulas
  mutate(models = names(logit_formulas), # columna con los nombres de las formulas
         expression = paste(logit_formulas), # columna con las expresiones de las formulas
         mod = map(logit_formulas, ~glm(.,family = 'binomial', data = default))) # Que estamos haciendo acá? Que vamos a encontrar en la columna?

Modelos simples

Probamos los primeros tres modelos, aquellos que tienen un único predictor. Usamos la función tidy para obtener los parámetros estimados para estos tres modelos.

models %>% 
  filter(models %in% c('bal','stud','inc')) %>%
  mutate(tidy = map(mod,tidy)) %>%  # Qué realizamos en este paso? Que va a tener esta columna?
  unnest(tidy, .drop = TRUE) %>% 
  mutate(estimate=round(estimate,5),
         p.value=round(p.value,4))

¿Son significativos?

¿Qué interpretación pueden darle a estos valores?

Modelo completo

Ahora probamos con un modelo que utiliza las tres predictoras

models %>% 
  filter(models == "full") %>%
  mutate(tidy = map(mod,tidy)) %>%
  unnest(tidy, .drop = TRUE) %>% 
  mutate(estimate=round(estimate,5),
         p.value=round(p.value,4))

¿Que cambios hay respecto a los 3 modelos individuales previos?

Evaluación de todos los modelos

Con map() agregamos la función glance para traernos información relevante para el diagnóstico del modelo.

Con unnest() podemos ver la evaluación de cada modelo. Por último ordenamos las modelos por el deviance.

models <- models %>% 
  mutate(glance = map(mod,glance))
models %>% 
  unnest(glance, .drop = TRUE) %>%
  mutate(perc_explained_dev = 1-deviance/null.deviance) %>% 
  select(-c(models, df.null, AIC, BIC)) %>% 
  arrange(deviance)

El modelo que utiliza las 3 variables es el que minimiza el deviance. Los 3 últimos modelos reducen muy poco el deviance respecto a la deviance nula.

Gráficos de evaluación

Realizamos los gráficos para el modelo completo y uno de los modelos con mayor deviance (student+income).

En este caso estamos:

1. Agregando las predicciones con augment con el parámetro type=“response” ¿Por qué hacemos esto? ¿Cuál es el valor por default de este parámetro?
2. Armando las curvas ROC con pROC

models <- models %>% 
  mutate(pred= map(mod,augment, type.predict = "response"))
models$pred[1]

$bal

NA

prediction_full <- models %>% 
  filter(models=="full") %>% 
  unnest(pred, .drop=TRUE)
roc_full <- roc(response=prediction_full$default, predictor=prediction_full$.fitted)
prediction_bad <- models %>% 
  filter(models=="stud_inc") %>% 
  unnest(pred, .drop=TRUE)
roc_bad <- roc(response=prediction_bad$default, predictor=prediction_bad$.fitted)

Violin plots

violin_full=ggplot(prediction_full, aes(x=default, y=.fitted, group=default,fill=factor(default))) + 
  geom_violin() +
  theme_bw() +
  guides(fill=FALSE) +
  labs(title='Violin plot', subtitle='Modelo completo', y='Predicted probability')
violin_bad=ggplot(prediction_bad, aes(x=default, y=.fitted, group=default, fill=factor(default))) + 
  geom_violin() + 
  theme_bw() +
  guides(fill=FALSE) +
  labs(title='Violin plot', subtitle='Modelo malo', y='Predicted probability')
plot_grid(violin_bad, violin_full)

¿Qué es lo que estamos viendo en ellos? (Especial atención al eje de ordenadas)

¿Cuál parece ser un punto de corte adecuado para cada modelo?

Curvas ROC

ggroc(list(full=roc_full, bad=roc_bad), size=1) + geom_abline(slope = 1, intercept = 1, linetype='dashed') + theme_bw() + labs(title='Curvas ROC', color='Modelo')

print(paste('AUC: Modelo completo', roc_full$auc))

[1] "AUC: Modelo completo 0.949391254006722"

print(paste('AUC: Modelo malo', roc_bad$auc))

[1] "AUC: Modelo malo 0.565799201537779"

¿Qué significa cada uno de los ejes?

Gráfico de Hosmer-Lemeshow

Hosmer_Lemeshow_plot <- function(dataset, predicted_column, class_column, bins, positive_value, color='forestgreen', nudge_x=0, nudge_y=0.05){
  "Realiza un grafico de Hosmer-Lemeshow para un dataset"
  
  "* dataset: conjunto de datos
   * predicted_column: columna con la probabilidad predicha
   * class_column: columna con la clase a predecir
   * possitive_value: valor de la clase a predecir
   * bins: cantidad de grupos del gráfico
   * color: color de los puntos
   * nudge_x: desplazamiento de la etiqueta en el eje x
   * nudge_y: desplazamiento de la etiqueta en el eje y"
  
  # Asignar los grupos a las observaciones de acuerdo a la probabilidad predicha
  dataset['group'] <- bin(dataset[predicted_column], nbins = bins, method = 'l', labels=c(1:bins))
  
  # Contar la cantidad de casos positivos por grupo
  positive_class <- dataset %>% filter(!!sym(class_column)==positive_value) %>% group_by(group) %>% count()
  
  # Obtener la media de las predicciones por grupo
  HL_df <- dataset %>% group_by(group) %>% summarise(pred=mean(!!sym(predicted_column)), count=n()) %>%
            inner_join(.,positive_class) %>%
            mutate(freq=n/count)
  
  # Gráfico 
  HM_plot <- ggplot(HL_df, aes(x=pred, y=freq)) + geom_point(aes(size=n), color=color) +
                geom_text(aes(label=n),nudge_y = nudge_y)+
                geom_abline(slope = 1, intercept = 0, linetype='dashed') + 
                theme_bw() +
                labs(title='Hosmer-Lemeshow', size='Casos', x="Probabilidad Predicha", y="Frecuencia observada")
  return(HM_plot)
}
Hosmer_Lemeshow_plot(prediction_full, '.fitted', 'default', 10, 1) + labs(subtitle="Modelo completo")

Hosmer_Lemeshow_plot(prediction_bad, '.fitted', 'default', 10, 1, color = "firebrick", nudge_y = 0.003) + scale_x_continuous(limits = c(0.02,.06)) + scale_y_continuous(limits = c(.02,.06)) + labs(subtitle="Modelo malo")

¿Qué vemos en estos gráficos?

¿Para qué valores parece existir una sobreestimación de la probabilidad? ¿Para cuáles subestimación?

Punto de corte

Hasta ahora hemos evaluado el modelo de manera general, pero el resultado final del modelo debe consistir en asignar al individuo una clase predicha. En nuestro caso debemos establecer un punto de corte según el cual vamos a separar a los indivuos en quienes defaultean y quienes no.

Probamos varios puntos de corte y graficamos el accuracy, la sensibilidad, la especificidad, el recall y la precision para cada uno de ellos.

Clases predichas / Clases	Positiva	Negativa
Positiva	True Pos	False Pos
Negativa	False Neg	True Neg

Recordemos que:

\(sensitivity = recall = \frac{TP}{TP+FN}\)

\(specificity = \frac{TN}{TN+FP}\)

\(precision = \frac{TP}{TP+FP}\)

prediction_metrics <- function(cutoff, predictions=prediction_full){
  table <- predictions %>% 
    mutate(predicted_class=if_else(.fitted>cutoff, 1, 0) %>% as.factor(),
           default= factor(default))
  
  confusionMatrix(table(table$predicted_class, table$default), positive = "1") %>%
    tidy() %>%
    select(term, estimate) %>%
    filter(term %in% c('accuracy', 'sensitivity', 'specificity', 'precision','recall')) %>%
    mutate(cutoff=cutoff)
  
}
cutoffs = seq(0.01,0.95,0.01)
logit_pred= map_dfr(cutoffs, prediction_metrics)%>% mutate(term=as.factor(term))
ggplot(logit_pred, aes(cutoff,estimate, group=term, color=term)) + geom_line(size=1) +
  theme_bw() +
  labs(title= 'Accuracy, Sensitivity, Specificity, Recall y Precision', subtitle= 'Modelo completo', color="")

¿Qué podemos observar en el gráfico?

¿Podemos definir un buen punto de corte? ¿Cuál sería?

¿Por qué la especificidad tiene ese comportamiento?

Dataset de testing

Seleccionamos el modelo completo, ya que es el que maximizaba el porcentaje de deviance explicada y en base a lo que vimos definimos un punto de corte en 0.25 (pueden probar otros)

sel_cutoff = 0.25
# Creamos el modelo
full_model <- glm(logit_formulas$full, family = 'binomial', data = default)
# Agregamos la predicciones al dataset de testeo
table= augment(x=full_model, newdata=test, type.predict='response') 
# Clasificamos utilizamos el punto de corte
table=table %>% mutate(predicted_class=if_else(.fitted>0.25, 1, 0) %>% as.factor(),
           default= factor(default))
# Creamos la matriz de confusión
confusionMatrix(table(table$default, table$predicted_class), positive = "1")

Confusion Matrix and Statistics

   
       0    1
  0 2839   65
  1   44   53
                                          
               Accuracy : 0.9637          
                 95% CI : (0.9564, 0.9701)
    No Information Rate : 0.9607          
    P-Value [Acc > NIR] : 0.21383         
                                          
                  Kappa : 0.4744          
 Mcnemar's Test P-Value : 0.05541         
                                          
            Sensitivity : 0.44915         
            Specificity : 0.98474         
         Pos Pred Value : 0.54639         
         Neg Pred Value : 0.97762         
             Prevalence : 0.03932         
         Detection Rate : 0.01766         
   Detection Prevalence : 0.03232         
      Balanced Accuracy : 0.71695         
                                          
       'Positive' Class : 1               
                                          

Segunda parte

Desbalanceo de la clase

Al explorar el dataset vimos que existía un fuerte desbalance de clase. Sòlo el 3% de las observaciones pertenecen a personas que defaultearon. Esto puede tener un efecto en las estimaciones del modelo y su clasificación final.

Existen dos maneras sencillas con las cuales podemos trabajar con una clase desbalanceada:

• Sobre-sampleo (oversampling) de la clase minoritaria
• Sub-sampleo (undersampling) de la clase mayoritaria

La función glm puede tomar como argumento una columna (weigths) de ponderadores para poder hacer esto. Podemos asignar pesos mayores a 1 a la clase minoritaria (oversampling) o menores a 1 a la clase mayoritaria (undersampling). En nuestro problema vamos a realizar un sobresampleo de la clase minoritaria.

# Creamos la columna de ponderadores
default <- default %>% mutate(wt= if_else(default==1,20,1))
# Creamos los modelos con la data 'balanceada'
balanced_models <- data_frame(logit_formulas) %>% # dataframe a partir del objeto formulas
  mutate(models = names(logit_formulas), # columna con los nombres de las formulas
         expression = paste(logit_formulas), # columna con las expresiones de las formulas
         mod = map(logit_formulas, ~glm(.,family = 'binomial', data = default, weights = wt))) #Pasamos la columna wt como ponderadores

Vemos las estimaciones de los parametros para el modelo completo. ¿Existen cambios?

Ahora veamos la evaluación de los modelos ¿Qué pasó con el porcentaje de deviance explicada? ¿Y con la nula?

Violin plots, Curvas ROC y AUCs

Realizamos los gráficos de violin, las curvas ROC y calculamos los AUC

[1] "AUC: Modelo completo 0.949423207532511"
[1] "AUC: Modelo malo 0.56572151061233"

¿Dónde se ven los cambios más notorios respecto a nuestros modelos anteriores que no tenían en cuenta el desbalance de la clase?

Punto de corte

Volvemos a realizar las pruebas para varios puntos de corte y graficamos el accuracy, la sensibilidad, la especificidad, el recall y la precision para cada uno de ellos.

¿Qué cambios vemos respecto al gráfico anterior?

Dataset de testing

Probamos en el dataset de testing nuestro modelo balanceado. No es necesario que le creemos pesos al dataset de testeo.

Confusion Matrix and Statistics

   
       0    1
  0 2287  617
  1    5   92
                                          
               Accuracy : 0.7927          
                 95% CI : (0.7778, 0.8071)
    No Information Rate : 0.7637          
    P-Value [Acc > NIR] : 8.196e-05       
                                          
                  Kappa : 0.1818          
 Mcnemar's Test P-Value : < 2.2e-16       
                                          
            Sensitivity : 0.12976         
            Specificity : 0.99782         
         Pos Pred Value : 0.94845         
         Neg Pred Value : 0.78753         
             Prevalence : 0.23625         
         Detection Rate : 0.03066         
   Detection Prevalence : 0.03232         
      Balanced Accuracy : 0.56379         
                                          
       'Positive' Class : 1