library(tidyverse)
library(tidymodels)
library(gtsummary)
library(modelsummary)
library(gt)

#options(dplyr.summarise.inform = FALSE)
options(scipen = 999)

Introducción

Operemos con Tidymodels para hacer un modelo de clasificación con el método de la regresión logística. Usaremos nuestra base de juguete de 200 datos ficticios, con la variable objetivo realiza_trabajo_domestico. Por ahora, tomaremos únicamente como variables explicativas las horas trabajadas en el mercado y el ingreso total familiar.

Visualicemos previamente estos datos.

base_juguete <- readRDS(file = "data/eut_juguete_clase4.RDS")


ggplot(base_juguete,aes(x = horas_trabajo_mercado,
                         y = ingreso_familiar,
                         shape = realiza_trabajo_domestico))+
  geom_point(size = 3)

Ahora veamos los pasos necesarios para correr el modelo. Luego, discutamos la interpretación de los coeficientes.

Seteo el tipo de modelo

En formato tidy, comienzo especificando el tipo de modelo. En este caso, logistic_reg(). De ser necesario puedo especificar con set_mode() el modo del modelo (qué tipo de predicción tiene el outcome, numérica o categórica) Luego, especifico cual es el sistema que utilizaré para estimar el modelo con la función set_engine() (en muchos casos responde al paquete que voy a usar para correr el modelo).

log_model <- logistic_reg() %>% #Defino el tipo de modelo
  set_mode("classification") %>%  #el modo (regresión o clasificación)
  set_engine("glm") #el motor sera de "Generalized linear models"
log_model
## Logistic Regression Model Specification (classification)
## 
## Computational engine: glm

Fiteo (entreno) el modelo

Tomo la especificación anterior y explicito qué variables usaré y de que dataset.

Importante: La variable objetivo realiza_trabajo_domestico debo convertirla en factor para correr una regresión logistica (es decir, que sea una variable con categorías). Por otra parte, si van a utilizar funciones que calculan métricas derivadas de la matriz de confusión es mucho muy importante que las categorías estén ordenadas con el valor positivo como primer nivel y el negativo como el segundo nivel

base_para_modelar <- base_juguete %>% 
  mutate(realiza_trabajo_domestico = factor(realiza_trabajo_domestico))

#Lo entreo con los datos 
log_fiteado <- log_model %>% 
  fit(realiza_trabajo_domestico ~ horas_trabajo_mercado+ingreso_familiar,
      data = base_para_modelar)

Veamos resultados de los parámetros estimados

Si bien podemos correr directo el objeto log_fiteado para ver los resultados, la función tidy() nos transforma los coeficientes estimados (y también su desvío estandar, z-statistic y p.valor) a un formato data_frame. Para visualizarlo de forma más prolija usamos aquí también la función gt() del paquete homónimo.

tidy(log_fiteado)
## # A tibble: 3 × 5
##   term                    estimate std.error statistic  p.value
##   <chr>                      <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)            7.01      1.12           6.24 4.47e-10
## 2 horas_trabajo_mercado -0.130     0.0224        -5.84 5.26e- 9
## 3 ingreso_familiar      -0.0000995 0.0000194     -5.13 2.86e- 7

El paquete modelsummary y su función homónima nos arroja también datos sobre la bondad de ajuste del modelo. Tenemos medidas como el Akaike information criterion (AIC) y el Bayesian Information Criterion (BIC). Una herramienta interesante es que podemos pasar una lista de modelos y comparar entre ellos. Solo para estos fines, entrenemos otro modelo que también incluye la variable de menores en el hogar.

log_fiteado2 <- log_model %>% 
  fit(realiza_trabajo_domestico ~ horas_trabajo_mercado+ingreso_familiar+menores_hogar,
      data = base_para_modelar)


modelsummary::modelsummary(
  list("logistica 2 predictores " = log_fiteado,
       "logistica 3 predictores" = log_fiteado2)
  )
logistica 2 predictores logistica 3 predictores
(Intercept) 7.006 5.406
(1.123) (1.225)
horas_trabajo_mercado -0.130 -0.131
(0.022) (0.026)
ingreso_familiar -0.000 -0.000
(0.000) (0.000)
menores_hogarSi 2.684
(0.553)
Num.Obs. 200 200
AIC 131.4 104.1
BIC 141.3 117.3
RMSE 0.31 0.28

