3. Estimar un modelo

Agenda

Parte I

  • Especificar un modelo
  • Estimar un modelo
  • Ver los resultados del modelo

Parte II

  • Dividir los datos
  • Especificar un modelo
  • Estimar un modelo
  • Evaluar el modelo
  • Estimar un modelo mejor

Estimar un modelo

Hay muchas formas de estimar modelos en R

  • lm(formula, data, ...)
  • stan_glm(formula, data, family= "gaussian",...)
  • glmnet(x=matrix, y=vector, family="gaussian",...)

Tip

lm(): Ajusta un modelo de regresión lineal a los datos.

stan_glm(): Ajusta un modelo de regresión utilizando el paquete Stan.

glmnet(): Ajusta un modelo de regresión utilizando la regularización Lasso o Ridge.

Especificar el modelo en tidymodels

Tidymodels ofrece una sintaxis general para estimar los modelos

https://www.tidymodels.org/

  • Especificar el tipo de modelo
    • linear_reg(), logistic_reg(), decition_tree()
  • Especificar el tipo de outcome
    • Regresión para outcomes continuos
    • Clasificación: multinomial, ordinal, binaria


library(tidymodels)

tidymodels_prefer()

linear_reg() |> 
  set_engine("lm") 

linear_reg() |> 
    set_engine("glmnet") 
    
linear_reg() |>
    set_engine("stan") 

Especificar el modelo en tidymodels

Tip

library(): Carga el paquete tidymodels en R.

tidymodels_prefer(): Establece el paquete tidymodels como el preferido para las funciones de modelado en R.

linear_reg(): Crea una especificación de modelo de regresión lineal en el marco de trabajo tidymodels.

set_engine(): Establece el motor de cálculo para un objeto de especificación de modelo. Se establece en “lm” para la primera llamada a linear_reg(), “glmnet” para la segunda llamada y “stan” para la tercera.

Estimar un modelo

Estimar un modelo de regresión logistica para evaluar la asociación entre el consumo de alcohol y los factores de riesgo

mini_datos |>
    select(PYALC,factores_de_riesgo) |>
    table1::table1( ~ ., data=_)


mini_datos |>
    select(PYALC,factores_de_riesgo) |>
    filter(!is.na(PYALC)) |> 
    table1::table1( ~ . | PYALC, data=_)

Explorar los datos visualmente

mini_datos |> 
  group_by(PYALC) |> 
  summarise(across("CRPAD":"SPRPI", \(x) mean(x, na.rm = TRUE))) |> 
  pivot_longer(-PYALC) |> 
  filter(!is.na(PYALC)) |> 
  ggplot(aes(value, fct_reorder(name, value), fill = PYALC)) +
  geom_col(alpha = 0.8, position = "dodge") +
  scale_x_continuous() +
  labs(x = "Promedio", y = NULL, fill = NULL)

1. Especificar el modelo

lm_model <-
    logistic_reg() %>% 
    set_engine("glm", family="binomial") |>
    set_mode("classification")


Especificar un modelo de regresión logística. Se establece el “paquete”, en este caso “glm” y opciones del paquete (la distribución de la variable de respuesta). Además, se establece el modo de modelado en “classification” para indicar que es una variable categórica. El resultado se almacena en el objeto lm_model.

2. Estimar el modelo

lm_results <- lm_model |>
    fit(PYALC ~ ., data = mini_datos)


Aquí se ajusta-estima el modelo utilizando la función fit(). Se especifica la variable de respuesta (PYALC) y los predictores (.) y se utiliza el objeto mini_datos como base de datos.

3. Ver los resultados del modelo

lm_results |> tidy() 

lm_results |> glance()


tidy(): Esta función se utiliza para convertir los resultados de un modelo en un marco de datos “ordenado” que muestra los coeficientes del modelo, los errores estándar, los valores t y los valores p para cada variable independiente.

glance(): Esta función se utiliza para resumir los resultados de un modelo como estadísticas globales del modelo, por ejemplo, R-cuadrado ajustado, el AIC y el BIC.


Tip

lm_results |> tidy(exp=TRUE, conf.int=TRUE)

