Introducción

Este documento hace parte del trabajo del curso de Analítica Predictiva de la Universidad Nacional de Colombia para la Maestría en Ingeniería y Especialización en Analítica. El alcance de este es la creación de modelos predictivos para la estimación de cantidad de accidentes en el valle del aburra a nivel semanal.

Carga de librerias

# Manipulación de datos
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
# Data modelling
library(recipes)
## 
## Attaching package: 'recipes'
## The following object is masked from 'package:stats':
## 
##     step
library(Metrics)
library(glmnet)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## Loaded glmnet 4.0-2
library(caret)
## Loading required package: lattice
## Loading required package: ggplot2
## 
## Attaching package: 'caret'
## The following objects are masked from 'package:Metrics':
## 
##     precision, recall
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(rpart)
library(rpart.plot)
library(e1071)

# Visualización
library(ggplot2)

# Formato tablas
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows

Carga de datos

df_train <- read.csv('../../data/processed/train_data.csv')
df_test <- read.csv('../../data/processed/test_data.csv')
special_dates_train <- read.csv('../../data/processed/special_date_monthly_2017.csv')
special_dates_test <- read.csv('../../data/processed/special_date_monthly_2018.csv')
kable(head(df_train)) %>%
  kable_styling(bootstrap_options = "striped", full_width = F) %>%
  scroll_box(width = "100%")
X DIA PERIODO CLASE DIRECCION DIRECCION_ENC CBML TIPO_GEOCOD GRAVEDAD BARRIO COMUNA DISENO DIA_NOMBRE MES LONGITUD LATITUD FECHA
1 1 2014 choque cr 63 cl 94 cr 063 094 000 00000 no ubicada solo danos tramo de via miercoles 1 -75.70382 6.221806 2014-01-01 19:00:00
2 1 2014 choque cl 30 cr 66 b cl 030 066 b 000 00000 1602 malla vial solo danos rosales belen interseccion miercoles 1 -75.58727 6.231716 2014-01-01 07:40:00
3 1 2014 choque cr 52 cl 97 cr 052 097 000 00000 0402 malla vial solo danos san isidro aranjuez interseccion miercoles 1 -75.56253 6.289907 2014-01-01 05:30:00
4 1 2014 choque tv 78 cl 65 tv 078 065 000 00000 0519 malla vial solo danos el progreso castilla tramo de via miercoles 1 -75.57365 6.275473 2014-01-01 13:50:00
5 1 2014 otro cr 63 cl 50 cr 063 050 000 00000 1101 malla vial solo danos carlos e. restrepo laureles estadio tramo de via miercoles 1 -75.57697 6.255457 2014-01-01 07:25:00
6 1 2014 choque cr 57 cl 51 cr 057 051 000 00000 1006 malla vial solo danos san benito la candelaria tramo de via miercoles 1 -75.57481 6.254322 2014-01-01 04:15:00

Armado del conjunto de datos

Debido a que se deben armar ciertas agregaciones a nivel de conjunto de datos antes de crear los datasets de entrenamiento, se incluye una etapa de armado del conjunto de variables, no está asociado directamente con la limpieza de los datos. La base que se construye a continuación serán los cimientos para el modelo semanal.

df <- rbind(
  df_train,
  df_test
)
special_dates_train <- special_dates_train %>% 
  dplyr::select(PERIODO, MES, DIA, DIA_FESTIVO, FECHA_ESPECIAL, FESTIVO_FECHA_ESPECIAL)
special_dates_test <- special_dates_test %>% 
  dplyr::select(PERIODO, MES, DIA, DIA_FESTIVO, FECHA_ESPECIAL, FESTIVO_FECHA_ESPECIAL)

special_dates <- rbind(
  special_dates_train,
  special_dates_test
)

df <- merge(df, special_dates, by = c("PERIODO", "MES", "DIA"))

kable(head(df)) %>%
  kable_styling(bootstrap_options = "striped", full_width = F) %>%
  scroll_box(width = "100%")
PERIODO MES DIA X CLASE DIRECCION DIRECCION_ENC CBML TIPO_GEOCOD GRAVEDAD BARRIO COMUNA DISENO DIA_NOMBRE LONGITUD LATITUD FECHA DIA_FESTIVO FECHA_ESPECIAL FESTIVO_FECHA_ESPECIAL
2014 1 1 1 choque cr 63 cl 94 cr 063 094 000 00000 no ubicada solo danos tramo de via miercoles -75.70382 6.221806 2014-01-01 19:00:00 1 1 1
2014 1 1 2 choque cl 30 cr 66 b cl 030 066 b 000 00000 1602 malla vial solo danos rosales belen interseccion miercoles -75.58727 6.231716 2014-01-01 07:40:00 1 1 1
2014 1 1 3 choque cr 52 cl 97 cr 052 097 000 00000 0402 malla vial solo danos san isidro aranjuez interseccion miercoles -75.56253 6.289907 2014-01-01 05:30:00 1 1 1
2014 1 1 4 choque tv 78 cl 65 tv 078 065 000 00000 0519 malla vial solo danos el progreso castilla tramo de via miercoles -75.57365 6.275473 2014-01-01 13:50:00 1 1 1
2014 1 1 5 otro cr 63 cl 50 cr 063 050 000 00000 1101 malla vial solo danos carlos e. restrepo laureles estadio tramo de via miercoles -75.57697 6.255457 2014-01-01 07:25:00 1 1 1
2014 1 1 6 choque cr 57 cl 51 cr 057 051 000 00000 1006 malla vial solo danos san benito la candelaria tramo de via miercoles -75.57481 6.254322 2014-01-01 04:15:00 1 1 1

