Modelado
Proyecto Final
Alfredo Lefranc Flores, Cynthia Raquel Valdivia Tirado, Rafael Sandoval Fernández, Marco Antonio Ramos Juárez y Francisco Velazquez Guadarrama
2021
Ambiente
El primer paso para el análisis partiendo de modelos es importar los datos procesados de la etapa anterior.Por ello, cargamos la base y revisamos sus variables. Notamos outliers en la variable de duración. Los truncamos al 99.9% para que no afecten las predicciones de nuestros modelos.
Para poder aumentar la eficiencia de los modelos, transformamos la base a formato ancho. Esto lo hacemos para dos variables con muchas categorías, las estaciones de salida y llegada.
# formato wide de start y end station
stations_wide <- trips_filter %>% select("seq_id", "strt_statn",
"end_statn") %>% as.data.table() %>% one_hot()
trips_filter <- inner_join(trips_filter, stations_wide, by = "seq_id")
trips_filter <- trips_filter %>% select(-c("strt_statn", "end_statn")) %>%
arrange(seq_id)Modelos
Primero separamos la base de validación y de entrenamiento. Usamos una división 80-20. Para esto decidimos estratificar por días y no por viajes. De esta manera tenemos un panorama completo de los días en entrenamiento, y evaluaremos en días nuevos en la base de validación. Esto nos parece un ejercicio más realista dados los objetivos del negocio.
Para tener un balance en nuestros sets de entrenamiento y validación, estratificamos por año, mes y día de la semana. Una vez hecha la estratificación, verificamos el balance entre entrenamiento y validación.
# para ello, nos fijamos en todos los días donde tenemos
# viajes
dates <- trips_filter %>% select(c("start_date", "start_month",
"start_weekday", "start_year")) %>% unique()
# estratificamos
treat_id <- treatment_assign(data = dates, share_control = 0.2,
n_t = 1, strata_varlist = c("start_month", "start_weekday",
"start_year"), missfits = "global", seed = 2021, key = "start_date")$data
# pegamos la variable que asigna tratamiento a la base
trips_filter <- left_join(trips_filter, treat_id, by = ("start_date"))
trips_filter <- trips_filter %>% select(-c("strata", "missfit"))Al hacer la separación, nos cercioramos de que exista balance entre ambas bases. Esto lo hacemos con prubas t y F para las variables numéricas.
# summary(trips_filter$treat) el share control se mantiene
# estable
# revisamos balance de clases no estratificadas
# pruebas t en las numéricas
is_num <- lapply(trips_filter, function(x) is.numeric(x))
trips_num <- trips_filter %>% select((names(is_num[is_num ==
TRUE])))
balance_t <- balance_table(data = sample_n(trips_num, 10000),
treatment = "treat")
balance_t$p_value1 <- balance_t$p_value1 %>% round(3)
# filtramos las variables donde se rechaza la hipótesis de
# igualdad al 10%
balance_t_filter <- balance_t %>% filter(p_value1 < 0.1)
datatable(balance_t_filter, options = list(pageLength = 5))37 variables muestran desbalance, dummies de estaciones y edad. La significancia conjunta no se rechaza. Por ende, consideramos que no es preocupante.
# pruebas F
balance_f <- balance_regression(data = sample_n(trips_num, 100000),
treatment = "treat")
kable(balance_f$F_test, align = "c", booktabs = T, digits = 2,
caption = "Regresión de balance", col.names = c("Estadístico",
"Valor")) %>% kable_styling(position = "center")| Estadístico | Valor |
|---|---|
| F-statistic | 1.19 |
| k | 291.00 |
| n-k-1 | 99708.00 |
| F_critical | 0.87 |
| p_value | 0.01 |
| R cuadrada | 0.00 |
Ahora procedemos a separar las bases en entrenamiento y validación. Acotamos la de entrenamiento para tener una muestra más pequeña. Convertimos las bases de covariables en matrices ralas.
train <- trips_filter[trips_filter$treat == 1, ]
y_train <- train$duration[train$treat == 1]
train <- train %>% select(-c("treat", "seq_id", "duration"))
sparse_train <- sparse.model.matrix(~. + 0, data = train)
# validacion
test <- trips_filter[trips_filter$treat == 0, ] %>% select(-c("treat",
"seq_id", "duration"))
sparse_test <- sparse.model.matrix(~. + 0, data = test)
y_test <- trips_filter[trips_filter$treat == 0, ] %>% select(c("seq_id",
"duration"))Lasso
Primero probamos un LASSO con Cross Validation.
lasso <- cv.gamlr(x = sparse_train, y = y_train, verb = T, nfold = 10)
seg100
5.554548 Guardamos las predicciones del modelo.
Random Forest
El segundo modelo que deseamos probar es el random forest. Probamos con distintas B, desde 100 hasta 700, siendo 500 el óptimo por reducir el error OOB.
Por razones de tiempo, corrimos el modelo sobre una submuestra de entrenamiento. Nos aseguramos que ésta respetara la estratificación antes descrita.
# ajuste del set de entrenamiento
dates <- train %>% select(c("start_date", "start_month", "start_weekday",
"start_year")) %>% unique()
# estratificamos
treat_id <- treatment_assign(data = dates, share_control = 0.75,
n_t = 1, strata_varlist = c("start_month", "start_weekday",
"start_year"), missfits = "global", seed = 2021, key = "start_date")$data
# pegamos la variable que asigna tratamiento a la base
train_cf <- left_join(cbind(train, y_train), treat_id, by = ("start_date"))
train_cf <- train_cf %>% select(-c("strata", "missfit"))
# subset de la base para los días elegidos
train_cf <- train_cf[train_cf$treat == 1, ]
y_train_cf <- train_cf$y_train[train_cf$treat == 1]
train_cf <- train_cf %>% select(-c("treat", "y_train"))
sparse_train_cf <- sparse.model.matrix(~. + 0, data = train_cf)Una vez hechos los ajustes, corremos el RF.
# para correr modelos en paralelo
cores <- detectCores()
cl <- makeCluster(cores)
inicio <- Sys.time()
# c(100,200,350,500,700)
random_forest <- map(c(500), function(z) ranger(y = y_train_cf,
x = sparse_train_cf, num.trees = z, mtry = ncol(train) %>%
sqrt() %>% floor(), min.node.size = 1, importance = "impurity",
status.variable.name = 1))
(tiempo <- Sys.time() - inicio)
stopCluster(cl)
error_prediccion <- tibble(trees = c(500), oob_error = map_dbl(random_forest,
~.$prediction.error))
# error medido en horas ggplot(error_prediccion, aes(trees,
# oob_error/60)) + geom_point() + geom_path() + theme_bw()df_imp <- data.frame(names = random_forest[[1]][["variable.importance"]] %>%
names(), importance = random_forest[[1]][["variable.importance"]]) %>%
arrange(-importance)
# predicciones OOS
y_test$pred_rf <- predict(random_forest[[1]], data = sparse_test)$predictionsPredicting.. Progress: 53%. Estimated remaining time: 27 seconds.y_test <- y_test %>% mutate(oos_rf = (duration - pred_rf)^2)XGB
Finalmente, probamos un XGB. Entrenamos el modelo y unimos las predicciones a la base de evaluación.
# Preparar la base de entrenamiento y de validación
dtrain <- xgb.DMatrix(sparse_train, label = y_train)
dtest <- xgb.DMatrix(sparse_test, label = y_test$duration)
watchlist <- list(train = dtrain, eval = dtest)
# Entrenamiento del modelo
param <- list(max_depth = 5, learning_rate = 0.06, objective = "reg:squarederror",
eval_metric = "rmse", subsample = 0.85, colsample_bytree = 0.7)
xgb_model <- xgb.train(params = param, dtrain, early_stopping_rounds = 10,
nrounds = 100, watchlist)# Predicción
y_test$pred_xgb <- predict(xgb_model, sparse_test)
y_test <- y_test %>% mutate(oos_xgb = (duration - pred_xgb)^2)Desempeño de modelos
Medidas tradicionales
Como medidas tradicionales, podemos compara ambos modelos midiendo el tamaño de los residuales. Al ser la variable objetico numérica, podemos emplear el R^2.
| \[R^2\] | Minutos residuales totales | Desviación en minutos | |
|---|---|---|---|
| RF | 0.116 | 22013847327 | 0.179 |
| Lasso | 0.073 | 23082505149 | 0.344 |
| XGB | 0.134 | 21568657965 | 0.142 |
Medidas alternativas
Predicciones y flujo de biciletas por estación
Para mostrar el desempeño de los modelos, modelaremos los flujos de bicicleta de nuestra base de validación pero ahora en lugar de usar la duración del viaje real, usaremos la predicha. Para ellos crearemos primero una base de datos que compare los flujos de viaje cada 10 minutos entre lo que en realidad pasó y lo que predice nuestro modelo.
En la siguiente tabla se muestran los flujos acumulados de bicicletas que salen y llegan (con datos reales y predecidos). Por motivos de costo computacional, solamente mostramos algunos días completos para algunas estaciones aleatorias. Por favor haga click primer primero en estación y luego en fecha para que la tabla esté ordenada de acuerdo a la hora:
# Deployment de las bases de datos en una app (selecciono de
# manera aleatoria algunas debido al peso computacional)
random_stations <- sample(3:133, 3, replace = FALSE)
output <- plantilla_llena %>% filter(strt_statn %in% random_stations)
# output<-sample_n(plantilla_llena, 1000)
output <- output %>% select(date, hour, strt_statn, acum_salida,
acum_llegada_real, acum_llegada_pred_rf, acum_llegada_pred_lasso,
acum_llegada_pred_xgb)
names(output) <- c("Fecha", "Hora", "Estación", "Salidas", "Llegadas(real)",
"RF", "Lasso", "XGB")
datatable(output, options = list(pageLength = 15))Error promedio a la hora de corte
Derivado de estas tablas, se nos ocurrió que una forma muy util de probar el desempeño de los modelos es a través del error del flujo en determinado punto del día. Esto pensando en que en las empresas de este tipo, hay un punto de corte en el que las estaciones se tienen que llenar o vaciar de bicicletas. Para este ejemplo, propusimos un punto de corte a las 12 de la noche, pero bien podría ser en la madrugada, o en la tarde. En la siguiente tabla se muestra el error promedio en el cálculo del flujo de bicis para la hora de corte.
Conclusiones
Con estos resultados es posible optimizar la gestión del servicio de bicicletas otorgado por la empresa Wheelie Wonka. Este pronóstico preciso sirve de guía para que los pasajeros puedan organizar mejor su hora de salida y las rutas de desplazamiento. Además, es beneficioso para el proveedor de servicios de bicicletas compartidas Wheelie Wonka en términos de mejorar la satisfacción de sus clientes y organizar con efectividad el horario de entrega y distribución de bicicletas y así eficientar sus recursos y disminuir pérdidas.
De manera más técnica, el modelo que resultó dar los mejores resultados tanto desde una perspectiva de análisis de resiudales tradicional como desde un problema aplicado de negocios fue el XGB.