Introducción

Se lee un archivo Excel con 129,373 registros. Se trata de una serie de tiempo de tres años, con el número de accidentes por día en vías.

Las variables tipo fecha usualmente permiten crear nuevas features de tipo categórico, necesarias para los modelos de predicción. El paquete timetk tiene una receta estilo recipes del paradigma tidymodels que realiza el proceso de crear las nuevas features. Sin embargo, hay features que no puede conocer dicha receta, por ejemplo, qué días son puente, semana santa o fin de año. Se añaden.

Machine Learning

Se observará que modeltime1 se acopla al paradigma de tidymodels.

Se usa la función time_series_split() del paquete modeltime para establecer el conjunto de entrenamiento y el de pruebas.

fechas <- fechas %>% arrange(fecha_hecho)
splits <- fechas %>% time_series_split(assess = '3 month', cumulative = TRUE)

En el parámetro assess se establece que los últimos tres meses son el conjunto de prueba.

cumulative = TRUE establece que toda la data anterior es el conjunto de entrenamiento.

Visualización:

splits  %>% 
  tk_time_series_cv_plan() %>%
  plot_time_series_cv_plan(fecha_hecho, cantidad , .interactive = FALSE)

La función tk_time_series_cv_plan() transforma el objeto en un data frame, preparándolo para su visualización.

La función plot_time_series_cv_plan() grafica la serie de tiempo usando las columnas fecha y cantidad.

Modelos tradicionales de series de tiempo

Para demostrar el adecuamiento de los modelos de machine learning a las series de tiempo, se calcula un pronóstico a partir de un modelo ARIMA tradicional.

fit_arima <- modeltime::arima_reg() %>%
  set_engine('auto_arima') %>%
  fit(cantidad ~ fecha_hecho, training(splits))

El comando set_engine especifica que se use la implementación de auto_arima de la librería forecast, la única que acepta arima_reg() a la fecha.

rmse <- sqrt(sum(fit_arima$fit$data$.residuals^2))

La raíz cuadrada del error cuadrático medio es: 637.980128.

festivos <- unique(c(fin_anio, inicio_anio, puentes, sem_santa))
fit_prophet <- modeltime::prophet_reg() %>%
  set_engine('prophet', holidays = festivos) %>%
  fit(cantidad ~ fecha_hecho, training(splits))

Si bien no es un modelo de series de tiempo tradicional, se presenta el algoritmo Prophet, que tampoco es de Machine Learning. Es un algoritmo de código abierto liberado por Facebook2. La idea general del modelo es similar a un generalized additive model (GAM) cuya novedad es que tiene en cuenta los días etiquetados como vacaciones como uno de los componentes, adicionalmente a la estacionalidad y la tendencia. Como es aplicado a series de tiempo, el GAM lo que realiza es una descomposición del modelo aditivo.

Ingeniería de características o preprocesamiento

La siguiente función realiza el pre procesamiento de las features para darle asideros adicionales a los modelos de Machine Learning que se utilizarán.

recipe_spec <- recipe(cantidad ~ fecha_hecho, training(splits)) %>%
  step_timeseries_signature(fecha_hecho) %>%
  step_fourier(fecha_hecho, period = 365, K = 5) %>%
  step_dummy(all_nominal())

La función step_timeseries_signature(), del paquete timetk, divide la variable fecha_hecho en múltiples features que pueden ayudar en los procesos de Machine Learning (por ejemplo, la separa en tres variables adicionales: días, meses y años). Funciona como uno de los ingredientes del paquete recipes.

La función step_fourier(), del mismo paquete, ajusta la serie de tiempo a una serie de fourier, lo cual ayuda al modelamiento si la serie es periódica.

Modelos de Machine Learning para series de tiempo.

Para evitar over-fitting a los datos, se aplican penalizaciones para suavizar el modelo y hacerlo más general.

Elastic_net

Se utilizará una regresión regularizada. La regresión puede ser lineal de orden n, o pueden ser B-splines. Esto significa que se abandona la idea de autocorrelaciones en la serie de tiempo de un dato respecto a datos anteriores. Es decir, se deja de considerar como una serie de tiempo, pero el objetivo es el mismo, hallar un modelo que pronostique la serie.

model_spec_glmnet <- linear_reg(penalty = 0.01, mixture = 0.5) %>% 
  parsnip::set_engine('glmnet')

El comando set_engine especifica que se use la implementación de la librería glmnet.

fit_wflow_glmnet <- workflows::workflow() %>%
  workflows::add_model(model_spec_glmnet) %>%
  workflows::add_recipe(recipe_spec %>% step_rm(fecha_hecho)) %>%
  fit(data = training(splits))

Obsérvese como el comando workflow permite ensamblar la receta del pre procesamiento con la especificación del modelo y entregársela a un comando fit() para que ajuste el modelo sobre unos datos específicos, en este caso, el split reservado para el entrenamiento.

El paquete vip permite especificar la importancia de cada variable:

## # A tibble: 52 x 3
##    Variable                 Importance Sign 
##    <chr>                         <dbl> <chr>
##  1 fecha_hecho_wday.lbl_2        20.3  POS  
##  2 fecha_hecho_month.lbl_06       9.92 POS  
##  3 fecha_hecho_month.lbl_05       4.92 POS  
##  4 fecha_hecho_month.lbl_03       4.41 POS  
##  5 fecha_hecho_month.lbl_09       3.59 POS  
##  6 fecha_hecho_mday7              3.19 POS  
##  7 fecha_hecho_month.lbl_08       2.90 POS  
##  8 fecha_hecho_month.lbl_04       2.82 POS  
##  9 fecha_hecho_wday.lbl_3         2.56 POS  
## 10 fecha_hecho_week.iso           1.19 POS  
## # ... with 42 more rows

Random Forest

Este es un potente algoritmo. También se puede utilizar sobre series de tiempo.

model_spec_rf <- rand_forest(trees = 500, min_n = 50) %>%
  set_engine('randomForest')
fit_wflow_rf <-
  workflows::workflow() %>%
  workflows::add_recipe(recipe_spec %>% step_rm(fecha_hecho)) %>%
  workflows::add_model(model_spec_rf) %>% 
  fit(training(splits))

La idea del Random Forrest es crear un ensamble de modelos en el que se eliminan algunos valores en el tiempo, o algunos términos, o ambos, para la especificación de cada modelo, obtener una variabilidad como resultado y tomar un modelo intermedio a partir del conjunto.

Prophet Boost

El algoritmo Prophet Boost combina Prophet con XGBoost. Primero modela la serie utilizando Prophet, y luego usa los residuales arrojados y las features de la receta en el modelo XGBoost.

model_spec_prophet_boost <- prophet_boost() %>%
  set_engine("prophet_xgboost", yearly.seasonality = TRUE) 

fit_wflow_prophet_boost <- workflows::workflow() %>%
  workflows::add_recipe(recipe_spec) %>% 
  workflows::add_model(model_spec_prophet_boost) %>% 
  fit(training(splits))

XGBoost es un método de baggin.

Administrando varios modelos

La función modeltime_table() del paquete modeltime organiza los modelos con IDs y permite mantenerlos organizados:

model_table <- modeltime_table(
  fit_arima, 
  fit_prophet,
  fit_wflow_glmnet,
  fit_wflow_rf,
  fit_wflow_prophet_boost) 

model_table
## # Modeltime Table
## # A tibble: 5 x 3
##   .model_id .model     .model_desc                      
##       <int> <list>     <chr>                            
## 1         1 <fit[+]>   ARIMA(0,1,1)(1,0,0)[7] WITH DRIFT
## 2         2 <fit[+]>   PROPHET                          
## 3         3 <workflow> GLMNET                           
## 4         4 <workflow> RANDOMFOREST                     
## 5         5 <workflow> PROPHET W/ XGBOOST ERRORS

La función modeltime_calibrate(), del mismo paquete, cuantifica el error y estima el intervalo de confianza.

calibration_table <- model_table %>%
  modeltime_calibrate(testing(splits))

calibration_table
## # Modeltime Table
## # A tibble: 5 x 5
##   .model_id .model     .model_desc                       .type .calibration_data
##       <int> <list>     <chr>                             <chr> <list>           
## 1         1 <fit[+]>   ARIMA(0,1,1)(1,0,0)[7] WITH DRIFT <NA>  <NULL>           
## 2         2 <fit[+]>   PROPHET                           <NA>  <NULL>           
## 3         3 <workflow> GLMNET                            Test  <tibble [90 x 4]>
## 4         4 <workflow> RANDOMFOREST                      Test  <tibble [90 x 4]>
## 5         5 <workflow> PROPHET W/ XGBOOST ERRORS         Test  <tibble [90 x 4]>

Obsérvese que se aplicó sobre testing(splits), para evaluar el modelo. Sólo puede verificar los tres modelos de Machine Learning.

Pronosticando

La función modeltime_forecast() realiza el trabajo de pronóstico.

# calibration_table %>%
calibration_table %>% 
  modeltime_forecast(actual_data = data.frame(fechas)) %>% 
  filter(.model_id %in% c(3:5)) %>% 
  plot_modeltime_forecast(.interactive = FALSE, .conf_interval_show = TRUE) +
  theme(legend.text=element_text(size=7)) + theme(legend.position = "right")

Evaluación

tabla_exactitud <- calibration_table %>%
  modeltime_accuracy() %>%
  table_modeltime_accuracy(.interactive = FALSE)
data.frame(tabla_exactitud[[1]]) %>% knitr::kable(row.names = FALSE)
.model_id .model_desc .type mae mape mase smape rmse rsq
3 GLMNET Test 9.88 31.14 0.77 16.87 13.35 0.78
4 RANDOMFOREST Test 11.17 35.14 0.87 16.51 14.33 0.78
5 PROPHET W/ XGBOOST ERRORS Test 7.56 11.01 0.59 11.35 10.11 0.85

Todas las métricas indican que el algoritmo Prophet Boost es el mejor de los tres modelos.

  1. https://cran.r-project.org/web/packages/modeltime/index.html↩︎

  2. https://facebook.github.io/prophet/docs/quick_start.html#r-api↩︎