Rafael Zambrano: Evaluación de predicciones mediante programación funcional en R

Rafael Zambrano

Evaluación de predicciones 📈

🗒️ En este posts se evaluaran las predicciones de 3 modelos mediante las métricas que estan disponibles en el paquete yardstick.¹ Sin embargo, se mostrará una situación particular para la cual actualmente no existe una función dentro del ecosistema de tidymodes² que permita realizar dicha evaluación 😫. Al menos, no pude encontrar una implementación en la documentación o en los distintos ejemplos publicados, de igual manera, esta será la excusa para crear funciones en R a través de programación funcional 🤩, los temas que se tratarán en esta publicación son los siguientes:

Ajuste y predicción rapida de modelos con tidymodels (no es el foco de este Post).
Evaluación de predicciones con yardstick .
Creación de funciones en R.
Uso de lapply.
Uso de mapply.
Uso de multilapply.
Uso de multieval, función creada y almacenada dentro del paquete sknifedatar.³

Libraries 📚

library(modeldata)
library(tidymodels)
library(tidyverse)

Se hará uso del conjunto de datos car_prices almacenado en el paquete modeldata⁴ , el objetivo es predecir el precio de los vehículos, para conocer más detalles del set de datos consultar la documentación del paquete.

data(car_prices)
rmarkdown::paged_table(car_prices)

Modelado 🚀

La partición de datos se hace con un split 80/20, se normalizan las variables a través de una receta de preprocesamiento, se ajustan 3 modelos predictivos y se realizan predicciones sobre los datos de prueba. 🔎El foco de este Post no es el ajuste o interpretación de modelos, por lo tanto, no es prioritario la ingeniería de características, la performance o la interpretación de los modelos. Para conocer en detalle el modelado dentro de tidymodels consular tidymodels.org.

##===============================================================
##                    Partición de datos                     =
##===============================================================
set.seed(123)
splits = initial_split(car_prices, prop = 0.8, strata = Price)

Receta = recipe(Price ~ ., data = training(splits)) %>%
  step_normalize(all_numeric(),- all_outcomes())

##===============================================================
##                    Definicion de modelos                     =
##===============================================================
mode = "regression"

model_knn = nearest_neighbor() %>%
  set_engine("kknn") %>%
  set_mode(mode)

model_dtr = decision_tree() %>%
  set_engine("rpart") %>%
  set_mode(mode)

model_xgb = boost_tree() %>%
  set_engine("xgboost") %>%
  set_mode(mode)

##===============================================================
##                          Workflows-                           =
##===============================================================
wf_knn = workflow() %>%
  add_recipe(Receta) %>%
  add_model(model_knn)

wf_dtr = workflow() %>%
  add_recipe(Receta) %>%
  add_model(model_dtr)

wf_xgb = workflow() %>%
  add_recipe(Receta) %>%
  add_model(model_xgb)

##================================================================
##                          Predicciones                         =
##================================================================
bind_rows(wf_knn %>% last_fit(splits) %>% collect_metrics() %>% mutate(modelo = "KNN"),
          wf_dtr %>% last_fit(splits) %>% collect_metrics() %>% mutate(modelo = "DCT"),
          wf_xgb %>% last_fit(splits) %>% collect_metrics() %>% mutate(modelo = "XGB"))

# A tibble: 6 x 5
  .metric .estimator .estimate .config              modelo
  <chr>   <chr>          <dbl> <chr>                <chr> 
1 rmse    standard    2251.    Preprocessor1_Model1 KNN   
2 rsq     standard       0.960 Preprocessor1_Model1 KNN   
3 rmse    standard    3923.    Preprocessor1_Model1 DCT   
4 rsq     standard       0.873 Preprocessor1_Model1 DCT   
5 rmse    standard    2137.    Preprocessor1_Model1 XGB   
6 rsq     standard       0.968 Preprocessor1_Model1 XGB

☑️ Para los tres modelos se evalúan las métricas rmse, rsq y mae sobre los datos de validación. Esto es realmente fácil de hacer a través de la función last_fit, permitiendo ajustar el modelo sobre el train y hacer las predicciones sobre el test, train y test están almacenados dentro del objeto splits.

Predicciones 🔮

Pero…¿Qué pasa cuando no es posible usar la funciones last_fit o collect_metrics?