tidy(exp=TRUE, conf.int=TRUE): Esta función se utiliza para convertir los resultados de un modelo en un marco de datos “ordenado” que muestra los coeficientes del modelo, los errores estándar, los valores t y los valores p para cada variable independiente. Los argumentos exp y conf.int se utilizan para incluir los intervalos de confianza y los exponentes en los resultados

Juego 2

  1. Con select, creen una nueva base de datos con los predictores que ustedes consideren relevantes.
  2. Estimen un modelo con esta nueva base de datos
  3. Calculen los odds ratio (glance())

Parte II

Para qué necesitamos los datos?

Necesitamos:

- Estimar parametros
- Seleccionar modelos
- Sintonizar los modelos (tunning)
- Evaluar los modelos

¿Cómo gastarnos los datos de una forma que sea eficiente para todos estos pasos? (validación empírica)

Primera idea: dividir los datos en dos.

Dividir los datos

  • Dividir los datos en dos conjuntos
    • Entrenamiento
      • La mayoría de los datos (¿70%, 80%?)
      • Aquí se ajusta el modelo
    • Prueba
      • Un pequeño conjunto de datos
      • Aquí se evaluará el modelo final

Los datos de prueba se utilizan una sola vez, si se utilizan más de una vez, se convierten en parte del proceso de entrenamiento.

La división de los datos se hace al nivel de unidad independiente de observación.

Contaminación del los datos de prueba (information leakage)

Dividir los datos

Crear dos bases: 80% y 20%

set.seed(1234)

mini_datos <- mini_datos |> drop_na()

datos_divididos <- initial_split(mini_datos, prop = 0.8, strata = "PYALC")

datos_entrenamiento <- training(datos_divididos)
datos_prueba <- testing(datos_divididos)
  • La opción strata es para que los datos de entrenamiento y prueba tengan la misma distribución de la variable PYALC

A veces la selección aleatoria de la muestra es problemática, por ejemplo cuando hay una componente de tiempo en los datos. En este caso, se puede usar la función initial_time_split().

Dividir los datos

set.seed(): se utiliza para establecer una semilla para la generación de números aleatorios en R. En este caso, se utiliza para establecer la semilla en 1234, lo que garantiza que los resultados sean reproducibles. Esto es especialmente importante cuando se trabaja con modelos de aprendizaje automático, ya que los resultados pueden variar según la semilla utilizada para la generación de números aleatorios.

drop_na(): Esta función se utiliza para eliminar filas con valores faltantes de un data frame o tibble. En este caso, se utiliza para eliminar filas con valores faltantes del objeto “mini_datos”.

initial_split(): Esta función se utiliza para dividir un data frame o tibble en conjuntos de entrenamiento y prueba. En este caso, se utiliza para dividir el objeto “mini_datos” en conjuntos de entrenamiento y prueba en una proporción de 80/20, y estratificando por la columna “PYALC”.

training(): Esta función se utiliza para extraer el conjunto de entrenamiento de un objeto creado con la función initial_split(). En este caso, se utiliza para extraer el conjunto de entrenamiento del objeto “datos_divididos”.

testing(): Esta función se utiliza para extraer el conjunto de prueba de un objeto creado con la función initial_split(). En este caso, se utiliza para extraer el conjunto de prueba del objeto “datos_divididos”.

Dividir los datos

datos_entrenamiento |>
    select(PYALC,factores_de_riesgo) |>
    table1::table1( ~ ., data=_)
datos_entrenamiento |>
  select(PYALC, factores_de_riesgo) |>
  filter(!is.na(PYALC)) |>
  table1::table1(~ . | PYALC, data = _)

1. Especificar el modelo

lm_model <-
  logistic_reg() %>%
  set_engine("glm", family = "binomial") |>
  set_mode("classification")

2. Estimar el modelo

lm_results <- lm_model |>
    fit(PYALC ~ ., data = datos_entrenamiento)

3. Ver los resultados del modelo

lm_results |> tidy()
lm_results |> tidy(exp=TRUE, conf.int=TRUE)

4. Evaluar el modelo

predecir_estos_valores <- datos_entrenamiento %>% 
    select(-PYALC) %>% 
    slice(1:10)