Aclaración: Dado que nos interesan más otro tipo de métricas, no veremos en detalle la estimación de medidas de bondad de ajuste como AIC y BIC. A los fines prácticos basta con tener noción de que:

  • Su utilidad es para comparar modelos (la interpretación del AIC o BIC de 1 modelo en si mismo es irrelevante)
  • Las formulas de AIC y BIC tienen la siguiente “pinta”: \(xIC= complejidad-ajuste\).
    Un termino positivo que actua como penalización a la complejidad del modelo (asociado a la cantidad de parametros o variables que utiliza). Un termino restando asociado a la capacidad de ajuste del modelo (vinculada a la máxima verosimilitud, es decir, cuan cerca le paso a los valores reales con mi predicción de probabilidad)
  • Al comparar dos modelos, son mejores los que tienen menor valor de AIC o BIC.
  • En el ejemplo anterior, a pesar de haber agregado complejidad en el 2do modelo (1 variable más), ganamos mucho en ajuste, por eso el AIC y BIC son más bajos.

Predicción

Aplicado sobre datos nuevos

Si queremos sólo obtener la predicción que nuestro modelo hará sobre un nuevo par de datos, la función predict() requiere como imput un data.frame que contenga columnas nombradas igual que los predictores de nuestro modelo y devuelve, en el mismo orden, las predicciones para la variable Y.

data_nueva<- data.frame(horas_trabajo_mercado = c(40,13),
                        ingreso_familiar = c(60000,10000))

predict(object = log_fiteado,new_data = data_nueva)
## # A tibble: 2 × 1
##   .pred_class
##   <fct>      
## 1 No         
## 2 Si

Una función que puede resultar todavía más úlil que predict() es augment(). Al pasarle un modelo entrenado y un dataset para realizar predicciones, esta última añade a dicho dataset el valor de la clase predicha y también las probabilidades estimadas que respaldan dicha predicción

augment(x = log_fiteado,new_data = data_nueva)
## # A tibble: 2 × 5
##   .pred_class .pred_No .pred_Si horas_trabajo_mercado ingreso_familiar
##   <fct>          <dbl>    <dbl>                 <dbl>            <dbl>
## 1 No            0.985    0.0150                    40            60000
## 2 Si            0.0132   0.987                     13            10000

Aplicado sobre el total de la base

Tomemos la función augment() para agregar a nuestro dataset base_para_modelar las probabilidades asignadas por nuestro modelo a cada uno de los calsos, y por ende, la clase predicha

base_con_pred<- augment(x = log_fiteado, base_para_modelar) %>% 
  mutate(prediccion_gw = ifelse(.pred_No > 0.8,yes = "No",no = "Si"))

base_con_pred
## # A tibble: 200 × 8
##    .pred_class .pred_No .pred_Si realiza_trabajo_domestico horas_trabajo_mercado
##    <fct>          <dbl>    <dbl> <fct>                                     <int>
##  1 No             0.805   0.195  No                                           56
##  2 Si             0.405   0.595  No                                           33
##  3 No             0.985   0.0147 No                                           36
##  4 No             0.965   0.0348 No                                           70
##  5 Si             0.288   0.712  No                                           43
##  6 Si             0.341   0.659  No                                           30
##  7 No             0.957   0.0431 No                                           39
##  8 No             0.858   0.142  No                                           51
##  9 No             0.937   0.0632 No                                           64
## 10 No             0.829   0.171  No                                           52
## # ℹ 190 more rows
## # ℹ 3 more variables: ingreso_familiar <dbl>, menores_hogar <chr>,
## #   prediccion_gw <chr>

Metricas de evaluación

Matriz de confusión

Sobre la base de los valores reales y los valores predichos podríamos tomar la matriz de confusión como un tabulado bivariado entre ambas columas

matriz_confusion <- conf_mat(data = base_con_pred,
                             truth = realiza_trabajo_domestico,
                             estimate = .pred_class)

matriz_confusion
##           Truth
## Prediction No Si
##         No 86 12
##         Si 14 88

¿Quien se anima a leer estos resultados?

Métricas derivadas de la matriz

Accuracy, Sensitividad/Recall, Specificity y Precision

accu <- accuracy(data = base_con_pred,
            truth = realiza_trabajo_domestico,
            estimate =  .pred_class)

accu
## # A tibble: 1 × 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary          0.87
  • El modelo predice correctamente un scales::percent(accu,accuracy = 0.01) de los casos.
sens <- sensitivity(base_con_pred,truth = realiza_trabajo_domestico,estimate =  .pred_class)

spec <- specificity(base_con_pred,truth = realiza_trabajo_domestico,estimate =  .pred_class)

prec <- precision(base_con_pred,truth = realiza_trabajo_domestico,estimate =  .pred_class)