🔎La función last_fit actúa sobre el objeto splits, luego aplicando la función collect_metrics se obtienen las métricas de evaluación. En el caso de querer realizar una comparación de modelos procedentes de distintos frameworks, evaluando predicciones de modelos de tidymodels y de marcos externos, no podemos aplicar estas funciones. Para este escenario puede consultarse una caso de uso en Zambrano (2020).

❗Otra situación que también es esquiva para estas funciones está relacionada con el ajuste sobre los splits, por ejemplo, dentro de un entorno empresarial es común que exista la necesidad de colocar el modelo en producción ⚙️, en este caso, el modelo final no debe ajustarse sobre los splits, en su lugar, debe hacerse sobre todo el dataset y estar listo para hacer predicciones para nuevas observaciones.

En estas situaciones, las funciones last_fit y collect_metrics no pueden ayudarnos (hasta donde pude investigar), debemos realizar el proceso manualmente para cada métrica 😭😭😭. A continuación, se ajustan los 3 modelos sobre todo el conjunto de datos, luego a modo de ejemplo se realizan las predicciones sobre el mismo dataset.

##================================================================
##                  Ajuste sobre todos los datos                 =
##================================================================

knn_predict = wf_knn %>% 
  fit(car_prices) %>% 
  predict(car_prices) %>% 
  rename(KNN = .pred)

dtr_predict = wf_dtr %>% 
  fit(car_prices) %>% 
  predict(car_prices) %>% 
  rename(DCT = .pred)

xgb_predict = wf_xgb %>% 
  fit(car_prices) %>% 
  predict(car_prices) %>% 
  rename(XGB = .pred)

predicciones = bind_cols(car_prices %>% select(Price), knn_predict, dtr_predict, xgb_predict)

predicciones

# A tibble: 804 x 4
    Price    KNN    DCT    XGB
    <dbl>  <dbl>  <dbl>  <dbl>
 1 22661. 21084. 19127. 20383.
 2 21725. 22265. 19127. 19888.
 3 29143. 30234. 29495. 29046.
 4 30732. 30845. 29495. 31397.
 5 33359. 33751. 29495. 33970.
 6 30315. 30874. 29495. 31166.
 7 33382. 33757. 29495. 33156.
 8 30251. 31994. 29495. 30269.
 9 30167. 31433. 29495. 31026.
10 27060. 27211. 29495. 27440.
# … with 794 more rows

Evaluación de pedicciones 🎯

Los datos observados (columna Price) y las predicciones de los 3 modelos se encuentran en un data frame, a través de la función rmse del paquete yardstick se calcula el rmse del modeloKNN.

rmse(data = predicciones,
     truth    = Price,
     estimate = KNN,
     na_rm    = TRUE)

# A tibble: 1 x 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       1293.

❗La función recibe como primer argumento el data frame y luego las columnas a evaluar, se replica el mismo proceso para la métrica mae sobre el modelo KNN y el rmse para el XGB.

mae(data = predicciones,
     truth    = Price,
     estimate = KNN,
     na_rm    = TRUE)

# A tibble: 1 x 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 mae     standard        853.

rmse(data = predicciones,
     truth    = Price,
     estimate = XGB,
     na_rm    = TRUE)

# A tibble: 1 x 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 rmse    standard       1554.

Perfecto 👌🏻👌🏻👌🏻. Ahora debemos calcular las métricas mae, rsq y rmse sobre los 3 modelos. ¿Cómo se podrían estimar? 🤔

Evaluación manual 😓😓😓

Son 3 métricas para 3 modelos, por lo tanto se tendrían que implementar 9 veces las funciones de evaluación 😨😨😨 . A continuación, se muestra el ejemplo de cómo se podría hacer esta evaluación de forma manual.

Metricas =
bind_rows(
##---------------KNN------------------------------------------------
  rmse(data = predicciones,
       truth    = Price,
       estimate = KNN,
       na_rm    = TRUE) %>%
    mutate(  modelo = "KNN"),
  
  mae(data = predicciones,
      truth    = Price,
      estimate = KNN,
      na_rm    = TRUE) %>%
    mutate(modelo = "KNN"),
  
    rsq(data = predicciones,
      truth    = Price,
      estimate = KNN,
      na_rm    = TRUE) %>%
    mutate(modelo = "KNN"),
##------------------DCT---------------------------------------------
    rmse(data = predicciones,
       truth    = Price,
       estimate = DCT,
       na_rm    = TRUE) %>%
    mutate(  modelo = "DCT"),
  
  mae(data = predicciones,
      truth    = Price,
      estimate = DCT,
      na_rm    = TRUE) %>%
    mutate(modelo = "DCT"),
  
    rsq(data = predicciones,
      truth    = Price,
      estimate = DCT,
      na_rm    = TRUE) %>%
    mutate(modelo = "DCT"),
##----------------XGB-----------------------------------------------
    rmse(data = predicciones,
       truth    = Price,
       estimate = XGB,
       na_rm    = TRUE) %>%
    mutate(  modelo = "XGB"),
  
  mae(data = predicciones,
      truth    = Price,
      estimate = XGB,
      na_rm    = TRUE) %>%
    mutate(modelo = "XGB"),
  
    rsq(data = predicciones,
      truth    = Price,
      estimate = XGB,
      na_rm    = TRUE) %>%
    mutate(modelo = "XGB")
)

