```{r setup, message=FALSE, warning=FALSE, include=FALSE}
knitr::opts_chunk$set(
	echo = TRUE,
	fig.width = 8,
	message = FALSE,
	warning = FALSE)
```

```{r librerias}
library(tidyverse)
library(tidymodels)
library(gtsummary)
library(modelsummary)
library(gt)

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

# Introducción     
Operemos con [Tidymodels](https://www.tidymodels.org/) 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. 

```{r}
base_juguete <- readRDS(file = "data/eut_juguete_clase4.RDS")


ggplot(base_juguete,aes(x = horas_trabajo_mercado,
                         y = ingreso_familiar,
                         color = 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). 

```{r}
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
```
## 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**

---

```{r}
base_para_modelar <- base_juguete %>% 
  mutate(realiza_trabajo_domestico = factor(realiza_trabajo_domestico))

#Lo entreno 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.
```{r}
tidy(log_fiteado)
```
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.
```{r}
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)
  )
```


----

*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.    

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

predict(object = log_fiteado,new_data = data_nueva)
```
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

```{r}
augment(x = log_fiteado,new_data = data_nueva)
```
### 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 
```{r}
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
```
## 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

```{r}
matriz_confusion <- conf_mat(data = base_con_pred,
                             truth = realiza_trabajo_domestico,
                             estimate = .pred_class)

matriz_confusion
```
>¿Quien se anima a leer estos resultados?

### Métricas derivadas de la matriz
Accuracy, Sensitividad/Recall, Specificity y Precision

```{r}
accu <- accuracy(data = base_con_pred,
            truth = realiza_trabajo_domestico,
            estimate =  .pred_class)

accu
```

- El modelo predice correctamente un `scales::percent(accu,accuracy = 0.01) ` de los casos.

```{r}
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)
```
    
 - 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.
```{r}
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()
```
<br>   

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.

```{r}
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    
```{r}
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     

```{r}
#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)
```


### Aplico predicciones a la base de testeo    

```{r}
base_test_con_pred <- augment(x = modelo_sobre_train,new_data = base_test)
base_test_con_pred
```
### 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

```{r}
matriz_confusion_test<- conf_mat(data = base_test_con_pred,
                                 truth =realiza_trabajo_domestico,
                                 estimate =.pred_class)

matriz_confusion_test


accuracy(data = base_test_con_pred,
                                 truth =realiza_trabajo_domestico,
                                 estimate =.pred_class) %>%
        bind_rows(specificity(base_test_con_pred,truth = realiza_trabajo_domestico,estimate =  .pred_class)) %>%
        bind_rows(precision(base_test_con_pred,truth = realiza_trabajo_domestico,estimate =  .pred_class))
```