Modelado semanal

Clase choque

Preprocesamiento de datos

Para el preprocesamiento de datos, se tiene en cuenta variables como rezagos de accidentes de la última semana de todas las clases y gravedades y las últimas 2 y 3 semanas de la clase de interes, al igual que accidentes en días como domingos y número de días festivos y/o especiales en la semana de interes.

df <- df %>%
  mutate(SEMANA=week(FECHA))
special_dates <- special_dates %>% 
  mutate(
    FECHA = as_date(paste(
      PERIODO, 
      formatC(MES, width = 2, flag = 0),
      formatC(DIA, width = 2, flag = 0),
      sep='-'
      ))
  ) %>% 
  mutate(
    SEMANA=week(FECHA)
  )

armado_dataset_semanal <- function(df, objective_var){
  # conjunto de datos semanal por clase
  accidentes_por_clase <- df %>% 
    group_by(PERIODO, SEMANA, CLASE) %>% 
    summarise(n_accidentes=n(), .groups="drop") %>% 
    spread(CLASE, n_accidentes) %>%
    mutate(PERIODO_LEAD=lead(PERIODO), SEMANA_LEAD=lead(SEMANA), rank = 1:length(PERIODO)) %>%
    filter(rank < max(rank)) %>%
    dplyr::select(-c(PERIODO, SEMANA, rank)) 
  
  # conjunto de datos semanal por gravedad
  accidentes_por_gravedad <- df %>% 
    group_by(PERIODO, SEMANA, GRAVEDAD) %>%
    summarise(n_accidentes=n(), .groups="drop") %>% 
    spread(GRAVEDAD, n_accidentes) %>%
    mutate(PERIODO_LEAD=lead(PERIODO), SEMANA_LEAD=lead(SEMANA), rank = 1:length(PERIODO)) %>%
    filter(rank < max(rank)) %>%
    dplyr::select(-c(PERIODO, SEMANA, rank))
  
  # rezagos de la clase
  rezagos_accidentes <- df %>% 
    filter(CLASE == objective_var) %>% 
    group_by(PERIODO, SEMANA) %>%
    summarise(n_accidentes=n(), .groups="drop") %>%
    mutate(
      t_minus_1=lag(n_accidentes, n = 1),
      t_minus_2=lag(n_accidentes, n = 2),
      t_minus_3=lag(n_accidentes, n = 3)
    ) %>% 
    dplyr::select(PERIODO, SEMANA, t_minus_2, t_minus_3) 
  
  # Accidentes en los domingos
  domingos <- df %>% 
    mutate(domingo = ifelse(DIA_NOMBRE == "domingo  ", 1, 0)) %>% 
    group_by(PERIODO, SEMANA) %>% 
    summarise(accidentes_domingo=sum(domingo), .groups="drop") %>% 
    mutate(PERIODO_LEAD=lead(PERIODO), SEMANA_LEAD=lead(SEMANA), rank = 1:length(PERIODO)) %>%
    filter(rank < max(rank)) %>% 
    dplyr::select(-c(PERIODO, SEMANA, rank))
  
  
  # fechas especiales
  fechas_especiales <- special_dates %>% 
    group_by(PERIODO, SEMANA) %>% 
    summarise(
      dias_festivos=sum(DIA_FESTIVO),
      dias_especiales=sum(FECHA_ESPECIAL),
      dias_festivos_especiales=sum(FESTIVO_FECHA_ESPECIAL),
      .groups="drop"
    )
  
  df_processed <- df %>%
    filter(CLASE == objective_var) %>%
    group_by(PERIODO, SEMANA) %>% 
    summarise(
      numero_accidentes = n(),
      .groups="drop"
    ) %>% 
    mutate(
      SEMANA=factor(SEMANA)
      ) %>% 
    merge(accidentes_por_clase,
          by.x = c("PERIODO", "SEMANA"), by.y = c("PERIODO_LEAD", "SEMANA_LEAD"), all.x = T) %>%
    merge(accidentes_por_gravedad, by.x = c("PERIODO", "SEMANA"), by.y = c("PERIODO_LEAD", "SEMANA_LEAD"), all.x = T) %>%
    merge(rezagos_accidentes, by = c("PERIODO", "SEMANA"), all.x = T) %>%
    merge(domingos, by.x = c("PERIODO", "SEMANA"), by.y = c("PERIODO_LEAD", "SEMANA_LEAD"), all.x = T) %>% 
    replace_na(list(incendio=0, muerto=0)) %>%
    filter(!is.na(t_minus_3)) %>%
    merge(fechas_especiales, by = c("PERIODO", "SEMANA")) %>% 
    arrange(PERIODO, SEMANA, .by_group = T)
  
  return(df_processed)
}