Metricas

# A tibble: 9 x 4
  .metric .estimator .estimate modelo
  <chr>   <chr>          <dbl> <chr> 
1 rmse    standard    1293.    KNN   
2 mae     standard     853.    KNN   
3 rsq     standard       0.983 KNN   
4 rmse    standard    2969.    DCT   
5 mae     standard    2332.    DCT   
6 rsq     standard       0.910 DCT   
7 rmse    standard    1554.    XGB   
8 mae     standard    1131.    XGB   
9 rsq     standard       0.976 XGB

Bien 🤨, en un data frame se almacena la evaluación de las 3 métricas para los 3 modelos, sin embargo, es un proceso manual, extenso, propenso a errores y no es reproducible para otro conjunto de predicciones

Evaluación funcional 🧐🧐🧐

🔎El último proceso de evaluación es un código repetitivo que se podría generar a través de una función, sin embargo, tiene un cierto grado de complejidad. El objetivo final será tener una función, que en una línea de código y a través de 3 parámetros, permite calcular distintas métricas disponibles en el paquete yardstick .

Además, con las ventajas de ser independiente del modelo, nombre de columna o paquete utilizado para generar las predicciones, únicamente se necesitan las predicciones almacenadas en un data frame 👌🏻.

🔷 Función 1: una métrica sobre una predicción

Partimos de una función sencilla, esta nos permitirá calcular una métrica para una predicción, a continuación su implementación.

multieval = function(data , observed, predictions, metrica){
  
  metrica(data = data,
       truth    = .data[[observed]],
       estimate = .data[[predictions]],
       na_rm    = TRUE) %>% 
    mutate(modelo = predictions)
}

multieval(data = predicciones, observed = "Price" ,predictions = "XGB", metrica = rmse)

# A tibble: 1 x 4
  .metric .estimator .estimate modelo
  <chr>   <chr>          <dbl> <chr> 
1 rmse    standard       1554. XGB

❗Básicamente recibe el data frame de predicciones, toma la columna observada, la predicción del modelo y calcula la métrica deseada.

🔷 Función 2: una métrica sobre “n” predicicones

Ahora se extenderá la función anterior para hacer el cálculo de una métrica sobre varias predicciones, la función conserva la mismas estructura.

multieval = function(data , observed, predictions, metrica){
  
  names(predictions) = predictions
  
  lapply(predictions, function(x){
    
    metrica(data = data,
            truth    = .data[[observed]],
            estimate = .data[[x]],
            na_rm    = TRUE) %>% 
      mutate(modelo = x)
    })
}

multieval(data = predicciones, observed = "Price" , predictions = c("XGB","KNN","DCT"), metrica = rmse) %>% 
  bind_rows()

# A tibble: 3 x 4
  .metric .estimator .estimate modelo
  <chr>   <chr>          <dbl> <chr> 
1 rmse    standard       1554. XGB   
2 rmse    standard       1293. KNN   
3 rmse    standard       2969. DCT

🔎En la mitad del código se observa que fue incorporada una función lapply, esta nos permite aplicar una función (la métrica de evaluación) sobre una lista o vector (en este caso el nombre de las columnas donde se encuentran las predicciones). Se observa el valor del rmse para los 3 modelos.

🔷 Función 3: “n” métricas sobre una prediccion

Siguiendo la estructura de la función anterior, se aplica la función lappy pero siendo la lista el conjunto de funciones de evaluación.

library(erer)
multieval = function(data , observed, predictions, metrica){
  
  lapply(metrica, function(x){
    
    x(data = data,
            truth    = .data[[observed]],
            estimate = .data[[predictions]],
            na_rm    = TRUE) %>% 
      mutate(modelo = predictions)
  }
  )
}

