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(): Estimar un modelo de regresión lineal.

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

glmnet(): Estimar 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(), decision_tree()
  • Especificar el tipo de outcome
    • Regresión: variables continuas
    • Clasificación: categórica,multinomial, ordinal

Especificar el modelo en tidymodels

library(tidymodels)

tidymodels_prefer()

linear_reg() |> 
  set_engine("lm") |> 
  set_mode("regression")

logistic_reg() |> 
    set_engine("glmnet") |> 
    set_mode("classification")

linear_reg() |>
    set_engine("stan")  |> 
    set_mode("regression")

Especificar el modelo en tidymodels

Tip

library(): cargar el paquete tidymodels.

tidymodels_prefer(): establecer el paquete tidymodels como el preferido.

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

logistic_reg() : Crear una especificación de modelo de regresión logística en el marco de tidymodels.

set_engine(): establecer el motor (paquete) de cálculo del modelo. Se utilizará “lm” para el primer modelo, “glmnet” para el segundo y “stan” para el tercero.

Estimar un modelo

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

Descriptivos


Calcular la tabla de descriptivos

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

Calcular la tabla de descriptivos por grupo (Sí y No ha consumido)

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

Explorar los datos visualmente

Hacer un gráfico para ver tendencias en los factores de riesgo

mini_datos |> 
 select(PYALC,factores_de_riesgo) |>
  pivot_longer(-PYALC) |> 
  ggplot(aes(value, fct_reorder(name, value), fill = PYALC)) +
  geom_boxplot()

Especificar el modelo

1. Especificar el modelo

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 especifica que es un problema de “classification”, debido a que la variable dependiente es categórica. Todo esto se almacena en el objeto lm_model.


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

2. Estimar el modelo

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


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

3. Obtener los resultados del modelo

Obtener los coeficientes (betas) del modelo

lm_results |> tidy() 

Obtener los índices de ajuste del modelo

lm_results |> glance()

3. Obtener los resultados del modelo

Tip

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.

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

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

Práctica

  1. Estime la asociación entre el consumo de marihuana en los últimos 30 días (MAR30DY), a partir de los factores de riesgo.

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 conjuntos. Uno para entrenar el modelo y otro para evaluarlo.

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.

Note

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

Esto no puede pasar: contaminación de los datos de prueba (information leakage)

Dividir los datos

Crear dos bases: entrenamiento 80% y prueba 20%

set.seed(1234)

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(): 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.

initial_split(): dividir los datos en conjuntos de entrenamiento y prueba. Divide el objeto “mini_datos” en conjuntos de entrenamiento y prueba en una proporción de 80/20, y estratificando por la columna “PYALC”.

training(): extraer el conjunto de entrenamiento de un objeto creado con la función initial_split(). Se utiliza para extraer el conjunto de entrenamiento del objeto “datos_divididos”.

testing(): extraer el conjunto de prueba de un objeto creado con la función initial_split(). Se utiliza para extraer el conjunto de prueba del objeto “datos_divididos”.

Dividir los datos

Descripivos de los datos de entrenamiento

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

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

1. Especificar el modelo

No es necesario volver a especificar el modelo, ya que ya lo hicimos en la parte I. Pero si lo hiciéramos, sería de la siguiente manera:

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

2. Estimar el modelo

Estimar el modelo con la función fit(). Se especifica la variable de resultado (PYALC) y los predictores (.). Se utiliza el objeto datos_entrenamiento como base de datos.

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

3. Obtener los resultados del modelo

Obtener los coeficientes (betas) del modelo en OR con sus intervalos de confianza.

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

4. Evaluar el modelo

4. Evaluar el modelo

Hay varias formas de evaluar un modelo. Algunas de las más comunes son: - Métricas de ajuste - Área bajo la curva (AUC) - Sensibilidad y especificidad - Matrices de confusión …

Dado que nuestro objevito es predecir el consumo de alcohol, vamos a utilizar las métricas que evaluen el modelo según su desempeño en la clasificación de observaciones (accuracy o exactitud).

4. Evaluar el modelo

Utilizar el modelo para predecir los valores de la variable dependiente en los datos de entrenamiento.

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

predict(lm_results, predecir_estos_valores)

Obtener las predicciones en forma de probabilidades

predict(lm_results, predecir_estos_valores, type="prob")

4. Evaluar el modelo

Obtener las predicciones y los valores observadods para calcular algunas métricas.

entrenamiento <- datos_entrenamiento %>% 
    select(-PYALC)

prediccion <- predict(lm_results, entrenamiento)

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

resultados_prueba

Métricas

Métricas

Exactitud (accuracy): 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.

Sensibilidad (sensitivity/recall): 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.

Especificidad (Specificity): 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.

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.

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.

Métricas

Calcular las métricas

Con los datos observados y los predichos, se pueden calcular las métricas del modelo.

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, event_level = "second")

Calcular las métricas

Más métricas

spec(resultados_prueba, truth = PYALC,
    estimate = .pred_class, event_level = "second")

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

recall(resultados_prueba, truth = PYALC,
    estimate = .pred_class, event_level="second")

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

Calcular las métricas

Calcular todas las métricas y guardarlas en una base de datos

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

lm_metrics  <- custom_metrics(resultados_prueba,
  truth = PYALC,
  estimate = .pred_class, 
  event_level = "second"
)  |> 
mutate(model= "Regresión Logística")

Tabla de Métricas

library(gt)

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

Área bajo la curva

El área bajo la curva (AUC) es una métrica que mide la capacidad de un modelo para distinguir entre dos clases. En el caso de un modelo de regresión logística, el AUC mide la capacidad del modelo para distinguir entre los que han consumido alcohol y los que no.

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_Sí ha consumido`,
  event_level = "second"
)

roc_curve(resultados_prueba_auc,
  truth = PYALC,
   `.pred_Sí ha consumido`,
  event_level = "second"
) |> autoplot()

Práctica

¿Cuál es la exactitud y la sensibilidad del modelo que estimaron en la práctica anterior?

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