predict(lm_results, predecir_estos_valores)
predict(lm_results, predecir_estos_valores, type="prob")
entrenamiento <- datos_entrenamiento %>% 
    select(-PYALC)

prediccion <- predict(lm_results, entrenamiento)

resultados_prueba <- cbind(prediccion, datos_entrenamiento) |> 
  select(.pred_class, PYALC, everything()) |> 
  tibble()

Métricas

Accuracy (Exactitud): Es una métrica que mide la proporción de predicciones correctas realizadas por un modelo en comparación con el total de predicciones realizadas. Es decir, la cantidad de veces que el modelo acertó sobre el total de datos que se le presentaron.

Sensitivity (Sensibilidad): Es la proporción de verdaderos positivos (TP) que son identificados correctamente por el modelo en relación con el total de verdaderos positivos y falsos negativos (FN). La sensibilidad mide la capacidad del modelo para detectar correctamente los casos positivos.

Specificity (Especificidad): Es la proporción de verdaderos negativos (TN) que son identificados correctamente por el modelo en relación con el total de verdaderos negativos y falsos positivos (FP). La especificidad mide la capacidad del modelo para detectar correctamente los casos negativos.

Métricas

Precision (Precisión): Es la proporción de verdaderos positivos (TP) que son identificados correctamente por el modelo en relación con el total de verdaderos positivos y falsos positivos (FP). La precisión mide la capacidad del modelo para no identificar falsamente un caso como positivo.

Recall (Recuperación): Es la proporción de verdaderos positivos (TP) que son identificados correctamente por el modelo en relación con el total de verdaderos positivos y falsos negativos (FN). El recall mide la capacidad del modelo para identificar todos los casos positivos.

F-measure (Puntuación F): Es una métrica que combina la precisión y el recall en una sola puntuación. El valor de la F-measure oscila entre 0 y 1, siendo 1 el valor óptimo.

Métricas

Kappa (Coeficiente Kappa): Es una medida de concordancia que compara la cantidad de acuerdos observados entre el modelo y las observaciones reales con la cantidad de acuerdos que se esperarían por casualidad. Un valor de kappa cercano a 1 indica una concordancia casi perfecta entre el modelo y las observaciones reales.

Matthews Correlation Coefficient (Coeficiente de Correlación de Matthews): Es una medida que se utiliza para evaluar la calidad de la clasificación binaria. El coeficiente de correlación de Matthews oscila entre -1 y 1, siendo 1 el valor óptimo. Un valor cercano a 1 indica una clasificación perfecta, mientras que un valor cercano a -1 indica una clasificación completamente incorrecta.

Calcular las métricas

conf_mat(resultados_prueba, truth = PYALC,
         estimate = .pred_class)

accuracy(resultados_prueba, truth = PYALC,
         estimate = .pred_class)

sens(resultados_prueba, truth = PYALC,
    estimate = .pred_class)

Calcular las métricas

spec(resultados_prueba, truth = PYALC,
    estimate = .pred_class)

precision(resultados_prueba, truth = PYALC,
    estimate = .pred_class)

recall(resultados_prueba, truth = PYALC,
    estimate = .pred_class)

kap(resultados_prueba, truth = PYALC,
    estimate = .pred_class)

Calcular las métricas

custom_metrics <- metric_set(accuracy, sens, spec, precision, recall, f_meas, kap, mcc)

custom_metrics(resultados_prueba,
  truth = PYALC,
  estimate = .pred_class
)

lm_metrics <- custom_metrics(resultados_prueba,
  truth = PYALC,
  estimate = .pred_class
) |>
  mutate(model = "lm")

Tabla de Métricas

library(gt)

custom_metrics(resultados_prueba,
  truth = PYALC,
  estimate = .pred_class
) |> gt()

Área bajo la curva

prediccion_auc <- predict(lm_results, entrenamiento, type = "prob")

resultados_prueba_auc <- cbind(prediccion_auc, datos_entrenamiento) |>
  tibble()

roc_auc(resultados_prueba_auc,
  truth = PYALC,
  `.pred_No ha consumido`
)

roc_curve(resultados_prueba_auc,
  truth = PYALC,
  `.pred_No ha consumido`
) |> autoplot()

¿Podemos encontrar un modelo que sea mejor que el anterior?