multieval(predicciones, "Price" , "XGB", metrica = listn(rmse,rsq,mae)) %>% 
  bind_rows()

# A tibble: 3 x 4
  .metric .estimator .estimate modelo
  <chr>   <chr>          <dbl> <chr> 
1 rmse    standard    1554.    XGB   
2 rsq     standard       0.976 XGB   
3 mae     standard    1131.    XGB

❗Para el modelo XGB se calculan las 3 métricas.

🧪 Combinación de funciones

🔎El objetivo es iterar sobre todas las combinaciones de predicciones y métricas en una misma función, la función lapply por sí sola únicamente permite iterar sobre un objeto 🤨. Ahora se probará con la función mapply, esta permite iterar sobre múltiples objetos, a continuación, se muestra la implementación iterando sobre métricas y predicciones simultáneamente.

multieval = function(data , observed, predictions, metrica){
  
  names(predictions) = predictions
  
  mapply(function(x, y ){
    
    x(data = data,
      truth    = .data[[observed]],
      estimate = .data[[y]],
      na_rm    = TRUE) %>%
      mutate(modelo = y)
    },
    metrica,predictions, SIMPLIFY = F)
  }

multieval(predicciones, "Price", c("XGB","KNN","DCT"), metrica = listn(rmse,rsq,mae)) %>% 
  bind_rows()

# A tibble: 3 x 4
  .metric .estimator .estimate modelo
  <chr>   <chr>          <dbl> <chr> 
1 rmse    standard    1554.    XGB   
2 rsq     standard       0.983 KNN   
3 mae     standard    2332.    DCT

😭😭😭 …se observan 3 métricas y 3 modelos, la salida debía tener 9 elementos, mapply aunque admite mas de un argumento para iterar, no permite realizar todas las combinaciones de ambos objetos, mapply no parece ser una solución suficiente para nuestro problema.

🔶 Función final

💡Indagando en distintos paquetes y en las documentaciones de las familias de funciones *apply(), no encontré una función similar a mapply que permitiera iterar sobre toda la combinación de elementos. Sin embargo, aunque puede construirse dicha función (🤔🧠🤓 se podría hacer un arreglo anidando dos funciones lapply ), en un repositorio público de Gist GitHub propietario de Alek Rutkowski’s, se encontró una función denominada mlapply, la misma actúa como un lapply múltiple y se adapta a la solución que se necesita 🤩.

A continuación, se muestra la estructura de la función, para visitar el repositorio original donde está almacenada puede consultarse el siguiente enlace, repositorio Gist GitHub de la función mlapply.

mlapply <- function(.Fun, ..., .Cluster=NULL, .parFun=parallel::parLapply) {
    `--List--` <-
        list(...)
    names(`--List--`) <-
        names(`--List--`) %>% 
        `if`(is.null(.),
             rep.int("", length(`--List--`)),
             .) %>% 
        ifelse(.=="", # for unnamed args in ...
               seq_along(.) %>% 
                   paste0(ifelse(.==1 | .>20 & .%%10==1, 'st', ""),
                          ifelse(.==2 | .>20 & .%%10==2, 'nd', ""),
                          ifelse(.==3 | .>20 & .%%10==3, 'rd', ""),
                          ifelse(.>3 & .<=20 | !(.%%10 %in% 1:3), 'th', "")) %>% 
                   paste("argument in mlapply's ..."),
               .)
    `--metadata--` <-
        data.frame(Name = paste0("`",names(`--List--`),"`"),
                   Len = lengths(`--List--`),
                   OriginalOrder = seq_len(length(`--List--`)),
                   stringsAsFactors=FALSE)
    eval(Reduce(function(previous,x)
        paste0('unlist(lapply(`--List--`$',x,',',
               'function(',x,')', previous,'),recursive=FALSE)'),
        x =
            `--metadata--` %>% 
            `[`(order(.$Len),) %>% 
            `$`(Name),
        init =
            `--metadata--` %>% 
            `[`(order(.$OriginalOrder),) %>% 
            `$`(Name) %>% 
            ifelse(grepl("argument in mlapply's ...",.,fixed=TRUE),
                   ., paste0(.,'=',.)) %>% 
            paste(collapse=',') %>%
            paste0('list(.Fun(',.,'))')) %>% 
            ifelse(.Cluster %>% is.null,
                   .,
                   sub('lapply(',
                       '.parFun(.Cluster,',
                       ., fixed=TRUE)) %>% 
            parse(text=.))
}

