library(dplyr)
library(rmarkdown)
library(Metrics)
library(plotly)
library(htmltools)Para realizar la predicción se hace uso del conjunto de datos cuyos valores de bicicletas disponibles es la agregación por horas:
Ademas, a partir de esta base de datos, se obtiene el número total de bornes de cada estación (lo que, cómo se verá, será necesario para realizar la predicción):
df_bicis_total <- df_merge %>% select('number_','avg_total')
df_bicis_total <- unique(df_bicis_total)
df_bicis_total <- df_bicis_total[!(df_bicis_total$number_ == 255 & df_bicis_total$avg_total == 13.5),]Como se analizará posteriormente, el orden de las filas es clave para que el programa cree una buena matriz de datos a la cual aplicar el modelo de regresión. En consecuencia, posterior conversión de la variable “fecha” a objetos fecha mediante la función “as.Date()”, se ordenó la base de datos por fecha y hora.
df_merge$fecha <- as.Date(df_merge$fecha, "%d/%m/%Y")
df_merge <- df_merge[order(df_merge$fecha, df_merge$hora_hora),]Una vez ordenados los datos, para la realización de este objetivo, hay que seleccionar únicamente las variables que son necesarias: el identificador de la estación, el número de bicis disponibles, el día de la semana, la hora y la fecha.
df_merge <- df_merge %>% select('number_','avg_av', 'día_semana', 'hora_hora', 'fecha')
colnames(df_merge) <- c('id_estacion', 'bicis_disponibles', 'dia_semana', 'hora', 'fecha')Posteriormente, la variable “hora” se convierte en una variable de tipo hora mediante la función “as.POSIXlt()” y, una vez hecho esto, se convierte a factor. Por su parte, las variables “día_semana” e “id_estación” también se convierten a factores, Finalmente, la variable “bicis_disponibles” se convierte a variable numérica y, cómo consecuencia de la comentada agrupación por horas, se trunca.
df_merge$hora <- as.POSIXlt(df_merge$hora, format="%H:%M:%S")$hour
df_merge$hora <- as.factor(df_merge$hora)
df_merge$dia_semana <- as.factor(df_merge$dia_semana)
df_merge$id_estacion <- as.factor(df_merge$id_estacion)
df_merge$bicis_disponibles <- as.numeric(df_merge$bicis_disponibles)
df_merge$bicis_disponibles <- trunc(df_merge$bicis_disponibles)Adicionalmente, se crea el conjunto de datos test con datos desde el 9 de febrero de 2023 (inclusive) hasta el 23 de febrero de 2023 (inclusive). Este conjunto no se utilizará para realizar la predicción (para ello se construirá una matriz) pero sí que será donde se almacenen las predicciones.
fecha_inicial_test <- as.Date("2023-02-09")
fecha_final_test <- as.Date("2023-02-23")
test_df <- df_merge[df_merge$fecha >= fecha_inicial_test & df_merge$fecha <= fecha_final_test, ]Para intentar predecir el número de bicis, se parte de una estructura de datos como la vista anteriormente. Sin embargo, para realizar la predicción, se decidió, después de probar otras variables con errores mayores y coeficientes de determinación menores, escoger como variables predictoras las bicis disponibles en las 3 horas anteriores y en la misma hora que se intenta predecir pero en el mismo día de la semana anterior (por ejemplo, si se quiere predecir el número de estaciones en un estación a las 14 horas, las variables predictoras son el número de bicis disponibles a las 13 horas, a las 12 horas, a las 11 horas y el número de bicis a las 14 horas del mismo día de la semana pero en la semana anterior)
Para pasar de una estructura de datos como vista anteriormente a la matriz de datos necesaria para realizar la predicción se utilizó la siguiente función, que recibe como parámetros la base de datos, un id de estación y una hora:
get_matrix <- function(df, id_estacion, hora){
datos_estacion <- df[df$id_estacion == id_estacion, ]
index_hora <- which(datos_estacion$hora == hora)
index_res <- c()
for (i in index_hora){
if (i < 4){
next
} else{
v = c(i - 1:3, i)
index_res = c(index_res, sort(v))
}
}
index_res = unique(unlist(index_res))
filas_horas <- datos_estacion[index_res, ]
df_fechas <- filas_horas
filas_horas <- subset(filas_horas, fecha != "2022-12-12")
filas_horas <- filas_horas[,-c(1)]
id <- rep(1:ceiling(nrow(filas_horas)/4), each = 4, length.out = nrow(filas_horas))
filas_horas <- cbind(id, filas_horas)
nuevos_datos <- reshape(filas_horas, idvar = "id", timevar = "hora", direction = "wide")
n <- nrow(nuevos_datos)
semana_anterior <- rep(NA, n)
nuevos_datos <- cbind(semana_anterior, nuevos_datos)
for (i in 1:nrow(nuevos_datos)){
fecha_act <- nuevos_datos[, ncol(nuevos_datos)][i]
fecha_ant = fecha_act - 7
datos_ant = df_fechas[df_fechas$hora == hora & df_fechas$fecha == fecha_ant,]
bicis_ant = datos_ant$bicis_disponibles
if (is.numeric(bicis_ant) && length(bicis_ant) > 0) {
nuevos_datos$semana_anterior[i] <- bicis_ant
}
}
return (nuevos_datos)
}Esta función realiza el siguiente proceso:
En primer lugar, se seleccionan únicamente las instancias para las cuales el id de la estación sea el mismo que el que se pasa como parámetro.
A continuación, se seleccionan las filas cuyo valor en el atributo hora sea el mismo que el valor del parámetro hora. Al estar ordenadas, basta con obtener las 3 filas posteriores para obtener las filas correspondientes a las 3 horas anteriores a la hora dada. El hecho de estar ordenadas, además, elimina el problema que podría suponer obtener las 3 horas anteriores de, por ejemplo, las 00:00.
Sobre esta nueva base de datos, se realizan 3 acciones:
Utilizando la función reshape, se crea una matriz donde cada fila se corresponde con un valor distinto de id, y las columnas son para cada hora: la fecha, la hora, el día de la semana y las bicis disponibles. En la siguiente función, se seleccionarán únicamente las columnas correspondientes a las bicis.
Para obtener los datos del mismo día de la semana a la misma hora, pero de la semana anterior, en primer lugar, se crea una nueva columna con valores faltantes. Una vez hecho esto, para cada fila de la matriz anterior:
Se calcula la fecha del mismo día de la semana, pero de la semana anterior, es decir, fecha – 7.
Con esto, se obtiene la cantidad de bicis disponibles en la fecha anterior a la hora que se pasa como parámetro. Si el valor es distinto de NA, se añade a la matriz
Pasando como parámetros la base de datos, id_estacion = 1 y hora = 0, el resultado de la función, posterior eliminación de algunas variables que no se considerarán, es el siguiente:
Una vez creada la matriz, el siguiente paso es realizar las predicciones, para lo cual se creó otra función, que recibe como parámetros el id de una estación y una hora determinada:
predict_bikes_test <- function(id_estacion, hora){
df = get_matrix(df_merge, id_estacion, hora)
total_estacion <- df_bicis_total[df_bicis_total$number_ == id_estacion,]$avg_total
fecha_inicial_train <- as.Date("2022-12-02")
fecha_final_test <- as.Date("2023-02-09")
train_matrix <- df[df[, ncol(df)] >= fecha_inicial_train & df[, ncol(df)] < fecha_final_test, ]
train_matrix <- train_matrix[,-c(2,4,5,7,8,10,11,13,14)]
names(train_matrix) = c("semana_anterior", "hora_3ant", 'hora_2ant','hora_ant','hora_pred')
model <- lm(hora_pred ~., data = train_matrix)
fecha_inicial_test <- as.Date("2023-02-09")
fecha_final_test <- as.Date("2023-02-23")
test_matrix <- df[df[, ncol(df)] >= fecha_inicial_test & df[, ncol(df)] <= fecha_final_test, ]
aux_indices_filas <- which(test_df$id_estacion == id_estacion & test_df$hora == hora)
test_matrix <- test_matrix[,-c(2,4,5,7,8,10,11,13,14)]
names(test_matrix) = c("semana_anterior", "hora_3ant", 'hora_2ant','hora_ant','hora_pred')
bicis = predict(model, newdata = test_matrix)
for(i in 1:length(bicis)) {
if(bicis[i] < 0) {
bicis[i] <- 0
} else if (bicis[i] > total_estacion){
bicis[i] <- total_estacion
}
}
res = list('indices' = aux_indices_filas, 'predicho' = bicis)
return (res)
}Esta función realiza el siguiente proceso:
Llama a la función anterior, obteniendo así una matriz cuyos datos se corresponden con la estación que recibe como parámetro. Las bicis a predecir es el número de bicis en la dicha estación a la hora que se pasa como parámetro.
Se divide la matriz obtenida en dos matrices de datos:
Una matriz de train, con datos desde el 2 de diciembre de 2022 (inclusive) hasta el 9 de febrero de 2023 (no inclusive)
Una matriz de test con datos desde el 9 de febrero de 2023 (inclusive) hasta el 23 de febrero de 2023 (inclusive)
Los datos de train se utilizan para entrenar un modelo de regresión múltiple, donde la variable respuesta son los datos de la hora que se pasa cómo parámetro y las variables predictoras son el resto de las variables (las bicicletas disponibles en las 3 horas anteriores y en el mismo día de la semana a la misma hora, pero en la semana anterior)
Los datos del test se pasan a la función predict() a la cual también se le pasa el modelo generado. Esta función devuelve un vector de la misma longitud que el número de filas de la matriz de test. Para cada uno de los elementos del vector, se comprueba que el elemento no sea mayor al número de bicis máximo de la estación (si esto ocurre el valor se sustituye por el número de bornes de la estación) y que no tenga un valor negativo (si esto ocurre, el valor se sustituye por 0)
Finalmente, del conjunto de datos de testeo (no la matriz) se obtiene las filas tales que el id de estación y la hora sean los valores que se pasan como parámetros. Con esto, se tienen tanto los valores predichos, como los índices de las filas correspondientes a las predicciones. El resultado final de la función es una lista con los valores y las filas
Una vez definidas las dos funciones, se procede con la obtención de las predicciones de los datos de testeo. Para cada fila, se llama a la función predictora (pasándole como parámetros el id de la estación y la hora correspondiente) que devuelve, como se ha dicho un vector de valores y otro de índices. Con estos dos objetos, se asignan a la variable predicho de las filas en cuestión los valores correspondientes.
test_df_multiple_1_semanas <- test_df
test_df_multiple_1_semanas$predicho <- rep(NA, nrow(test_df_multiple_1_semanas))
for (i in 1:nrow(test_df_multiple_1_semanas)){
r = predict_bikes_test(test_df_multiple_1_semanas$id_estacion[i], test_df_multiple_1_semanas$hora[i])
if (is.na(test_df_multiple_1_semanas[r$indices[1], ]$predicho)){
filas <- test_df_multiple_1_semanas[r$indices, ]
filas$predicho <- r$predicho
test_df_multiple_1_semanas[r$indices, ] <- filas
} else if (sum(is.na(test_df_multiple_1_semanas$predicho)) == 0){
break
}
}Para comprobar la calidad de la predicción llevada a cabo (y por tanto de los 276 modelos generados), se han utilizado medidas como el rmse, el mae y el coeficiente de determinación. Sobre los conjuntos de test, dichos indicadores toman los siguientes valores:
rmse = rmse(test_df_multiple_1_semanas$bicis_disponibles, test_df_multiple_1_semanas$predicho)
rmse## [1] 1.766471mae = mae(test_df_multiple_1_semanas$bicis_disponibles, test_df_multiple_1_semanas$predicho)
mae## [1] 1.133816r2 = cor(test_df_multiple_1_semanas$bicis_disponibles, test_df_multiple_1_semanas$predicho)^2
r2## [1] 0.9170628En general, estos valores son considerablemente buenos. Sin embargo, para medir si realmente la predicción puede considerarse aceptable, decidimos crear un modelo que predijera mediante medias (la media de una estación, un día de la semana a una hora concreta).
fecha_inicial_train <- as.Date("2022-12-02")
fecha_final_test <- as.Date("2023-02-09")
train_df <- df_merge[df_merge$fecha >= fecha_inicial_train & df_merge$fecha < fecha_final_test, ]
predict_bikes_test <- function(id_estacion, dia_semana, hora){
datos_train = train_df[train_df$id_estacion == id_estacion & train_df$dia_semana == dia_semana & train_df$hora == hora,]
media_bicis <- mean(datos_train$bicis_disponibles, na.rm = TRUE)
filas_test = which(test_df$id_estacion == id_estacion & test_df$dia_semana == dia_semana & test_df$hora == hora)
res = list('indices' = filas_test, 'predicho' = media_bicis)
return(res)
}Los mismos indicadores para este modelo toman los siguientes valores:
rmse = rmse(test_df_multiple_media$bicis_disponibles, test_df_multiple_media$predicho)
rmse## [1] 5.072797mae = mae(test_df_multiple_media$bicis_disponibles, test_df_multiple_media$predicho)
mae## [1] 3.982729r2 = cor(test_df_multiple_media$bicis_disponibles, test_df_multiple_media$predicho)^2
r2## [1] 0.3382811Concluimos, por tanto, que el modelo es aceptable en cuanto a las predicciones que realiza
Para analizar un caso concreto, a continuación, se muestra un gráfico que muestra la evolución del número de bicicletas disponibles (reales en azul y predicha en naranja) el día 2023-02-14 en la estación 267.
est_info <- test_df_multiple_1_semanas[test_df_multiple_1_semanas$id_estacion == '267' & test_df_multiple_1_semanas$fecha == '2023-02-14', ]
trace_real <- est_info$bicis_disponibles
trace_predicho <- est_info$predicho
data <- data.frame(hora=c(0:23), trace_real, trace_predicho)
fig <- plot_ly(data, x = ~hora, y = ~trace_real, name = 'Real', type = 'scatter', mode = 'lines+markers')
fig <- fig %>% add_trace(y = ~trace_predicho, name = 'Predicho', mode = 'lines+markers')
fig <- fig %>%
layout(hovermode = "x unified", title = 'Evolución de las bicis disponibles en la estación 267 el 2023-02-14', plot_bgcolor = "#e5ecf6",
xaxis = list(title = 'Hora'), yaxis = list(title = 'Bicis disponibles'))
figAdicionalmente, se ha creado una página web en la cual los usuarios pueden consultar cuantas bicis habrá disponibles en una estación, un día de la semana, a una hora concreta. Debido a que para poder utilizar el modelo de regresión creado, sería necesario disponer de datos en tiempo real, el modelo que subyace a la web es el mencionado modelo basado en medias. Por el momento, la web puede ser consultada en este enlace y dejamos como propuesta de futuro la consideración de datos en tiempo real, lo que permitirá utilizar el modelo de regresión y que, por tanto, proporcionará predicciones con un error menor.
El código que crea esta aplicación es el siguiente:
# Cargamos las librerias
library(shiny)
library(readxl)
library(leaflet)
library(dplyr)
library(stringr)
# Cargamos los datos
load('df_merge.RDa')
load('df_bicis_total.RDa')
estaciones = read_excel("valenbisi.xlsx")
estaciones <- estaciones %>%
mutate(longitud = as.numeric(longitud),
latitud = as.numeric(latitud))
fecha_inicial_train <- as.Date("2022-12-02")
fecha_final_test <- as.Date("2023-02-09")
train_df <- df_merge[df_merge$fecha >= fecha_inicial_train & df_merge$fecha < fecha_final_test, ]
fecha_inicial_test <- as.Date("2023-02-09")
fecha_final_test <- as.Date("2023-02-23")
test_df <- df_merge[df_merge$fecha >= fecha_inicial_test & df_merge$fecha <= fecha_final_test, ]
# Cargamos la función get_matrix() y predict_bikes_test()
source("get_matrix.R")
source("predict_bikes_mean.R")
# Creamos la interfaz de usuario
ui <- fluidPage(
titlePanel("Mapa de las estaciones de Valenbisi"),
leafletOutput("map"),
titlePanel("Estimador de bicicletas disponibles"),
sidebarLayout(
sidebarPanel(
selectInput("estacion", "Selecciona la estación:",
choices = unique(df_merge$id_estacion),
selected = 1),
selectInput("dia_semana", "Selecciona el dia de la semana:",
choices = c('Lunes', 'Martes', 'Miercoles','Jueves','Viernes','Sábado','Domingo'),
selected = 1),
sliderInput("hora", "Selecciona la hora:",
min = 0, max = 23, value = 12, step = 1)
),
mainPanel(
mainPanel(style = "width: 100%;",
wellPanel(
h4("Número estimado de bicicletas disponibles:"),
fluidRow(
column(6,
wellPanel(
tags$h5("Predicción", class = "text-center"),
tags$div(verbatimTextOutput("prediccion"), class = "text-center")
)
),
column(6,
wellPanel(
tags$h5("Desviación típica", class = "text-center"),
tags$div(verbatimTextOutput("desv"), class = "text-center")
)
)
),
tags$h5(verbatimTextOutput("texto"), class = "text-center")
)
)
)
)
)
# Creamos el servidor
server <- function(input, output, session) {
values <- reactiveValues(estacion = "", popup_text = "")
output$map <- renderLeaflet({
leaflet() %>%
addTiles() %>%
addMarkers(data = estaciones, lng = ~longitud, lat = ~latitud, popup = ~paste("<b>", name, "</b></br>", actionButton("selectlocation", "Selecciona esta estación",
onclick = 'Shiny.onInputChange(\"button_click\", this.parentNode.parentNode.firstChild.innerHTML); Shiny.onInputChange(\"estacion\", estaciones$number_[which(estaciones$name==name)])')),
icon = makeIcon(iconUrl = "icono_bici.png", iconWidth = 20, iconHeight = 20))
})
# Función para actualizar la estación
update_estacion <- function() {
values$popup_text <- input$popup_text
values$estacion <- input$estacion
}
# Observador para obtener el texto del popup
observeEvent(input$button_click, {
values$popup_text <- input$button_click
estacion <- input$estacion
estacion <- str_extract(values$popup_text, "(?<=<b>\\s)\\d+(?=_)")
estacion <- gsub("^0+", "", estacion)
updateSelectInput(session, "estacion", selected = estacion)
update_estacion()
})
# Definimos la salida de la predicción
output$prediccion <- renderText({
# Obtenemos la predicción
dias = c('Lunes', 'Martes', 'Miercoles','Jueves','Viernes','Sábado','Domingo')
abrev = c('L','M','X','J','V','S','D')
for (i in 1:length(dias)){
if (input$dia_semana == dias[i]){
dia = abrev[i]
break
}
}
res <- predict_bikes_mean(input$estacion, dia, input$hora)
prediccion = res$media
return(round(prediccion, 2))
})
output$desv <- renderText({
# Obtenemos la predicción
dias = c('Lunes', 'Martes', 'Miercoles','Jueves','Viernes','Sábado','Domingo')
abrev = c('L','M','X','J','V','S','D')
for (i in 1:length(dias)){
if (input$dia_semana == dias[i]){
dia = abrev[i]
break
}
}
res <- predict_bikes_mean(input$estacion, dia, input$hora)
error = res$desviación
return(round(error, 2))
})
output$texto <- renderText({
dias = c('Lunes', 'Martes', 'Miercoles','Jueves','Viernes','Sábado','Domingo')
abrev = c('L','M','X','J','V','S','D')
for (i in 1:length(dias)){
if (input$dia_semana == dias[i]){
dia = abrev[i]
break
}
}
res <- predict_bikes_mean(input$estacion, dia, input$hora)
prediccion = res$media
error = res$desviación
paste0("En la estación ", input$estacion, " a las ", input$hora,":00", " de un ", input$dia_semana, " se estima que hay ", round(prediccion, 2), " +- ", round(error, 2), " bicis disponibles.")
})
}
# Lanzamos la aplicación
shinyApp(ui = ui, server = server)