Se procede a realizar la construcción del conjunto de datos para la clase choque y hacer todo el pipeline de preprocesamiento para implementar el escalado y estandarización al conjunto de entrenamiento y validación.

df_preprocessed <- armado_dataset_semanal(df=df, objective_var = "choque")
df_choque_train <- df_preprocessed %>% filter(PERIODO!=2018) %>% dplyr::select(-c(PERIODO))
df_choque_test <- df_preprocessed %>% filter(PERIODO==2018) %>% dplyr::select(-c(PERIODO))

model.choque.recipe <- recipe(numero_accidentes ~ ., 
                              data = df_choque_train)
model.choque.steps <- model.choque.recipe %>% 
  step_log(all_outcomes()) %>% 
  step_center(all_predictors(), -SEMANA) %>% 
  step_scale(all_predictors(), -SEMANA)

model.choque.prepared <- prep(model.choque.steps, training = df_choque_train)
df_choque_train <- bake(model.choque.prepared, df_choque_train)
df_choque_test <- bake(model.choque.prepared, df_choque_test)

Modelamiento

Regresión lineal simple

Como primer intento, vale la pena evaluar un modelo de regresión lineal simple con stepwise variable selection para descartar variables que no entreguen valor.

modelo.choque <- lm(numero_accidentes ~ ., data=df_choque_train)
step.model.choque <- stepAIC(modelo.choque, direction = "both", trace = F, steps = 2000)

model.choque.summary <- summary(step.model.choque)
model.choque.summary
## 
## Call:
## lm(formula = numero_accidentes ~ SEMANA + choque + volcamiento + 
##     herido + muerto + dias_festivos + dias_especiales + dias_festivos_especiales, 
##     data = df_choque_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.31931 -0.04272  0.00204  0.04418  0.65037 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               6.091776   0.105166  57.925  < 2e-16 ***
## SEMANA2                   0.120890   0.098511   1.227 0.221695    
## SEMANA3                   0.239065   0.104498   2.288 0.023559 *  
## SEMANA4                   0.211269   0.115874   1.823 0.070267 .  
## SEMANA5                   0.221504   0.114468   1.935 0.054874 .  
## SEMANA6                   0.257169   0.115558   2.225 0.027553 *  
## SEMANA7                   0.421244   0.124817   3.375 0.000942 ***
## SEMANA8                   0.245216   0.118518   2.069 0.040271 *  
## SEMANA9                   0.279817   0.122708   2.280 0.024005 *  
## SEMANA10                  0.390873   0.123614   3.162 0.001899 ** 
## SEMANA11                  0.267459   0.123258   2.170 0.031598 *  
## SEMANA12                  0.427632   0.128646   3.324 0.001117 ** 
## SEMANA13                  0.277290   0.111819   2.480 0.014259 *  
## SEMANA14                  0.280529   0.113682   2.468 0.014731 *  
## SEMANA15                  0.286883   0.111210   2.580 0.010857 *  
## SEMANA16                  0.232139   0.111059   2.090 0.038295 *  
## SEMANA17                  0.317643   0.112979   2.812 0.005595 ** 
## SEMANA18                  0.305886   0.119138   2.567 0.011228 *  
## SEMANA19                  0.422618   0.118738   3.559 0.000499 ***
## SEMANA20                  0.299458   0.121294   2.469 0.014685 *  
## SEMANA21                  0.201198   0.117870   1.707 0.089915 .  
## SEMANA22                  0.266244   0.121587   2.190 0.030098 *  
## SEMANA23                  0.269983   0.116036   2.327 0.021328 *  
## SEMANA24                  0.287876   0.117162   2.457 0.015155 *  
## SEMANA25                  0.305695   0.119857   2.551 0.011766 *  
## SEMANA26                  0.223210   0.117207   1.904 0.058785 .  
## SEMANA27                  0.239772   0.110271   2.174 0.031254 *  
## SEMANA28                  0.248458   0.114295   2.174 0.031297 *  
## SEMANA29                  0.329809   0.115071   2.866 0.004757 ** 
## SEMANA30                  0.278489   0.115904   2.403 0.017502 *  
## SEMANA31                  0.451454   0.119779   3.769 0.000236 ***
## SEMANA32                  0.354359   0.138585   2.557 0.011558 *  
## SEMANA33                  0.329505   0.121372   2.715 0.007414 ** 
## SEMANA34                  0.242358   0.119862   2.022 0.044969 *  
## SEMANA35                  0.341818   0.118529   2.884 0.004511 ** 
## SEMANA36                  0.303472   0.129757   2.339 0.020676 *  
## SEMANA37                  0.318107   0.122775   2.591 0.010521 *  
## SEMANA38                  0.372854   0.122930   3.033 0.002857 ** 
## SEMANA39                  0.328007   0.123207   2.662 0.008615 ** 
## SEMANA40                  0.328143   0.127634   2.571 0.011121 *  
## SEMANA41                  0.208450   0.124516   1.674 0.096212 .  
## SEMANA42                  0.257197   0.115013   2.236 0.026822 *  
## SEMANA43                  0.239461   0.116181   2.061 0.041032 *  
## SEMANA44                  0.515072   0.118772   4.337 2.65e-05 ***
## SEMANA45                  0.222118   0.129668   1.713 0.088797 .  
## SEMANA46                  0.251365   0.115131   2.183 0.030579 *  
## SEMANA47                  0.246092   0.112833   2.181 0.030750 *  
## SEMANA48                  0.274417   0.116180   2.362 0.019469 *  
## SEMANA49                  0.342478   0.115395   2.968 0.003495 ** 
## SEMANA50                  0.245617   0.116066   2.116 0.035991 *  
## SEMANA51                  0.321246   0.116408   2.760 0.006513 ** 
## SEMANA52                  0.180701   0.123848   1.459 0.146655    
## SEMANA53                 -1.907555   0.107141 -17.804  < 2e-16 ***
## choque                    0.060202   0.019959   3.016 0.003009 ** 
## volcamiento               0.010564   0.008407   1.257 0.210888    
## herido                   -0.035709   0.019263  -1.854 0.065747 .  
## muerto                   -0.011589   0.008745  -1.325 0.187102    
## dias_festivos            -0.064942   0.016685  -3.892 0.000149 ***
## dias_especiales          -0.069392   0.023987  -2.893 0.004390 ** 
## dias_festivos_especiales  0.057739   0.027524   2.098 0.037615 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09639 on 149 degrees of freedom
## Multiple R-squared:  0.9441, Adjusted R-squared:  0.922 
## F-statistic: 42.69 on 59 and 149 DF,  p-value: < 2.2e-16
plot_estimation_errors <- function(y_true, y_pred, title){
  min_value <- min(y_true, y_pred) - min(y_true, y_pred) * 0.1
  max_value <- max(y_true, y_pred) + max(y_true, y_pred) * 0.1
  
  plot(x=y_true,y=y_pred,
       ylab="Predicciones",xlab="Observados",
       xlim=c(min_value, max_value),ylim=c(min_value, max_value),
       las=1, cex=1, pch=16,
       main=title)
  abline(a=0,b=1,lwd=2,col="black", lty=2)
  R_vl<-cor(y_pred, y_true)
  R_vl<-format(R_vl, digits = 3, nsmall = 3)
  rmse_vl<- rmse(actual = y_true, predicted = y_pred)
  rmse_vl<-format(rmse_vl,digits = 3,nsmall = 2)
  grid()
  legend("topleft",legend=paste0(c("Correlación: ","RMSE: "),c(R_vl,rmse_vl)), bty="n")
}
y_pred <- exp(predict(step.model.choque, df_choque_train))
y_true <- exp(df_choque_train$numero_accidentes)