❗Un poco extensa y con tópicos de programación funcional bastante avanzados Alek Rutkowski’s construyó una función que se adapta a nuestro problema. ⏯️ A continuación, se muestra la implementación final de la función multieval, solo se necesita suministrar el data frame de predicciones y especificar dentro del data frame cuál es la columna observada y cuales son las predicciones 👌🏻. Finalmente se especifica una lista de métricas a evaluar.

library(sknifedatar)
library(erer)

multieval = function(data , observed, predictions, metrica){

  names(predictions) = predictions

  mlapply(function(x, y ){

    x(data = data,
      truth    = .data[[observed]],
      estimate = .data[[y]],
      na_rm    = TRUE) %>%
      dplyr::mutate(model = y)
  },
  metrica,predictions
  ) %>%
    dplyr::bind_rows() %>%
    dplyr::arrange(.data$.metric)
}

evaluation <- multieval(predicciones, "Price", c("XGB","KNN","DCT"), metrica = listn(rmse,rsq,mae))

evaluation

# A tibble: 9 x 4
  .metric .estimator .estimate model
  <chr>   <chr>          <dbl> <chr>
1 mae     standard    1131.    XGB  
2 mae     standard     853.    KNN  
3 mae     standard    2332.    DCT  
4 rmse    standard    1554.    XGB  
5 rmse    standard    1293.    KNN  
6 rmse    standard    2969.    DCT  
7 rsq     standard       0.976 XGB  
8 rsq     standard       0.983 KNN  
9 rsq     standard       0.910 DCT

Perfecto 😎😎 la función hace justo lo que buscamos 💪🏻. Se recomienda suministrar la lista de métricas a través de la función listn del paquete erer en lugar de la función list del base de R. Esto simplemente por un tema visual, ya que al usar list tendríamos que nombrar cada métrica que se desee evaluar ,por ejemplo, “list(rmse = rmse, mae = mae, ….),” de lo contrario el elemento de la lista no tendría un nombre asignado.

❗Finalmente se aplicaron algunos conceptos de programación funcional en R, utilizando las funciones lapply y mapply, este es un tema central dentro del lenguaje R, para un mayor entendimiento se recomienda consultar a Wickham (2019) . También se presentó el funcionamiento de la función multieval, que puede ser de utilidad para evaluar predicciones en la situaciones planteadas al inicio del artículo.

📌Como nota final, a través del paquete purrr⁵ se pueden encontrar implementaciones más simples y robustas de la familia de funciones *apply().

🌐 Siéntase libre de comentar , sugerir o compartir cualquier idea.

Muchas gracias 👏🏻👏🏻👏🏻

Contacto ✉

Rafael Zambrano, Linkedin, Twitter, Github, Blogpost.

Acknowledgments

Estudio desarrollado mediante el lenguaje estadístico R (R Core Team 2020) .

Henry, Lionel, and Hadley Wickham. 2020. Purrr: Functional Programming Tools. https://CRAN.R-project.org/package=purrr.

Kuhn, Max. 2020. Modeldata: Data Sets Used Useful for Modeling Packages. https://CRAN.R-project.org/package=modeldata.

Kuhn, Max, and Davis Vaughan. 2020. Yardstick: Tidy Characterizations of Model Performance. https://CRAN.R-project.org/package=yardstick.

Kuhn, Max, and Hadley Wickham. 2020. Tidymodels: A Collection of Packages for Modeling and Machine Learning Using Tidyverse Principles. https://www.tidymodels.org.

R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Wickham, Hadley. 2019. Advanced r. Chapman; Hall/CRC. https://doi.org/10.1201/9781351201315.

Zambrano, Rafael. 2020. “Rafael Zambrano: Predicción de Delitos En CABA.” https://rafael-zambrano-blog-ds.netlify.app/posts/2020-12-22-prediccin-de-delitos-en-caba/.

Zambrano, Rafael, and Karina Bartolome. 2021. Sknifedatar: Swiss Knife of Data for r. https://github.com/rafzamb/sknifedatar.

Kuhn and Vaughan (2020)↩︎
Kuhn and Wickham (2020)↩︎
Zambrano and Bartolome (2021)↩︎
Kuhn (2020)↩︎
Henry and Wickham (2020)↩︎

Comment on this article Share:

Evaluación de predicciones mediante programación funcional en R