bind_rows(accu,sens,spec,prec)
## # A tibble: 4 × 3
##   .metric     .estimator .estimate
##   <chr>       <chr>          <dbl>
## 1 accuracy    binary         0.87 
## 2 sensitivity binary         0.86 
## 3 specificity binary         0.88 
## 4 precision   binary         0.878
  • De los casos positivos (que efectivamente realizan trabajo doméstico), el modelo predice bien un 87%.
  • De los casos negativos (que no realizan trabajo doméstico), el modelo predice bien un 88%.
  • De los casos que son clasificados como positivo 87.7%. lo son.

Veamos esto de forma gráfica para tratar de entender como está operando nuestro modelo:
En este gráfico utilizamos el parametro shape ara distinguir el valor real de la variable objetivo (“No” o “Si”) y el color para mostrar las probabildiades de pertenecer a la clase “Si” estimadas por nuestro modelo.

ggplot(base_con_pred,aes(x = horas_trabajo_mercado,
                         y = ingreso_familiar,
                         shape = realiza_trabajo_domestico,
                         color = .pred_Si))+
  geom_point(size = 3)+
  scale_color_viridis_c()


En este otro utilizamos el color no para mostrar las probabildiades, sino para ver la clase predicha por nuestro modelo (lo que tenia probabilidad mayor a 0,5 queda como “Si”, los menores a 0,5 como “No”). Ver como la “decision boundary” parece ser líneal.

ggplot(base_con_pred,aes(x = horas_trabajo_mercado,
                         y = ingreso_familiar,
                         shape = realiza_trabajo_domestico,
                         color = .pred_class))+
  geom_point(size = 3)+
  scale_color_viridis_d()+
  theme_minimal()

El ejercicio anterior tiene un problema. Utilizamos para evaluar la performance del modelo los mismos datos que utilizamos para entrenarlo. Eso es potencialmente peligroso por posible OVERFITTING.

Train-test split

Hagamos un modelo nuevo con un split de la base original en test y train. Aclaración: Si queremos resultados replicables debemos setear una semilla con set.seed(), para guardar el mecanismo de pseudo-aleatorización que realiza la computadora.

Split train-test

set.seed(18/12/2022)
base_spliteada <-  initial_split(base_para_modelar,prop = 0.8)
base_train <-  training(base_spliteada)
base_test <-  testing(base_spliteada)

Fiteo

#log_model <- logistic_reg() %>% # Ya lo corri arriba este paso 
#  set_engine("glm")

modelo_sobre_train <- log_model  %>% 
  fit(realiza_trabajo_domestico ~ horas_trabajo_mercado+ingreso_familiar,
      data = base_train)

tidy(modelo_sobre_train)
## # A tibble: 3 × 5
##   term                    estimate std.error statistic      p.value
##   <chr>                      <dbl>     <dbl>     <dbl>        <dbl>
## 1 (Intercept)            7.58      1.39           5.43 0.0000000553
## 2 horas_trabajo_mercado -0.139     0.0276        -5.04 0.000000463 
## 3 ingreso_familiar      -0.0000935 0.0000197     -4.76 0.00000196

Aplico predicciones a la base de testeo

base_test_con_pred <- augment(x = modelo_sobre_train,new_data = base_test)
base_test_con_pred
## # A tibble: 40 × 7
##    .pred_class .pred_No .pred_Si realiza_trabajo_domestico horas_trabajo_mercado
##    <fct>          <dbl>    <dbl> <fct>                                     <int>
##  1 No             0.972  0.0277  No                                           36
##  2 No             0.964  0.0357  No                                           70
##  3 Si             0.244  0.756   No                                           43
##  4 Si             0.248  0.752   No                                           30
##  5 No             0.825  0.175   No                                           51
##  6 No             0.931  0.0686  No                                           64
##  7 Si             0.209  0.791   No                                           34
##  8 No             0.990  0.00994 No                                           66
##  9 No             0.978  0.0218  No                                           45
## 10 Si             0.268  0.732   No                                           38
## # ℹ 30 more rows
## # ℹ 2 more variables: ingreso_familiar <dbl>, menores_hogar <chr>

Matriz de confusión

Sobre la base de los valores reales y los valores predichos podríamos tomar la matriz de confusión como un tabulado bivariado entre ambas columas

matriz_confusion_test<- conf_mat(data = base_test_con_pred,
                                 truth =realiza_trabajo_domestico,
                                 estimate =.pred_class)

matriz_confusion_test
##           Truth
## Prediction No Si
##         No 22  1
##         Si  6 11