train_rmse <- rmse(y_true, y_pred)

plot_estimation_errors(
  y_true = y_true, 
  y_pred = y_pred, 
  title="Rendimiento de modelo semanal choque ~ Entrenamiento"
  )

De primera instancia, se observa un muy buen ajuste de los datos de entrenamiento respecto al modelo de regresión lineal, sin embargo, para evaluar su capacidad de generalizar es necesario evaluarlo en el conjunto de validación.

y_pred <- exp(predict(step.model.choque, df_choque_test))
y_true <- exp(df_choque_test$numero_accidentes)

test_rmse <- rmse(y_true, y_pred)

test_train_ratio <- round(((test_rmse / train_rmse) - 1)*100, 3)

plot_estimation_errors(
  y_true = y_true, 
  y_pred = y_pred, 
  title="Rendimiento de modelo mensual semanal ~ Validación"
  )

Para el conjunto de validación, se puede ver un buen ajuste del modelo, sin embargo, respecto al rmse de entrenamiento se tiene una variación del 33.066% lo cual genera indicios de sobrentrenamiento. Procedemos a comparar directamente las funciones de distribución de probabilidad de los errores de entrenamiento vs los de validación.

plot_errors_density <- function(model, train, y_train, test, y_test, title){
  # data extraction
  y_pred <- exp(predict(model, train))
  y_true <- exp(y_train)
  train_error <- y_true - y_pred
  y_pred <- exp(predict(model, test))
  y_true <- exp(y_test)
  test_error <- y_true - y_pred
  #errors dataframe
  errors_df <- data.frame(
    errors = c(train_error, test_error),
    error_type = c(
      rep("train", length(train_error)),
      rep("test", length(test_error))
    )
  )
  # plotting
  ggplot(errors_df, aes(x=errors, color=error_type)) +
    geom_density() + 
    theme_classic() + 
    ggtitle(title) + 
    theme(plot.title = element_text(hjust = 0.5))
}

plot_errors_density(
  model = step.model.choque, 
  train = df_choque_train, 
  y_train = df_choque_train$numero_accidentes,
  test = df_choque_test, 
  y_test = df_choque_test$numero_accidentes,
  title="Errores de entrenamiento y validación"
  )

No se ven diferencias significativas, lo que si es claro es que el modelo en el conjunto de validación tiende a sobreestimar algunas observaciones por lo que se ve que la media no está centrada en cero sino que se ve levemente desplazada hacia la izquierda al igual que la función de distribución de probabilidad una mayor concentración de la masa a la izquierda.
Debido a esto, procedemos a utilizar métodos de regularización para realizar una selección de variables más adecuada.

Regresión Lasso

train_lasso_model <- function(df_train, df_test, min_lambda=-2){
  # Construcción del dataset
  x_train <- model.matrix(numero_accidentes~. , df_train)[,-1]
  y_train <- df_train$numero_accidentes
  
  # malla de lambdas a probar
  lambda_seq <- 10^seq(2, min_lambda, by = -.1)
  
  # modelado
  set.seed(42)
  cv.lasso <- cv.glmnet(x_train, y_train,
                        alpha = 1,
                        lambda = lambda_seq, 
                        nfolds = 5)
  
  lasso.model <- glmnet(x_train, y_train, alpha = 1, lambda = cv.lasso$lambda.min)  
  
  y_train_pred <- predict(lasso.model, x_train)
  train_rmse <- rmse(exp(y_train), exp(y_train_pred))
  
  
  x_test <- model.matrix(numero_accidentes~. , df_test)[,-1]
  y_test <- df_test$numero_accidentes
  y_test_pred <- predict(lasso.model, x_test)
  
  test_rmse <- rmse(exp(y_test), exp(y_test_pred))
  
  response=list(
    "training"=list(
      "data"=x_train,
      "response"=y_train,
      "predictions"=y_train_pred,
      "rmse"=train_rmse
    ),
    "validation"=list(
      "data"=x_test,
      "response"=y_test,
      "predictions"=y_test_pred,
      "rmse"=test_rmse
    ),
    "model"=lasso.model
  )
  return(response)
}
lasso.choque.results <- train_lasso_model(df_train = df_choque_train, df_test = df_choque_test)

plot_estimation_errors(
  y_true = exp(lasso.choque.results$training$response), 
  y_pred = exp(lasso.choque.results$training$predictions), 
  title="Rendimiento de modelo mensual choque ~ Entrenamiento"
  )

plot_estimation_errors(
  y_true = exp(lasso.choque.results$validation$response), 
  y_pred = exp(lasso.choque.results$validation$predictions), 
  title="Rendimiento de modelo mensual choque ~ Validación"
  )

test_train_ratio <- round((lasso.choque.results$validation$rmse / lasso.choque.results$training$rmse -1 )*100, 3)

Para el modelo de lasso en la clase choque, la variación del rmse de validación respecto a entrenamiento es de un 14.202 lo cual es cercano al 15%, por lo que procedemos a comparar las fpd de los errores de entrenamiento y validación.

plot_errors_density(
  model = lasso.choque.results$model, 
  train = lasso.choque.results$training$data, 
  y_train = lasso.choque.results$training$response,
  test = lasso.choque.results$validation$data, 
  y_test = lasso.choque.results$validation$response,
  title="Errores de entrenamiento y validación"
  )

Se puede evidenciar que ambas distribuciones son muy similares, sin embargo una leve diferencia en las medias. De esto, nos queda que el modelo más apropiado es el de lasso para la clase choque.

Clase atropello

Regresión Lasso

df_preprocessed <- armado_dataset_semanal(df=df, objective_var = "atropello")
df_atropello_train <- df_preprocessed %>% filter(PERIODO!=2018) %>% dplyr::select(-c(PERIODO))
df_atropello_test <- df_preprocessed %>% filter(PERIODO==2018) %>% dplyr::select(-c(PERIODO))

model.atropello.recipe <- recipe(numero_accidentes ~ ., 
                              data = df_atropello_train)
model.atropello.steps <- model.atropello.recipe %>% 
  step_log(all_outcomes()) %>% 
  step_center(all_predictors(), -SEMANA) %>% 
  step_scale(all_predictors(), -SEMANA)

model.atropello.prepared <- prep(model.atropello.steps, training = df_atropello_train)
df_atropello_train <- bake(model.atropello.prepared, df_atropello_train)
df_atropello_test <- bake(model.atropello.prepared, df_atropello_test)
lasso.atropello.results <- train_lasso_model(df_train = df_atropello_train, df_test = df_atropello_test, min_lambda = -3)

plot_estimation_errors(
  y_true = exp(lasso.atropello.results$training$response), 
  y_pred = exp(lasso.atropello.results$training$predictions), 
  title="Rendimiento de modelo mensual atropello ~ Entrenamiento"
  )

plot_estimation_errors(
  y_true = exp(lasso.atropello.results$validation$response), 
  y_pred = exp(lasso.atropello.results$validation$predictions), 
  title="Rendimiento de modelo mensual atropello ~ Validación"
  )

test_train_ratio <- round((lasso.atropello.results$validation$rmse / lasso.atropello.results$training$rmse -1 )*100, 3) 
plot_errors_density(
  model = lasso.atropello.results$model, 
  train = lasso.atropello.results$training$data, 
  y_train = lasso.atropello.results$training$response,
  test = lasso.atropello.results$validation$data, 
  y_test = lasso.atropello.results$validation$response,
  title="Errores de entrenamiento y validación"
  )

Para el modelo de atropellos, se ve una deficiencia en la calidad de las predicciones, en especial para el conjutno de validación, donde es claro un sobreestimación de las observaciones, sin embargo, en terminos de RMSE, es un modelo que no posee una gran variación respecto al conjunto de entrenamiento.

Clase caida ocupante

Regresión Lasso

df_preprocessed <- armado_dataset_semanal(df=df, objective_var = "caida ocupante")
df_caida_ocupante_train <- df_preprocessed %>% filter(PERIODO!=2018) %>% dplyr::select(-c(PERIODO))
df_caida_ocupante_test <- df_preprocessed %>% filter(PERIODO==2018) %>% dplyr::select(-c(PERIODO))

model.caida_ocupante.recipe <- recipe(numero_accidentes ~ ., 
                              data = df_caida_ocupante_train)
model.caida_ocupante.steps <- model.caida_ocupante.recipe %>% 
  step_log(all_outcomes()) %>% 
  step_center(all_predictors(), -SEMANA) %>% 
  step_scale(all_predictors(), -SEMANA)

model.caida_ocupante.prepared <- prep(model.caida_ocupante.steps, training = df_caida_ocupante_train)
df_caida_ocupante_train <- bake(model.caida_ocupante.prepared, df_caida_ocupante_train)
df_caida_ocupante_test <- bake(model.caida_ocupante.prepared, df_caida_ocupante_test)
lasso.caida_ocupante.results <- train_lasso_model(df_train = df_caida_ocupante_train, df_test = df_caida_ocupante_test, min_lambda = -3)

plot_estimation_errors(
  y_true = exp(lasso.caida_ocupante.results$training$response), 
  y_pred = exp(lasso.caida_ocupante.results$training$predictions), 
  title="Rendimiento de modelo mensual caida_ocupante ~ Entrenamiento"
  )

plot_estimation_errors(
  y_true = exp(lasso.caida_ocupante.results$validation$response), 
  y_pred = exp(lasso.caida_ocupante.results$validation$predictions), 
  title="Rendimiento de modelo mensual caida_ocupante ~ Validación"
  )

test_train_ratio <- round((lasso.caida_ocupante.results$validation$rmse / lasso.caida_ocupante.results$training$rmse -1 )*100, 3) 
plot_errors_density(
  model = lasso.caida_ocupante.results$model, 
  train = lasso.caida_ocupante.results$training$data, 
  y_train = lasso.caida_ocupante.results$training$response,
  test = lasso.caida_ocupante.results$validation$data, 
  y_test = lasso.caida_ocupante.results$validation$response,
  title="Errores de entrenamiento y validación"
  )

Para caida ocupante, la variación del RMSE de validación respecto a entrenamiento es de un 11.725% lo cual es un buen indicativo de la capacidad de generalización del modelo. Adicional, las FDP de los errores de entrenamiento y validación no se diferencian mucho por lo que se ve que el modelo tiene un buen ajuste.

Clase volcamiento

Regresión Lasso

df_preprocessed <- armado_dataset_semanal(df=df, objective_var = "volcamiento")
df_volcamiento_train <- df_preprocessed %>% filter(PERIODO!=2018) %>% dplyr::select(-c(PERIODO))
df_volcamiento_test <- df_preprocessed %>% filter(PERIODO==2018) %>% dplyr::select(-c(PERIODO))

model.volcamiento.recipe <- recipe(numero_accidentes ~ ., 
                              data = df_volcamiento_train)
model.volcamiento.steps <- model.volcamiento.recipe %>% 
  step_log(all_outcomes()) %>% 
  step_center(all_predictors(), -SEMANA) %>% 
  step_scale(all_predictors(), -SEMANA)

model.volcamiento.prepared <- prep(model.volcamiento.steps, training = df_volcamiento_train)
df_volcamiento_train <- bake(model.volcamiento.prepared, df_volcamiento_train)
df_volcamiento_test <- bake(model.volcamiento.prepared, df_volcamiento_test)
lasso.volcamiento.results <- train_lasso_model(df_train = df_volcamiento_train, df_test = df_volcamiento_test, min_lambda = -2)

plot_estimation_errors(
  y_true = exp(lasso.volcamiento.results$training$response), 
  y_pred = exp(lasso.volcamiento.results$training$predictions), 
  title="Rendimiento de modelo mensual volcamiento ~ Entrenamiento"
  )

plot_estimation_errors(
  y_true = exp(lasso.volcamiento.results$validation$response), 
  y_pred = exp(lasso.volcamiento.results$validation$predictions), 
  title="Rendimiento de modelo mensual volcamiento ~ Validación"
  )

test_train_ratio <- round((lasso.volcamiento.results$validation$rmse / lasso.volcamiento.results$training$rmse -1 )*100, 3) 
plot_errors_density(
  model = lasso.volcamiento.results$model, 
  train = lasso.volcamiento.results$training$data, 
  y_train = lasso.volcamiento.results$training$response,
  test = lasso.volcamiento.results$validation$data, 
  y_test = lasso.volcamiento.results$validation$response,
  title="Errores de entrenamiento y validación"
  )

A pesar de que el modelo de volcamiento no posea una variación tan alta -2.724%, es evidente que el modelo se encuentra subestimando los datos, y el rango de las predicciones no es acorde al rango lo la variabilidad natural del fenomeno en estudio, es necesario profundizar para este modelo.

Clase otros

Regresión Lasso

df_preprocessed <- armado_dataset_semanal(df=df, objective_var = "otro")
df_otro_train <- df_preprocessed %>% filter(PERIODO!=2018) %>% dplyr::select(-c(PERIODO))
df_otro_test <- df_preprocessed %>% filter(PERIODO==2018) %>% dplyr::select(-c(PERIODO))

model.otro.recipe <- recipe(numero_accidentes ~ ., 
                              data = df_otro_train)
model.otro.steps <- model.otro.recipe %>% 
  step_log(all_outcomes()) %>% 
  step_center(all_predictors(), -SEMANA) %>% 
  step_scale(all_predictors(), -SEMANA)

model.otro.prepared <- prep(model.otro.steps, training = df_otro_train)
df_otro_train <- bake(model.otro.prepared, df_otro_train)
df_otro_test <- bake(model.otro.prepared, df_otro_test)
modelo.otro <- lm(numero_accidentes ~ ., data=df_otro_train)
step.model.otro <- stepAIC(modelo.otro, direction = "both", trace = F, steps = 2000)

model.otro.summary <- summary(step.model.otro)
model.otro.summary
## 
## Call:
## lm(formula = numero_accidentes ~ SEMANA + incendio + otro + volcamiento + 
##     herido + t_minus_2 + t_minus_3 + dias_festivos + dias_especiales + 
##     dias_festivos_especiales, data = df_otro_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.55054 -0.08653 -0.00171  0.08862  0.69727 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               4.27814    0.17600  24.308  < 2e-16 ***
## SEMANA2                   0.27741    0.18098   1.533 0.127461    
## SEMANA3                   0.39366    0.20096   1.959 0.052013 .  
## SEMANA4                   0.30154    0.19419   1.553 0.122624    
## SEMANA5                   0.25222    0.19126   1.319 0.189323    
## SEMANA6                   0.51760    0.19573   2.645 0.009069 ** 
## SEMANA7                   0.40548    0.21382   1.896 0.059873 .  
## SEMANA8                   0.35929    0.19863   1.809 0.072520 .  
## SEMANA9                   0.30940    0.20533   1.507 0.133998    
## SEMANA10                  0.53931    0.20793   2.594 0.010456 *  
## SEMANA11                  0.24996    0.20918   1.195 0.234034    
## SEMANA12                  0.48539    0.21875   2.219 0.028025 *  
## SEMANA13                  0.21568    0.19019   1.134 0.258631    
## SEMANA14                  0.16919    0.19322   0.876 0.382673    
## SEMANA15                  0.29590    0.18722   1.581 0.116141    
## SEMANA16                  0.19678    0.18842   1.044 0.298027    
## SEMANA17                  0.24500    0.19287   1.270 0.205988    
## SEMANA18                  0.34719    0.19626   1.769 0.078968 .  
## SEMANA19                  0.43957    0.20284   2.167 0.031843 *  
## SEMANA20                  0.33075    0.20083   1.647 0.101710    
## SEMANA21                  0.48797    0.19694   2.478 0.014352 *  
## SEMANA22                  0.58663    0.21010   2.792 0.005931 ** 
## SEMANA23                  0.32410    0.20649   1.570 0.118670    
## SEMANA24                  0.35902    0.20254   1.773 0.078369 .  
## SEMANA25                  0.27361    0.20370   1.343 0.181290    
## SEMANA26                  0.18691    0.19944   0.937 0.350209    
## SEMANA27                  0.26044    0.18444   1.412 0.160042    
## SEMANA28                  0.44733    0.19392   2.307 0.022462 *  
## SEMANA29                  0.34981    0.19725   1.773 0.078223 .  
## SEMANA30                  0.26717    0.19965   1.338 0.182908    
## SEMANA31                  0.52721    0.20121   2.620 0.009711 ** 
## SEMANA32                  0.17184    0.22679   0.758 0.449839    
## SEMANA33                  0.37244    0.20675   1.801 0.073694 .  
## SEMANA34                  0.35525    0.20577   1.726 0.086374 .  
## SEMANA35                  0.45314    0.20266   2.236 0.026858 *  
## SEMANA36                  0.29861    0.22246   1.342 0.181577    
## SEMANA37                  0.34551    0.20510   1.685 0.094198 .  
## SEMANA38                  0.34223    0.21014   1.629 0.105547    
## SEMANA39                  0.40422    0.20613   1.961 0.051768 .  
## SEMANA40                  0.31852    0.21830   1.459 0.146672    
## SEMANA41                  0.30822    0.20775   1.484 0.140064    
## SEMANA42                  0.30943    0.19507   1.586 0.114841    
## SEMANA43                  0.33303    0.19601   1.699 0.091419 .  
## SEMANA44                  0.55917    0.20088   2.784 0.006082 ** 
## SEMANA45                  0.25143    0.21831   1.152 0.251306    
## SEMANA46                  0.49746    0.19567   2.542 0.012045 *  
## SEMANA47                  0.07649    0.19573   0.391 0.696528    
## SEMANA48                  0.14629    0.19553   0.748 0.455536    
## SEMANA49                  0.20191    0.18869   1.070 0.286344    
## SEMANA50                  0.23886    0.18985   1.258 0.210331    
## SEMANA51                  0.35023    0.18725   1.870 0.063419 .  
## SEMANA52                  0.30142    0.19945   1.511 0.132857    
## SEMANA53                 -1.61770    0.18143  -8.916 1.74e-15 ***
## incendio                  0.01935    0.01379   1.403 0.162665    
## otro                      0.04265    0.02007   2.125 0.035261 *  
## volcamiento               0.03446    0.01431   2.409 0.017245 *  
## herido                   -0.04758    0.02823  -1.685 0.094068 .  
## t_minus_2                 0.02649    0.01840   1.440 0.152097    
## t_minus_3                 0.03968    0.01835   2.163 0.032190 *  
## dias_festivos            -0.05970    0.02935  -2.034 0.043722 *  
## dias_especiales          -0.11317    0.04211  -2.688 0.008023 ** 
## dias_festivos_especiales  0.16453    0.04819   3.414 0.000827 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1687 on 147 degrees of freedom
## Multiple R-squared:  0.8287, Adjusted R-squared:  0.7576 
## F-statistic: 11.65 on 61 and 147 DF,  p-value: < 2.2e-16
plot_estimation_errors(
  y_true = exp(df_otro_train$numero_accidentes), 
  y_pred = exp(predict(step.model.otro, df_otro_train)), 
  title="Rendimiento de modelo mensual otro ~ Entrenamiento"
  )

plot_estimation_errors(
  y_true = exp(df_otro_test$numero_accidentes), 
  y_pred = exp(predict(step.model.otro, df_otro_test)), 
  title="Rendimiento de modelo mensual otro ~ Validación"
  )

plot_errors_density(
  model = step.model.otro, 
  train = df_otro_train, 
  y_train = df_otro_train$numero_accidentes,
  test = df_otro_test, 
  y_test = df_otro_test$numero_accidentes,
  title="Errores de entrenamiento y validación"
  )

Claramente el modelo de otros con una regresión lineal muestra indicios de sobreentrenamiento, vamos a intentar con una regresión lasso.

lasso.otro.results <- train_lasso_model(df_train = df_otro_train, df_test = df_otro_test, min_lambda = -2)

plot_estimation_errors(
  y_true = exp(lasso.otro.results$training$response), 
  y_pred = exp(lasso.otro.results$training$predictions), 
  title="Rendimiento de modelo mensual otro ~ Entrenamiento"
  )

plot_estimation_errors(
  y_true = exp(lasso.otro.results$validation$response), 
  y_pred = exp(lasso.otro.results$validation$predictions), 
  title="Rendimiento de modelo mensual otro ~ Validación"
  )

test_train_ratio <- round((lasso.otro.results$validation$rmse / lasso.otro.results$training$rmse -1 )*100, 3) 
plot_errors_density(
  model = lasso.otro.results$model, 
  train = lasso.otro.results$training$data, 
  y_train = lasso.otro.results$training$response,
  test = lasso.otro.results$validation$data, 
  y_test = lasso.otro.results$validation$response,
  title="Errores de entrenamiento y validación"
  )

Aunque la variación del RMSE de entrenamiento a validación para la regresión lasso es menor que la regresión lineal, claramente el modelo está sobreestimando varias de las observaciones.
Se hace necesario evaluar nuevos modelos para está clase de accidentes.

Conclusiones

Los modelos semanales tuvieron un buen desempeño a nivel general, a excepeción de los modelos de volcamiento y otros, los cuales sufrieron de tener rangos muy cerrados de estimación, por lo que no tienen en cuenta toda la variabilidad total del fenomeno de estudio. Es importante explorar nuevas variables que le permitan a los modelos capturar la variabilidad de estás clases.