TFM - Jorge Fernández Becerra
2023-09-19
Introducción
Para la realización de este trabajo contamos con un dataset que recoge la demanda de bicicletas por horas en Seúl a lo largo de un año. En el dataset se incluyen, además, datos meteorológicos e información relativa a la fecha. Con estos datos, se va a intentar solucionar un problema de predicción de la demanda, es decir, un problema de regresión. Se ha considerado este tipo de problema ya que saber cuántas bicicletas se van a demandar puede ayudar a la hora de tomar decisiones como retirar bicicletas temporalmente para su mantenimiento o tener tan solo las bicicletas necesarias para cubrir la demanda estando almacenadas el resto evitando así que sufran las inclemencias del clima y alargando su vida útil.
Decidido, el tipo de aplicación queda decidir con cuánto tiempo de antelación queremos predecir la demanda. Una predicción a demasiado corto plazo puede no ser útil, ya que no proporciona el tiempo suficiente para adaptarse. Por otra parte, intentar predecir a demasiado largo plazo dificultaría las predicciones.Teniendo todo esto en cuenta se ha decidido crear un modelo que prediga la demanda con 24h de antelación.
Queda también decidir la granularidad de la predicción. En el dataset los datos están agregados por horas y aunque podríamos agruparlos por días o semanas para predecir las demandas diarias o semanales considero que es más útil para el servicio de alquiler conocer la predicción de la demanda por horas, ya que le permite saber cuando empezar a retirar y reponer las bicicletas.
En conclusión, el objetivo de este trabajo será crear un modelo que prediga la demanda de bicicletas por horas en el servicio de alquiler de bicicletas de Seúl con una antelación de 24 horas.
Otra cuestión, es como utilizar los datos meteorológicos.En el dataset tenemos los datos meteorológicos para cada hora. Sin embargo, cuándo se utilice el modelo no tendrá los datos meteorológicos del día siguiente en todo caso tendrá las predicciones. Si bien es cierto que a día de hoy los modelos meteorológicos, especialmente en el corto plazo, son altamente precisos sigue existiendo cierto grado de incertidumbre. Está incertidumbre puede afectar negativamente al modelo cuando se implemente en escenarios reales, especialmente si el modelo se vuelve altamente sensible a los datos meteorológicos. Dicho de otra manera, estaríamos convirtiendo el problema en una cascada de decisiones donde un error en la primera etapa (predicción del clima) se propaga a la segunda etapa (predicción de la demanda). Además, con los datos disponibles no podemos valorar como afecta esta incertidumbre al modelo (para ello necesitaríamos un histórico de las predicciones). Por todo ello, he decidido no incluir datos meteorológicos del día en cuestión, aunque si que se pueden introducir datos meteorológicos de días pasados, ya que al ser observaciones reales y no predicciones puede proporcionar información valiosa sin introducir incertidumbre adicional.
Transformación Inicial
En primer lugar vamos a transformar los datos añadiendo una columna con el día de la semana y otra con el mes. A continuación crearemos un dataset que agrupe la demanda por días, para poder analizar la influencia de los distintos eventos relacionados con las fechas. De esta forma evitamos las distorsiones que puedan introducir en el análisis las fluctuaciones de la demanda a los largo del días. En este dataset omitiremos los días en los que el servicio no está en funcionamiento.
#Vamos a cargar el dataset con el que trabajaremos
SeoulBikeData <- read.csv("Datos Originales.csv", header=TRUE, fileEncoding="latin1")
#cambiamos los nombres de las columnas
colnames(SeoulBikeData) = c("Fecha", "Num_Alquileres", "Hora", "Temperatura", "Humedad", "Vel_viento",
"Visibilidad", "Punto_rocio", "Rad_solar", "Cant_lluvia", "Cant_nieve",
"Estacion", "Festivo", "Dia_de_func")
#Aplicamos los cambios que hemos definido en el análisis exploratorio
SeoulBikeData <- SeoulBikeData %>%
mutate(Fecha = as.Date(Fecha, format="%d/%m/%Y")) %>% #Convertimos la columna "Fecha" en tipo "Date"
select(-Punto_rocio) #Eliminamos la variable de "Punto_rocio" ya que estaba altamente correlacionada con la temperatura
# Creamos los vectores con el orden de los días y los meses
dias = c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo")
meses = c("enero", "febrero", "marzo", "abril", "mayo", "junio", "julio", "agosto",
"septiembre", "octubre", "noviembre", "diciembre")
# Añadimos columna de días y de meses
SeoulBikeData <- SeoulBikeData %>%
mutate(
Dia_Semana = factor(weekdays(Fecha), levels = dias),
Mes = factor(months(Fecha), levels = meses)
)
# Creamos un dataset con los datos agrupados por días
Demanda_Diaria <- SeoulBikeData %>%
group_by(Fecha, Dia_de_func, Festivo, Dia_Semana, Mes) %>%
summarise(Num_Alquileres = sum(Num_Alquileres, na.rm = TRUE), .groups = "drop") %>%
mutate(Dia_Semana = factor(Dia_Semana, levels = dias)) %>%
filter(Dia_de_func == "Yes")Análisis Descriptivo y transformación de datos
Se ha decido combinar la sección el análisis descriptivo y la transformación de los datos. Esto se debe a que, a medida que se transformen los datos evaluaremos su relevancia mediante técnicas descriptivas e inferenciales. Por lo tanto, se ha considerado que lo más natural y coherente es juntar las dos partes.
Vamos a segmentar el análisis descriptivo según la naturaleza del dato en las siguientes secciones:
- Análisis relacionado con fechas y temporalidades
- Análisis de eventos climáticos.
- Análisis de la demanda
Analisis relación de la demanda con fechas y temporalidades
En este punto analizaremos los atributos relacionados con el tiempo cronológico, en concreto, analizaremos las siguientes variables:
- Hora
- Día de la semana
- Día festivo
- Meses del año
Dado que estamos intentando predecir la demanda horaria, vamos a comenzar analizando este atributo y, a partir de este, expandiremos el análisis al resto de variables. Para evitar que valores atípicos distorsionen los datos vamos a representar la mediana de unidades demandadas por hora, ya que la mediana es una medida robusta a valores atípicos. Para poder tener información de la distribución de la demanda por horas vamos a graficar también el rango intercuartílico.
#Creamos dataset para obtener los estadísticos relacionados con las horas
Resumen_horas <- SeoulBikeData %>%
filter(Dia_de_func == "Yes") %>%
group_by(Hora) %>%
summarise( mediana = median(Num_Alquileres),
p25 = quantile(Num_Alquileres, 0.25),
p75 = quantile(Num_Alquileres, 0.75))
#Creamos el gráfico
ggplot(data = Resumen_horas , aes(x = Hora, y = mediana)) +
geom_line(aes(y = mediana, color = "Mediana"), linewidth = 1) +
geom_ribbon(aes(ymin = p25, ymax = p75, fill = "Rango Intercuartílico"), alpha = 0.2) +
ggtitle("Mediana de unidades demandadas por hora") +
labs(y = "Uds. de bicicletas demandadas", x = "", color = "Leyenda") +
scale_color_manual(values = c("Mediana" = "blue"), name="") +
scale_fill_manual(values = c("Rango Intercuartílico" = "blue"), name = "" ) +
theme_minimal() +
theme(legend.position = "bottom")
#Borramos el dataset para ahorrar espacio en memoria
rm(Resumen_horas)Lo primero que destaca es que la demanda, tiene dos picos a lo largo del día, probablemente debido a la entrada y salida del trabajo. También vemos que, tras el primer pico de la demanda, el rango intercuartílico es bastante amplio, lo que indica bastante disparidad de la demanda a la misma hora en distintos días, probablemente causadas por el resto de variables que todavía no hemos estudiado.
Vamos a añadir una dimensión más al gráfico anterior, los días de la semana, de tal forma que analizaremos como es el patrón horario de la demanda en los distintos días de la semana.
#Creamos dataset para obtener los estadísticos relacionados con las horas
# Creamos dataset para obtener los estadísticos relacionados con las horas y días de la semana
Resumen_horas_dias <- SeoulBikeData %>%
filter(Dia_de_func == "Yes") %>%
group_by(Hora, Dia_Semana) %>%
summarise(mediana = median(Num_Alquileres), .groups = "drop")
# Creamos el gráfico
ggplot(data = Resumen_horas_dias , aes(x = Hora, y = mediana)) +
geom_line(aes(color = Dia_Semana), linewidth = 1) +
ggtitle("Ev. de la demanda de bicicletas por horas según el día de la semana") +
labs(y = "Mediana de bicicletas demandadas", x = "", color = "Día de la Semana") +
theme_minimal() +
theme(legend.position = "bottom")
#Borramos el dataset para ahorrar espacio en memoria
rm(Resumen_horas_dias)Observamos una distinción significativa en los patrones horarios entre los fines de semana y los días laborables. Aunque dentro de los días laborables las variaciones en la demanda son mínimas, es notable que, si bien la evolución de la demanda a lo largo del día es parecida para sábados y domingos, la mediana de la demanda los domingos es ligeramente menor.
Veamos ahora como es el patrón de demanda horario para los días festivos comparándolo con el que, a priori, debería ser el día más similar en comportamiento, los domingos. Para ello representaremos de nuevo la demanda y los rangos intercuartílicos. Para el análisis gráfico tendremos solo en cuenta los datos a partir de 2018, ya que, como vimos en el análisis exploratorio, de los 17 festivos 4 tenían lugar en EN diciembre, lo que podría sesgar los datos.
# Creamos dataset para obtener los estadísticos relacionados con las horas y días de la semana
Resumen_dom_fest <- SeoulBikeData %>%
mutate(Dia_Semana = as.character(Dia_Semana)) %>%
mutate(Dia_Semana = ifelse(Festivo == "Holiday", "Festivo", Dia_Semana)) %>%
filter(Dia_de_func == "Yes", Dia_Semana %in% c("domingo", "Festivo"), Mes != "diciembre") %>%
group_by(Hora, Dia_Semana) %>%
summarise(
mediana = median(Num_Alquileres),
p25 = quantile(Num_Alquileres, 0.25),
p75 = quantile(Num_Alquileres, 0.75),
.groups = "drop",)
# Creamos el gráfico
ggplot(data = Resumen_dom_fest , aes(x = Hora, y = mediana)) +
geom_ribbon(aes(ymin = p25, ymax = p75, fill = Dia_Semana), alpha = 0.3) + # Agregamos el ribbon transparente
geom_line(aes(color = Dia_Semana), linewidth = 1) +
ggtitle("Patrón demanda horaria domingos y festivos") +
labs(y = "Uds. de bicicletas demandadas", x = "", color = "Día de la Semana", fill = "Día de la Semana") +
theme_minimal() +
theme(legend.position = "bottom")
#Borramos el dataset para ahorrar espacio en memoria
rm(Resumen_dom_fest)Vemos que el patrón horario es muy similar, cabria preguntarse si podemos unir las dos características creando un atributo “domingo/festivo” para los días que sean domingo o festivo, de esta forma solucionaríamos el problema del poco muestreo de los festivos. Para comprobar si las distribuciones de los datos de ambos días son iguales vamos aplicar el test de Kolmogorov-Smirnov, cuya hipótesis nula es que las dos muestras provienen de la misma distribución. La elección de este test se debe a que es un test no paramétrico y no necesita que la distribución de datos sea normal ya que, como hemos visto anteriormente, la demanda horaria no sigue una distribución normal. Además, este test resulta adecuado para muestras pequeñas. A pesar de que en el análisis gráfico hemos excluido los festivos de diciembre, en este contraste de hipótesis los vamos a incluir, para poder aumentar el muestreo y que sea lo más representativo posible. Vamos a aplicar un test para cada hora y representar el p-valor obtenido:
#Creamos dataframe vacio donde se anotarán los resultados
resultados_ks <- data.frame(Hora=integer(), KS_statistic=double(), p_value=double())
# Iteramos sobre cada hora y realizamos el ks.test
for (hora in 0:23) {
muestra_festivos <- SeoulBikeData %>% filter(Festivo == "Holiday", Hora == hora)
muestra_domingos <- SeoulBikeData %>% filter(Dia_Semana == "domingo" & Festivo != "Holiday" , Hora == hora)
test <- ks.test(muestra_festivos$Num_Alquileres, muestra_domingos$Num_Alquileres)
resultados_ks <- rbind(resultados_ks, data.frame(Hora=hora, KS_statistic=test$statistic, p_value=test$p.value))
}
# Creamos un grafico para visualizar los resultados
ggplot(resultados_ks, aes(x=Hora, y=p_value)) +
geom_line(color="#1D8077", linewidth=1) +
geom_point(color="#18A396") +
geom_hline(yintercept = 0.05, linetype = "dashed", color = "#0C720F") +
annotate("text", x = Inf, y = 0.055, label = "Acepto Ho: las distribuciones no son iguales", hjust = 1.95, vjust = 0, color = "#49AB33") +
labs(title="Resultados p-valor por hora",
x="Hora",
y="p-valor") +
scale_y_continuous(limits = c(0, NA)) +
theme_minimal()
#Borramos los datasets utilizados
rm(resultados_ks, muestra_domingos, muestra_festivos, test)El p-valor es mayor que 0.05 en todas las horas, lo que sugiere que no podemos rechazar la Hipótesis Nula de igualdad de distribuciones. Me gustaría señalar que, al ser la muestra de días festivos reducida, corremos un mayor riesgo de cometer un error de Tipo II, es decir, no rechazar la Hipótesis nula cuando en realidad no es falsa. Aun así, se ha decidido transformar la variable día de la semana cambiando todos los días que sean domingo o festivo o como domingo/festivo, ya que elimina el problema del bajo muestreo de festivos y reduce la dimensionalidad del modelo.
SeoulBikeData <- SeoulBikeData %>%
mutate(Festivo = ifelse(Dia_Semana == "domingo", "Holiday", Festivo)) %>%
mutate(Tipo_Dia = ifelse(
Dia_de_func == "No", "No Funcional", ifelse(
Festivo == "Holiday", "domingo/festivo", ifelse(
Dia_Semana == "sábado", "sábado", "laborable"
))))Veamos ahora si existen diferencias en la distribución horaria de la demanda dependiendo del mes del año. Para ello, representaremos la mediana de la demanda horaria por cada mes y hora. Dado que queremos analizar si las distribución de la demanda es similar, normalizaremos los datos al máximo valor de la mediana para cada mes:
# Calcula la mediana de Num_Alquileres por Mes y Hora
resumen_mediana <- SeoulBikeData %>%
group_by(Mes, Hora) %>%
summarise(Mediana_Alquileres = median(Num_Alquileres, na.rm = TRUE), .groups = "drop")
# Calcula el valor máximo de la mediana para cada Mes
resumen_max <- resumen_mediana %>%
group_by(Mes) %>%
summarise(Max_Mediana = max(Mediana_Alquileres, na.rm = TRUE), .groups = "drop")
# Une los datos de la mediana con los valores máximos por Mes
mediana_normalizada <- merge(resumen_mediana, resumen_max, by = "Mes")
# Calcula la mediana normalizada
mediana_normalizada$Mediana_Normalizada <- mediana_normalizada$Mediana_Alquileres / mediana_normalizada$Max_Mediana
# Grafica los datos normalizados
ggplot(mediana_normalizada, aes(x = Hora, y = as.factor(Mes), fill = Mediana_Normalizada)) +
geom_tile() +
scale_fill_gradient2(low = "#002BFF", mid= "white", high = "#23E700", midpoint = 0.5) +
labs(title = "Dist. demanda por mes y hora",
x = "Hora",
y = "Mes",
fill = "Mediana Normalizada") +
theme_minimal() +
theme(legend.position = "bottom")
#Borramos los dataframes creados para ahorrar especio en memoria
rm(resumen_mediana, resumen_max, mediana_normalizada)Las horas pico tienden a respetarse en todos los meses, especialmente la 2ª hora pico, que es la mayor. Además, en los meses más cálidos (de junio a octubre) la demanda se mantiene relativamente estable tras la 2ª hora pico. El resto de meses, se produce el efecto contrario, hay una mayor demanda antes de la hora pico y durante los meses más cálidos esas horas tienen menos demanda. Probablemente esto se deba a que, con el calor, la gente prefiera no usar el servicio en las horas más cálidas del día. Dado la limitación de datos del conjunto de datos, un año de históricos, se ha decidido no incluir variables “mes” ni “Estacion”, ya que no tenemos forma de evaluar si las rutinas mensuales o estacionales de un año especifico se repiten el resto de años.
Analisis descriptivo de la demanda
En primer lugar, vamos cual ha sido la evolución de la demanda diaria de bicicletas a lo largo de todo el histórico.
#Creamos dataframe con la media móvil
Demanda_Med_Mov <- Demanda_Diaria %>%
arrange(Fecha) %>%
mutate(
media_movil_30d = zoo::rollmean(Num_Alquileres, 30, fill = NA, align = "right")
)
#Creamos gráfico
ggplot(data = Demanda_Med_Mov , aes(x = Fecha, y = Num_Alquileres)) +
geom_line(color = "grey") +
geom_line(aes(y = media_movil_30d, color = "30 días"), linewidth = 1, na.rm = TRUE) +
ggtitle("Ev. demanda de bicicletas") +
labs(color = "Media Móvil") +
theme_minimal() +
theme(legend.position = "bottom")
#Borramos dataframes
rm(Demanda_Med_Mov)Lo primero que podemos destacar de la evolución histórica es la inestabilidad de la demanda ya que, aunque hay una tendencia general, se aprecian súbitos bajones en la demanda. Mas adelante estudiaremos que factores influyen en estos descensos. Por otra parte, también se aprecia que en los últimos días del histórico la demanda es mucho mayor que al principio.
Una importante fuente de información de la que podría disponer el modelo de ML es la demanda que ha habido de bicicletas un determinado intervalo de horas previas. Recordemos que el objetivo del modelo es predecir la demanda que va a haber con 24h de antelación, por lo que el nº de horas previos mínimo que deberemos considerar serán 24h. Para el cálculo de estás horas previas excluiremos los días en los que el servicio no estuvo en funcionamiento, es decir, si analizamos la demanda registrada 24h antes nos referimos a 24h horas efectivas de servicio, que no siempre coincidirá con un período real de 24h. Para establecer que período temporal utilizaremos un gráfico que represente la correlación entre la demanda y la demanda que hubo a X horas previas. Destacar que este enfoque tiene una importante debilidad, la correlación tan solo captura la relación líneal pero podría ser posible que hubiera otro tipo de correlación que estuvieramos obviando.
# Creamos un dataframe vacío para guardar las correlaciones
Estudio_Corr <- data.frame()
EstudioDemanda <- SeoulBikeData %>% filter( Dia_de_func == "Yes")
# Bucle para calcular la correlación a las distintas horas
for(i in 24:96) {
EstudioDemanda$Demanda_Mov <- lag(EstudioDemanda$Num_Alquileres, i)
# Calculamos la correlación
correlacion <- cor(EstudioDemanda$Num_Alquileres, EstudioDemanda$Demanda_Mov, use="complete.obs")
# Insertamos la correlación
Estudio_Corr <- rbind(Estudio_Corr, data.frame(Hora=i, Correlacion=correlacion))
}
# Obtemos las coordenadas de 24,48,72 y 96 horas antes
coordenadas <- Estudio_Corr %>% filter( Hora %in% c( 24, 48, 72, 96))
# Creamos el gráfico para ver como evolucionan las correlaciones
ggplot(Estudio_Corr, aes(x=Hora, y=Correlacion)) +
geom_line() +
geom_point(data=coordenadas, aes(x=Hora, y=Correlacion), color="blue", size=3) +
geom_text(data=coordenadas, aes(x=Hora, y=Correlacion, label=sprintf("%.0f", Hora)), vjust=-1) +
labs(title="Correlación de la demanda móvil",
x="Horas",
y="Coef. de correlación") +
scale_y_continuous(limits=c(min(Estudio_Corr$Correlacion, na.rm=TRUE), 0.8)) +
theme_minimal()
#Borramos los dataframes creados
rm(Estudio_Corr, EstudioDemanda, coordenadas)El análisis revela que la máxima autocorrelación de la demanda se observa con un desfase de 24 horas, presentando un coeficiente cercano a 0.8. Adicionalmente, la autocorrelación exhibe un patrón cíclico con un periodo de 24 horas, lo cual es coherente con la naturaleza de los datos ya que tiene mas sentido la demanda de bicicletas correlacione mas fuerte con datos registrados a la misma hora.
Antes de decidir de incluir esta variable, vamos a tener otro factor en consideración. Como veíamos previamente, hay diferencias notables entre el patrón de demanda horario de los días laborables, de los sábados y de los domingos/festivos. Debido a esto, vamos a crear una nueva característica que recoja la demanda que había 24 horas funcionales antes pero para cada tipo de día se van a tener en cuenta solo las horas funcionales de cada tipo de día. Por ejemplo, para la demanda del Lunes a las 15h se tendrá en cuenta la demanda del viernes a las 15h, no del domingo. Vamos a comparar cual es el grado de correlación con la demanda para cada hora de esta nueva característica, a la que llamaremos “Demanda 24h ajustada”.
# Creamos un datas con las demandas 24h previas
Estudio_corr_horas <- SeoulBikeData %>%
filter(Dia_de_func == "Yes") %>%
mutate(Demanda_24h = lag(Num_Alquileres, 24)) %>%
group_by(Tipo_Dia) %>%
mutate(Demanda_24h_aj = ifelse(Tipo_Dia == lag(Tipo_Dia, 24), lag(Num_Alquileres, 24), NA)) %>%
ungroup()
# Creamos un dataset con las correlaciones para cada tipo de demanda
correlaciones <- Estudio_corr_horas %>%
group_by(Hora) %>%
summarize(
correlacion = cor(Demanda_24h, Num_Alquileres, use = "complete.obs"),
correlacion_aj = cor(Demanda_24h_aj, Num_Alquileres, use = "complete.obs")
)
# Guardamos el valor promedio de las correlaciones para presentarlas en el gráfico
avg_24h <- mean(correlaciones$correlacion, na.rm = TRUE)
avg_24h_aj <- mean(correlaciones$correlacion_aj, na.rm = TRUE)
# Creamos el gráfico
ggplot(correlaciones, aes(x = Hora)) +
geom_line(aes(y = correlacion, colour = "Correlación sin ajuste"), linewidth = 0.75) +
geom_line(aes(y = correlacion_aj, colour = "Correlación con ajuste"), linewidth = 0.75) +
geom_hline(aes(yintercept = avg_24h, colour = "Correlación sin ajuste"), linetype = "dashed") +
geom_hline(aes(yintercept = avg_24h_aj, colour = "Correlación con ajuste"), linetype = "dashed") +
annotate("text", x = Inf, y = avg_24h - 0.01, label = sprintf("Prom.: %.2f", avg_24h), hjust = 1, color = "#0058D7") +
annotate("text", x = Inf, y = avg_24h_aj - 0.01, label = sprintf("Prom.: %.2f", avg_24h_aj), hjust = 1, color = "#F09800") +
labs(title = "Corr.entre Demanda 24h con y sin ajuste y la demanda",
x = "Hora",
y = "Coeficiente de correlación") +
scale_colour_manual(values = c("#0058D7", "#F09800"),
name = "",
breaks = c("Correlación sin ajuste", "Correlación con ajuste")) +
theme_minimal() +
theme(legend.position = "bottom")
#Borramos dataframes y variables utilizados
#rm(Estudio_corr_horas, correlaciones, avg_24h, avg_24h_aj)Vemos que el promedio de correlación de la “Demanda 24h ajustada” es ligeramente superior a la “Demanda 24 sin ajustar”, sin embargo, esto no es el factor mas relevante a la hora de elegir una modalidad u otra, sino el hecho de que en las horas correspondientes al primer pico de la demanda tenga una mayor correlación el la “Demanda 24 ajustada”. Si incorporáramos la “Demanda 24 sin ajustar” estaríamos dificultando al modelo una buena predicción precisamente en uno de los puntos más críticos para la logística.
# Creamos la variable que recoja la demanda hace 24h
SeoulBikeData <- SeoulBikeData %>%
group_by(Tipo_Dia) %>%
mutate(Demanda_24h_aj = ifelse(Tipo_Dia == lag(Tipo_Dia, 24), lag(Num_Alquileres, 24), NA)) %>%
ungroup()Otra posible fuente de información para la construcción del modelo, mas allá de la demanda en una hora en concreto, es la demanda acumulada en un determinado período de horas. Hay que mencionar, que este período no puede tener en cuenta las 24h previas a la demanda, ya que no tendrá esa información a la hora de realizar la predicción. Vamos a realizar una análisis similar al que realizamos con la demanda a 24h y analicemos como es la correlación en cada ventana de horas:
# Creamos un dataframe vacío para guardar las correlaciones
Estudio_Corr <- data.frame()
EstudioDemanda <- SeoulBikeData %>%
filter(Dia_de_func == "Yes") %>%
group_by(Tipo_Dia) %>%
mutate( demanda_acum_24 = rollapply(Num_Alquileres, width = 24, FUN = sum, align = "right", fill = NA) ) %>%
ungroup()
# Bucle para calcular la correlación a las distintas horas
for(i in 36:120) {
EstudioDemanda <- EstudioDemanda %>%
group_by(Tipo_Dia) %>%
mutate( Demanda_Acumulada = rollapply(Num_Alquileres, width = i, FUN = sum, align = "right", fill = NA) - demanda_acum_24 ) %>%
ungroup()
# Calculamos la correlación
correlacion <- cor(EstudioDemanda$Num_Alquileres, EstudioDemanda$Demanda_Acumulada, use="complete.obs")
# Insertamos la correlación
Estudio_Corr <- rbind(Estudio_Corr, data.frame(Hora=i, Correlacion=correlacion))
}
# Obtemos las coordenadas del valor mas alto
coordenadas <- Estudio_Corr %>% top_n(1, Correlacion)
# Creamos el gráfico para ver como evolucionan las correlaciones
ggplot(Estudio_Corr, aes(x=Hora, y=Correlacion)) +
geom_line() +
geom_point(data=coordenadas, aes(x=Hora, y=Correlacion), color="blue", size=3) +
geom_text(data=coordenadas, aes(x=Hora, y=Correlacion, label=sprintf("%.0f", Hora)), vjust=0, hjust= -0.7) +
labs(title="Correlación de la demanda móvil acumulada",
x="Horas",
y="Correlación") +
theme_minimal()
#Borramos los dataframes creados
rm(Estudio_Corr, correlacion)En este caso la mayor correlación se produce al acumular la demanda entre 24h y 53h antes, por lo que incorporaremos esta característica al conjunto de datos para entrenar al modelo.
#Creamos demanda acumulada
SeoulBikeData <- SeoulBikeData %>%
group_by(Tipo_Dia) %>%
mutate(Demanda_Acum_24_53h =
rollapply(Num_Alquileres, width = 53, FUN = sum, align = "right", fill = NA) -
rollapply(Num_Alquileres, width = 24, FUN = sum, align = "right", fill = NA) ) %>%
ungroup()
#Creamos un array con los nombres de las columnas relacionadas con la demanda
Nombres_col_demanda = c("Demanda_24h_aj", "Demanda_Acum_24_53h")Para terminar este punto, vamos a describir gráficamente la correlación entre las dos variables que hemos creado y la demanda:
#Gráfico distribución demanda 24h
p1 <- ggplot(data = SeoulBikeData, aes(x=Num_Alquileres, y = Demanda_24h_aj)) +
scale_fill_distiller(palette=4, direction=-1) +
geom_point(size=1, alpha=0.1, color="#046ECB") +
stat_density_2d(aes(fill = ..level..), geom = "polygon", colour="white") +
ggtitle("Corr. Demanda y Demanda 24h aj.") +
theme_minimal() +
theme(legend.position = "none")
#Gráfico distribución demanda acum 24-53h
p2 <- ggplot(data = SeoulBikeData, aes(x=Num_Alquileres, y = Demanda_Acum_24_53h)) +
scale_fill_distiller(palette=4, direction=-1) +
geom_point(size=1, alpha=0.1, color="#046ECB") +
stat_density_2d(aes(fill = ..level..), geom = "polygon", colour="white") +
ggtitle("Corr. Demanda y Demanda Acum. 24-53h") +
theme_minimal() +
theme(legend.position = "none")
p1 + p2
Analisis descriptivo de la relación del clima y la demanda
En primer lugar, vamos a representar gráficamente la relación de la demanda con las distintas variables climáticas. Para ello vamos a trabajar con un dataset en el que solo estén los datos pertenecientes a las 18h de los días laborables. Se ha decidido filtrar los datos para aislar de la demanda los efectos del día y la hora, que como hemos visto, tienen una gran influencia. Por otra parte, se ha elegido las 18h por ser el pico de mayor demanda y los días laborables porque representan la mayoría de datos disponibles. La principal desventaja de introducir este filtro es que podemos pasar por alto información de como afecta el clima a distintos tipos de días. Aun así, teniendo en cuenta que la idea general de esta parte del análisis es poder tener una visión general del efecto del clima en la demanda, se puede aceptar este inconveniente.
#Recogemos las variables climáticas y las estaciones en vectores para facilitar el trabajo
variables_clima <- c("Temperatura", "Humedad", "Vel_viento", "Visibilidad", "Rad_solar","Cant_lluvia", "Cant_nieve")
Estaciones = c("Winter", "Spring", "Summer", "Autumn")
#Asignamos los niveles del factor a la variable "Estacion"
SeoulBikeData$Estacion <- factor(SeoulBikeData$Estacion, levels = Estaciones)
#Definimos los colores para las estaciones, para tener el mismo criterio en todos los gráficos
colores_estaciones <- c("Winter" = "#84CEFF",
"Spring" = "#77dd77",
"Summer" = "#FC9F20",
"Autumn" = "#B86B1A")
#Creamos el dataset con los datos que usaremos en los graficos
SeoulBikeData_estudio_clima <- SeoulBikeData %>%
filter(Dia_de_func == "Yes", Tipo_Dia == "laborable", Hora == 18) %>%
select(c(Num_Alquileres, variables_clima, Estacion))
#Creamos los gráficos
SeoulBikeData_estudio_clima %>%
gather(-Num_Alquileres, -Estacion, key = "var", value = "value") %>%
ggplot(aes(x = value, y = Num_Alquileres)) +
facet_wrap(~ var, scales = "free") +
geom_point(aes(color = Estacion) ,alpha = 0.6, size = 0.75) +
stat_smooth(method = "loess", formula = "y ~ x") +
scale_color_manual(values = colores_estaciones) +
theme_minimal() +
theme(legend.position = "bottom", axis.title.x=element_blank())
Vamos a analizar los de la correlación entre las variables climáticas y la demanda de bicicletas:
Cantidad de lluvia: Los datos parecen indicar una relación binaria en los que, en caso de que llueva, la demanda es menos que en caso de que no llueva. Esto tiene sentido económico, ya que ante la lluvia lo esperable es que los posibles usuarios del servicio de bicis busquen un transporte alternativo. Aún así, haremos un examen más detallado para conocer con mayor precisión la magnitud de esta relación.
Cantidad de nieve: Se observa algo similar que en “Cantidad de lluvia”, ya que los datos de las horas que nieva presentan una distribución mas baja que en los datos en los que no nieva. En este caso, todos los días que nieva son en invierno, que tiene una demanda más baja que el resto del año, por lo que no queda del todo claro cuál es la incidencia de la haber nevado en la demanda. Por lo tanto, realizaremos un análisis teniendo en cuenta solo los datos de invierno.
Humedad: En primer lugar, hay que aclarar que la tendencia creciente que vemos entre el 0% y 50% se produce porque en invierno, que como vimos tenía una menor demanda, se concentra la mayor parte de datos con baja humedad lo que distorsiona la linea de tendencia. En el resto de estaciones del año se observa una relación negativa entre el 50% y el 100% de humedad. Las conclusiones no son claras, por lo que haremos un examen mas detallado separando los datos por estación del año para obtener unas conclusiones más claras.
Radiación solar: Hay una clara relación positiva, en todas las estaciones del año la demanda aumenta cuando hay una mayor radiación solar, lo cuál tiene sentido económico, ya que una mayor irradiación solar indica una mayor cantidad de luz solar lo que puede estimular la demanda de bicicletas.
Temperatura: También tiene una relación bastante clara, con una relación positiva hasta los 25º C momento en el cual la relación es negativa. Lo cuál tiene sentido dentro del negocio, ya que es razonable pensar que a mayor temperatura mas apetecible será el uso de la bicicleta pero, a partir de cierto punto, el calor será excesivo y empezará a desincentivar la demanda.
Velocidad del viento: No hay una relación clara entre la velocidad del viento y la demanda, por lo que examinaremos con detalle como es la relación en cada estación para poder obtener mas información acerca de la relación entre la velocidad del viento y la demanda.
Visibilidad: En todas las estaciones se observa como a medida que se reduce la visibilidad disminuye la demanda. Algo que tiene sentido económico ya que es lógico que una menor visibilidad desincentive el uso de la bicicleta, ya que las condiciones son mas inseguras.
Tal como comentamos anteriormente, vamos a examinar en mayor detalle como afecta la lluvia y la nieve a la demanda. En los gráficos de distribución observábamos que cuando llovía o nevaba los valores eran mas bajos, aunque al acumularse un gran número de puntos cuando el valor del eje X era 0, es decir cuando no llueve/nieva, nos perdíamos detalles de como eran las distribuciones. Por lo tanto, vamos a crear dos gráficos “boxplot” para poder ver las diferencias de las distribuciones entre los días que llueve y los que no llueve y entre los días que nieva y los que no nieva. Dado que, como hemos visto en análisis anterior, todos los días que nieva es invierno en el caso de la segunda comparativa utilizaremos solo datos de invierno.
SeoulBikeData_estudio_clima <- SeoulBikeData_estudio_clima %>%
mutate( lluvia = ifelse( Cant_lluvia > 0, "Llueve", "No llueve" ),
nieve = ifelse( Cant_nieve > 0, "Nieva", "No nieva")
)
p1 <- ggbetweenstats(
data = SeoulBikeData_estudio_clima,
x = lluvia,
y = Num_Alquileres,
results.subtitle = FALSE
)
SeoulBikeData_estudio_clima <- SeoulBikeData_estudio_clima %>%
mutate( lluvia = ifelse( Cant_lluvia > 0, "Llueve", "No llueve" ),
nieve = ifelse( Cant_nieve > 0, "Nieva", "No nieva")
)
SeoulBikeData_clima_winter <- SeoulBikeData_estudio_clima %>% filter( Estacion == "Winter")
p2 <- ggbetweenstats(
data = SeoulBikeData_clima_winter,
x = nieve,
y = Num_Alquileres,
pairwise.comparisons = TRUE,
pairwise.display = "all",
results.subtitle = FALSE,
bf.message = TRUE
)
print(p1 + p2)
Tanto la distribución de los datos como la demanda y el promedio parecen
dejarnos claros que la demanda es menor cuando llueve y cuando nieva a
los días en los que no lo hace. Para asegurarnos que esta diferencia es
significativa vamos a utilizar un contraste no parámetrico,
concretamente el test de suma de rangos de Wilkoxon-Mann-Whitney para
comprobar si la mediana de los días que llueva/nieva es menor que los
días que no lo hace.
#Realizamos el contraste de hipótesis para los días de lluvia
demanda_no_lluvia <- SeoulBikeData_estudio_clima$Num_Alquileres[SeoulBikeData_estudio_clima$lluvia == "No llueve"]
demanda_lluvia <- SeoulBikeData_estudio_clima$Num_Alquileres[SeoulBikeData_estudio_clima$lluvia == "Llueve"]
resultado_lluvia <- wilcox.test( demanda_lluvia, demanda_no_lluvia, alternative = "less")
#Realizamos el contraste de hipótesis para los días de nieve
demanda_no_nieve <- SeoulBikeData_estudio_clima$Num_Alquileres[SeoulBikeData_estudio_clima$nieve == "No nieva"]
demanda_nieve <- SeoulBikeData_estudio_clima$Num_Alquileres[SeoulBikeData_estudio_clima$nieve =="Nieva"]
resultado_nieve <- wilcox.test(demanda_nieve, demanda_no_nieve, alternative = "less")
#Imprimimos los p valores por pantalla
print( paste("Hip. Nula: mediana días con lluvia mayor o igual que días sin lluvia, p valor = ", round( resultado_lluvia$p.value, 10 )))## [1] "Hip. Nula: mediana días con lluvia mayor o igual que días sin lluvia, p valor = 6e-10"print( paste("Hip. Nula: mediana días con nieve mayor o igual que días sin nieve, p valor = ", round(resultado_nieve$p.value, 10 )))## [1] "Hip. Nula: mediana días con nieve mayor o igual que días sin nieve, p valor = 8.67e-08"#Eliminamos las variables que no vamos a necesitar
rm(demanda_no_lluvia, demanda_lluvia, resultado_lluvia, demanda_no_nieve,demanda_nieve, resultado_nieve)En ambos conjuntos de datos el p valor es prácticamente cero. Esto nos lleva a rechazar la hipótesis nula, la cual sugiere que durante los días de lluvia/nieve, la demanda de bicicletas es igual o mayor que en días sin dichas condiciones climáticas. En su lugar, nuestros datos sugieren una mayor probabilidad de que haya una menor demanda de bicicletas en días con lluvia o nieve. A pesar de que el muestreo es pequeño (17 y 12 datos) y normalmente deberíamos tomar con precaución los resultados, el hecho de que los resultados tengan un sentido lógico (que disminuya la demanda de bicis con lluvia y nieve) y que en la distribución de los datos que observamos en los gráficos se observe como la demanda es claramente menor nos permiten, en este caso, dar por válido el resultado.
Anteriormente, destacamos nuestra incertidumbre respecto a cómo la humedad y la velocidad del viento influían en la demanda de bicicletas. En lo que se refiere a la humedad, la mayor concentración de puntos de baja humedad en invierno, que es la época de menor demanda, distorsionaba los datos. Por otro lado, con la velocidad del viento, no pudimos identificar una relación evidente con la demanda. Sin embargo, creemos que al desglosar estos datos por estaciones, podríamos descubrir patrones que antes eran difíciles de ver. Para investigar esto más a fondo, utilizaremos gráficos de dispersión y segmentaremos los datos por estación. Esto nos permitirá observar las relaciones de estas variables climáticas a lo largo de diferentes épocas del año.
#Vamos a crear una función que nos facilite la creación de este tipo de gráficos en todo el trabajo
crear_grafico_correlacion <- function(data_frame, x, y, color) {
#Creamos dataframe para obtener los coeficientes de correlación
correlaciones <- data_frame %>%
group_by(!!sym(color)) %>%
summarize(r = cor(!!sym(x), !!sym(y)))
#Creamos el gráfico
graf <- ggplot(data_frame, aes_string(x = x, y = y, color = color)) +
geom_point(data = data_frame %>% select(-!!sym(color)), colour = "grey85", show.legend = FALSE) +
geom_point(show.legend = FALSE) +
stat_smooth(show.legend = FALSE, formula = 'y ~ x', method = "lm") +
scale_color_manual(values = colores_estaciones) +
facet_wrap(as.formula(paste0("~", color))) +
geom_text(data = correlaciones, aes(label = sprintf("r = %.2f", r), x = Inf, y = Inf), hjust = "right", vjust = "top",
fontface = "bold") +
guides(color= FALSE)
return(graf)
}#Creamos los gráficos de dispersión demanda-humedad
p1 <- crear_grafico_correlacion(SeoulBikeData_estudio_clima, "Humedad", "Num_Alquileres", "Estacion") +
ggtitle("Relacion Humedad y Demanda")
#Creamos los gráficos de dispersión demanda-velocidad del viento
p2 <- crear_grafico_correlacion(SeoulBikeData_estudio_clima, "Vel_viento", "Num_Alquileres", "Estacion") +
ggtitle("Relacion Velocidad del viento y demanda") +
xlim(0,6) +
ylab(NULL)
p1 + p2
En el caso de la humedad con esta separación por estación del año podemos observar mas claramente la relación con la demanda. De acuerdo a a los gráficos de dispersión, podemos observar como, cuando supera el 50% de humedad, esta variable comienza a tener una relación negativa con la demanda, especialmente en verano. Una posible explicación es que la alta humedad aumenta la sensación térmica en días calurosos y la disminuye en días fríos, desincentivando la demanda del servicio de alquiler de bicis.
Por su parte, en el caso del viento no se aprecia ningún tipo de relación directa en ninguno de los meses. Aunque pudiera parecer que podríamos eliminar esta variable a la hora de construir el modelo de machine learning, no vamos a hacerlo, ya que podría existir interacciones con otras variables que no hemos contemplado.
Hemos analizado gráficamente la relación entre las variables meteorológicas y la demanda de bicicletas. Sin embargo, el objetivo de este trabajo es crear un modelo que ayude a predecir la demanda con 24 horas de antelación, por lo tanto, tal como se justificó en la introducción, no podemos utilizar las condiciones climáticas del mismo día y hora de la predicción. Por lo tanto, debemos entrenar nuestro modelo con datos climatológicos con al menos 24 horas de antelación. Aunque hay diversas formas de transformar y utilizar estos datos, en este estudio nos centraremos en dos: el valor registrado 24 horas antes y el valor promedio entre 24-48 horas antes. Se han elegido estos dos datos porque lo esperable es, que la información mas reciente tenga mas influencia en el modelo. Queda abierta la posibilidad, en futuras ampliaciones del trabajo, de analizar la incorporación de otras ventanas temporales que puedan mejorar el modelo.
for (var in variables_clima) {
# Crear el nombre de la nueva columna para datos desplazados
new_col_name_24h_antes <- paste0(var, "_24h_antes")
# Desplazar los datos 24 filas hacia abajo
SeoulBikeData[[new_col_name_24h_antes]] <- lag(SeoulBikeData[[var]], 24)
# Crear el nombre de la nueva columna para promedio de las últimas 24 sesiones
new_col_name_prom_24_48h <- paste0(var, "_prom_24_48h")
# Desplazar datos 48 filas hacia abajo para comenzar desde el registro 48
shifted_data <- lag(SeoulBikeData[[var]], 48)
# Calcular el promedio de las sesiones entre 24 y 48 registros atrás usando rollmean de zoo
SeoulBikeData[[new_col_name_prom_24_48h]] <- rollmean(shifted_data, 24, align = "right", fill = NA)
}Analisis de la multicolinealidad
Antes de entrenar el modelo con todas las variables que hemos mencionado a lo largo del documento, debemos analizar si hay variables que estén creando una alta multicolinealidad. Para identificar aquellas variables que produzcan una alta multicolinealidad vamos a utilizar dos métricas:
Factor de Inflación de la Varianza (VIF, por sus siglas en inglés), es una métrica que se utiliza para detectar si hay multicolinealidad en modelos de regresión. Normalmente se considera que, si el valor está por encima de 10, se está creando una alta multicolinealidad.
Coeficiente de correlación: Esta métrica la utilizaremos de forma complementaria al VIF y nos servirá para identificar aquellos pares que pueden estar causando la multicolinealidad.
# Obtenemos los nombres de las nuevas columnas
Nombres_col_clima <- c(sapply(variables_clima, function(var) paste0(var, "_24h_antes")),
sapply(variables_clima, function(var) paste0(var, "_prom_24_48h")))
# Calcular la matriz de correlación para esas columnas
correlation_matrix <- cor(SeoulBikeData[, Nombres_col_clima], use = "pairwise.complete.obs")
# Filtrar las correlaciones absolutas mayores o iguales a 0.5, excluyendo la diagonal y pares repetidos
upper_tri <- which(upper.tri(correlation_matrix, diag = FALSE), arr.ind = TRUE)
high_correlations <- which(abs(correlation_matrix[upper_tri]) >= 0.5, arr.ind = FALSE)
# Crear un dataframe de las correlaciones altas
high_correlations_df <- data.frame(
Variable1 = rownames(correlation_matrix)[upper_tri[high_correlations, 1]],
Variable2 = colnames(correlation_matrix)[upper_tri[high_correlations, 2]],
Correlación = correlation_matrix[upper_tri][high_correlations]
) %>%
arrange(desc(abs(Correlación)))
#Calculamos el modelo der regresión lineal
model <- lm(Num_Alquileres ~ ., data = SeoulBikeData[, c("Num_Alquileres", Nombres_col_clima, "Demanda_24h_aj", "Demanda_Acum_24_53h")])
# Calcula el VIF para cada variable en el modelo
vif<- vif(model)
# Convertir vif a un dataframe para poder representarlo
vif_df <- data.frame(Variables = names(vif), VIF = as.numeric(vif), stringsAsFactors = FALSE) %>%
arrange(VIF) %>%
mutate(Variables = factor(Variables, levels = Variables))
# Crear el gráfico con valores de VIF
p1 <- ggplot(data = vif_df, aes(x = Variables, y = VIF)) +
geom_bar(stat = "identity", fill = "#1585E8", color="#004C90", alpha=0.6) +
coord_flip() +
theme_minimal() +
labs(title = "Valores VIF para las variables", y = "Valor de VIF", x = "Variables") +
geom_text(aes(label=round(VIF, 1)), hjust=-0.2) +
geom_hline(aes(yintercept = 1, linetype="Baja"), color="green") +
geom_hline(aes(yintercept = 3, linetype="Media"), color="orange") +
geom_hline(aes(yintercept = 10, linetype="Alta"), color="red") +
scale_linetype_manual(values=c("dashed", "dashed", "dashed"),
breaks = c("Baja", "Media", "Alta"),
guide = guide_legend(title = "Niveles de Multicolinealidad",
override.aes = list(color = c("green", "orange", "red")))) +
theme(legend.position="top")
# Ordenar el dataframe por correlación absoluta de mayor a menor
high_correlations_df <- high_correlations_df %>% arrange(Correlación)
# Convertir la combinación de Variable1 y Variable2 a un factor con niveles ordenados por correlación absoluta
high_correlations_df$VariablePair <- factor(high_correlations_df$Variable1 %>% paste("-", high_correlations_df$Variable2),
levels = high_correlations_df$Variable1 %>% paste("-", high_correlations_df$Variable2))
#Creamos gráfico con los pares de correlaciones más altos
p2 <- ggplot(data = high_correlations_df, aes(y = VariablePair, x = Correlación)) +
geom_bar(stat = "identity", fill = "#1585E8", color="#004C90", alpha=0.6) +
geom_text(aes(label=round(Correlación, 2)), hjust=1.2, vjust=0.5) +
theme_minimal() +
labs(title = "Correlaciones más altas entre pares de variables", x = "Valor de Correlación", y = "Pares de Variables") +
scale_x_continuous(limits = c(-1, 1)) + # Establecer los límites del eje X
geom_vline(aes(xintercept = 0), color = "black") + # Línea vertical negra en X = 0
geom_vline(aes(xintercept = 0.9), color = "red") + # Línea vertical roja en X = 0.9
geom_vline(aes(xintercept = -0.9), color = "red") # Línea vertical roja en X = -0.9
# Mostramos los gráficos
p1
p2
#Eliminamos las variables que ya no necesitamos
rm( upper_tri, high_correlations, high_correlations_df, model, vif)Tan solo hay un par de variables en el que podamos identificar un problema de alta colinealidad: “Temperatura_24_antes” y “Temperatura_prom_24_48h”, ambas con un VIF mayor a 10 puntos (11,6 y 12). Además, este par de variables tiene un coeficiente de correlación de 0.93, por lo que queda claro que es esta alta correlación lo que está causando la alta multicolinealidad y deberemos eliminar una de estas variables. Optaremos por eliminar la variable “Temperatura_prom_24_48h”, ya que puede ser menos intuitiva para el negocio. Al conservar la variable “Temperatura_24_antes”, facilitamos la comprensión del modelo para las partes interesadas.
# Eliminamos la variable del conjunto de datos y del vector con los nombres de las variables climáticas
SeoulBikeData$Temperatura_prom_24_48h <- NULL
Nombres_col_clima = setdiff(Nombres_col_clima, "Temperatura_prom_24_48h")Creación del modelo de predicción de Machine Learning
Para poder tener un conjunto datos apto para entrenar y validar un modelo que prediga la demanda de bicletas vamos a realizar las siguientes transformaciones:
Eliminamos los días es lo que el servicio no estuvo en funcionamiento, ya que ya sabemos que en estos días no va a haber demanda y no es necesario predecirlos. De esta forma simplificamos el modelo y evitamos distorsiones innecesarias.
Utilizaremos la técnica de one-hot encoding para los diferentes tipos de días: días laborables y sábados. Para evitar la multicolinealidad, hemos omitido la categoría de domingos/festivos. Por lo tanto, cuando ambas variables dummy (días laborables y sábados) sean 0, se entenderá que es un domingo o un día festivo.
Eliminamos las filas donde alguna variable no tenga los datos necesarios para su calculo, ya que algunas variables del modelo se derivan del promedio o acumulado de horas anteriores, por ejemplo, la demanda acumulada entre 24 y 53 horas previas. Si no se han registrado suficientes horas previas, estos valores no pueden ser calculados adecuadamente. Se ha optado por esta alternativa, ya que el número de datos que excluimos del dataset no es significativo (1 día y medio respecto a los 365 días).
Transformamos la hora en dos dimensiones usando funciones trigonométricas, de esta forma hacemos que horas como las 23:00 y las 0:00 estén mas cerca entre sí. Las dos variables que creemos solo se utilizarán para los modelos que asuman relaciones lineales. En los modelos que no asuman relaciones lineales, no utilizaremos estas dos variables y utilizaremos la hora original, ya que así evitamos que el modelo pierda interpretabilidad.
A continuación, dividiremos el dataset en dos subconjuntos: entrenamiento y evaluación, con el 75% y 25% de los datos respectivamente. Para garantizar que cada subconjunto represente adecuadamente las diferentes estaciones del año y eventos meteorológicos específicos como la lluvia y la nieve, optamos por una división estratificada. Esta decisión se basa en el análisis previo de los datos, donde observamos como la demanda era mucho baja en invierno que en el resto de meses del año y que, aunque poco frecuentes, la lluvia y la nieve tenían una clara influencia en la demanda.
Antes de crear un modelo de ML, vamos a crear tres modelos de referencia o “baselines” que nos permitan tener un punto de referencia para entender cuanto mejora la predicción un modelo más complejo. Los baselines que vamos a crear son:
Modelo Mediana: Predice la demanda basada en la mediana de unidades alquiladas para una hora específica y tipo de día en cuestión. Esta mediana se calcula a partir de los datos históricos disponibles para ese intervalo horario y tipo de día.
Modelo basado en demanda 24h antes: predice la demanda como el valor de la demanda que hubo para esa misma hora y tipo de día (valor que recoge la variable “Demanda_24h_aj”). Es decir, si queremos predecir la demanda de un lunes a las 14:00 tomará le valor del viernes (último dia laborable) a esa misma hora. La creación de este modelo se justifica por la alta correlación que vimos que había entre la demanda 24h antes y la demanda.
Modelo de regresión lineal: Este modelo incorpora todas las características (o variables independientes) con las que planeamos entrenar los modelos de ML más sofisticados. Aunque es más complejo que los dos anteriores, sigue siendo un baseline útil, ya que nos brinda una visión inicial sobre cómo las características se relacionan con la demanda.
Para evaluar el comportamiento tanto de los baselines como de los futuros modelos de machine learning utilizaremos dos métricas:
R-cuadrado: Indica la proporción de la variabilidad en la variable de respuesta que se explica por el modelo. Un valor cercano a 1 indica que el modelo explica una gran proporción de la variabilidad, mientras que un valor cercano a 0 sugiere lo contrario.
Error Absoluto Medio ( MAE, por sus siglas en inglés): Es el promedio de los errores absolutos entre las predicciones y las observaciones reales.
Error Absoluto Medio Porcentual (MAPE, por sus siglas en inglés): Es una métrica que cuantifica el error en términos porcentuales, lo cual puede ser útil para comprender el error en relación con la magnitud de las observaciones que se están prediciendo,
Error Absoluto Mediano Porcentual (MedAPE, por sus siglas en inglés): de calculo similar al MAPE pero utilizando la mediana, ya que el MAPE es muy sensible a los valores bajos y, de esta forma, lo complementamos con una métrica menos sensible a estos valores.
Aunque se podrían haber utilizado otras métricas como, por ejemplo, el error cuadrático medio, se han elegido estas tres porque ademas de ser representativas de la calidad de la predicción son fácilmente entendidas por el negocio, ya que están expresadas en términos porcentuales o en la misma unidad que la variable a predecir.
calcular_metricas <- function(real, pred) {
# Verificar que los vectores tienen la misma longitud
if (length(real) != length(pred)) {
stop("Los vectores deben tener la misma longitud")
}
# Calcular R-cuadrado
ss_tot = sum((real - mean(real))^2)
ss_res = sum((real - pred)^2)
r2 = 1 - (ss_res/ss_tot)
# Calcular MAE
mae = mean(abs(real - pred))
# Calcular MAPE
mape = mean(abs((real - pred) / real)) * 100
# Calcular MedAPE
medape = median(abs((real - pred) / real)) * 100
return(c(R2 = r2, MAE = mae, MAPE = mape, MedAPE = medape))
}Veamos los resultados de los tres baselines de acuerdo a las métricas que hemos definido:
#Obtenemos la demanda real
demanda_real <- validacion$Num_Alquileres
#Añadimos las métricas de los baselines
res_metricas <-
rbind( "Mediana" = calcular_metricas(demanda_real, pred_mediana),
"Demanda_24h" = calcular_metricas(demanda_real, pred_demanda24h),
"Lineal" = calcular_metricas(demanda_real, pred_lm)) %>%
as.data.frame()
#Creamos una tabla para representar los datos
res_metricas %>%
gt(rownames_to_stub = TRUE, auto_align = TRUE) %>%
fmt_number(
columns = "R2",
decimals = 2) %>%
fmt_number(
columns = setdiff(names(res_metricas), "R2"),
decimals = 0
) %>%
data_color(
columns = c("R2"),
fn = scales::col_numeric(
palette = c("red", "white", "green"),
domain = NULL
)
) %>%
data_color(
columns = setdiff(names(res_metricas), "R2"), # Todas las columnas excepto "R2"
fn = scales::col_numeric(
palette = c("green", "white", "red"),
domain = NULL
)
)| R2 | MAE | MAPE | MedAPE | |
|---|---|---|---|---|
| Mediana | 0.36 | 396 | 212 | 43 |
| Demanda_24h | 0.60 | 208 | 127 | 17 |
| Lineal | 0.68 | 229 | 155 | 25 |
De los tres baselines, el que peor resultados obtiene es el basado en las Medianas, por lo tanto, podremos eliminar este ya que con los otros dos nos basta para compararlos con los modelos que creemos. De acuerdo al valor del R-cuadrado, el baseline basado en un modelo lineal es el que mejor explica la variabilidad de la demanda, ya que es el que mayor coeficiente tiene (0.66). Por su parte, el baseline basado en Demanda 24h es el que menos se desvía de la predicción, de acuerdo al MAE y el MAPE se desvia de media 207 unidades y un 127% respecto al valor real. El alto valor del MAPE se puede deber a la alta sensibilidad a valores bajos, ya que el MedAPE tiene un desviación signficativamente menor, del 17%.
#Eliminamos el modelo basado en la mediana, ya que es el que peor se comporto en todas las métricas
res_metricas <- res_metricas %>%
subset(row.names(res_metricas)!= "Mediana")Después de establecer los modelos de baseline, nuestro siguiente paso es la creación de modelos predictivos iniciales. Estos modelos se configurarán con parámetros básicos o predeterminados en su primera iteración. Utilizaremos las métricas de evaluación previamente definidas para medir su desempeño. El modelo o modelos que demuestren un rendimiento superior en esta etapa preliminar se someterán a una optimización más intensiva, que incluirá una búsqueda exhaustiva de hiperparámetros.
La elección de los algoritmos se basarán en las características mas relevantes que hemos observado del conjunto de datos:
Dado que contamos con solo un año de datos históricos, tenemos una muestra limitada de ciertos eventos, como los días de lluvia o nieve, que influyen significativamente en la demanda.
La demanda varia por hora y tipo de día (laborable, sábado o domingo/festivo) tienen una alta incidencia en la demanda.
Hemos identificado relaciones no lineales entre algunas variables. Por ejemplo, la demanda aumenta con la temperatura hasta los 25ºC; sin embargo, al superar este umbral, la relación se vuelve negativa.
Contamos con 13 variables relacionadas con el clima (todas variables continuas), 2 vinculadas a la demanda anterior (tambien variables continuas) y 3 asociadas con el tiempo cronológico para entrenar el modelo ( Hora y las variables dummy de “sábado” y “laborable”).
Para la creación de los modelos preliminares elegiremos los siguientes algoritmos de machine learning:
Random Forest: Hemos elegido este modelo debido a la capacidad de manejar relaciones no lineales gracias a su estructura en árboles y la resisentencia al sobreajuste.Además,proporciona métricas de importancia de características, lo que facilita la interpretación del modelo.
SVM: Seleccionamos este modelo ya que puede capturar relaciones no lineales mediante el uso de kernels y suele funcionar muy bien con conjuntos de datos de tamaño moderado, como el nuestro.
XGBoost: Este modelo, al igual que el Random Forest, se basa en una estructura de árbol, lo que nos permite capturar relaciones no lineales. A priori, XGBoost es considerado una mejora respecto al Random Forest gracias a su regularización y optimización enfocada.
#Establecemos la semilla
set.seed(seed)
# RANDOM FOREST
# Entrenar el modelo de Random Forest
modelo_rf <- randomForest(formula_modelo_no_lin,
data = entrenamiento,
ntree = 25, # reducido para evaluación inicial
mtry = sqrt(ncol(entrenamiento[columnas_seleccionadas_lin])),
importance = TRUE)
# Asignamos las predicciones al modelo
pred_rand_forest <- predict(modelo_rf, newdata = validacion)
#SVM
# Entrenar el modelo SVM
modelo_svm <- svm(formula_modelo_no_lin, data = entrenamiento, kernel = "radial", cost = 1)
# Asignamos las predicciones al modelo
pred_svm <- predict(modelo_svm, newdata = validacion)
# XGBOOST
#creamos los datasets especificos para entrenar XGBoost
entrenamiento_xgb <- entrenamiento %>% ungroup() %>% select(Num_Alquileres, all_of(columnas_seleccionadas_no_lin) )
validacion_xgb <- validacion %>% ungroup()%>% select(Num_Alquileres, all_of(columnas_seleccionadas_no_lin) )
#Creamos las matrices que utilizaremos en el modelo XGBoost
dtrain <- xgb.DMatrix(data = as.matrix(entrenamiento_xgb[, -1]), label = entrenamiento_xgb[[1]])
dvalid <- xgb.DMatrix(data = as.matrix(validacion_xgb[, -1]))
# Parámetros para XGBoost
params <- list(
booster = "gbtree",
objective = "reg:squarederror",
eta = 0.3, # incrementado para evaluación inicial
gamma = 1,
max_depth = 3, # árbol más pequeño para evaluación inicial
min_child_weight = 1,
subsample = 0.8,
colsample_bytree = 0.8 # incrementado para usar más columnas
)
# Entrenar el modelo
num_rounds <- 50 # reducido para evaluación inicial
modelo_xgb <- xgb.train(params = params, data = dtrain, nrounds = num_rounds)
# Realizar predicciones con el modelo
pred_xgboost <- predict(modelo_xgb, newdata = dvalid)Veamos los resultados que han obtenido los modelos:
res_metricas <- res_metricas %>%
rbind( "Random Forest" = calcular_metricas(demanda_real, pred_rand_forest),
"SVM" = calcular_metricas(demanda_real, pred_svm),
"XGBoost" = calcular_metricas(demanda_real, pred_xgboost) )
res_metricas %>%
gt(rownames_to_stub = TRUE, auto_align = TRUE) %>%
fmt_number(
columns = "R2",
decimals = 2) %>%
fmt_number(
columns = setdiff(names(res_metricas), "R2"),
decimals = 0
) %>%
data_color(
columns = c("R2"),
fn = scales::col_numeric(
palette = c("red", "white", "green"),
domain = NULL
)
) %>%
data_color(
columns = setdiff(names(res_metricas), "R2"), # Todas las columnas excepto "R2"
fn = scales::col_numeric(
palette = c("green", "white", "red"),
domain = NULL
)
)| R2 | MAE | MAPE | MedAPE | |
|---|---|---|---|---|
| Demanda_24h | 0.60 | 208 | 127 | 17 |
| Lineal | 0.68 | 229 | 155 | 25 |
| Random Forest | 0.85 | 149 | 93 | 15 |
| SVM | 0.82 | 154 | 101 | 15 |
| XGBoost | 0.80 | 182 | 126 | 21 |
#Eliminamos el modelo basado en la mediana, ya que es el que peor se comporto en todas las métricas
res_metricas <- res_metricas %>%
subset(!(row.names(res_metricas) %in% c("Demanda_24h", "Lineal", "SVM", "XGBoost")))El algoritmo que mejor resultado ha obtenido, en todas las métricas, es Random Forest. Por lo tanto, tal como expusimos anteriormente, vamos a buscar la mejor combinación de parámetros. Para ello, en primer lugar vamos a utilizar la función TuneRF() para optmizar el número de variables a considerar en cada división (parámetro mtry) para 10,25, 50 y 100(ntree). Hemos acotado el tamaño del número de árboles a 100 para evitar un tiempo excesivo de computación, aunque si los resultados sugieren lo contrario, deberemos ampliar el número.
Posteriormente, evaluaremos el comportamiento en el conjunto de validación de cada conjunto de parámetros (mejor mtry para cada ntree) y escogeremos para la siguiente fase de optimización el que mejor puntuación haya obtenido en las métricas que hemos definido previamente. Se han elegido optimizar primero estas dos variables porque son las que más influencia suelen tener en los modelos de Random Forest.
En la siguiente fase, buscaremos el valor optimo para el tamaño minímo que tiene que tener un nodo para seguir dividiendo (nodesize). Para ello, evaluareamos el comportamiento del modelo escogido en el punto anterior con tamaños de nodesize entre 1 y 10.
Veamos en primer lugar las métricas de los distintos nº de árboles con el nº de variables en cada división óptimo:
res_metricas %>%
gt(rownames_to_stub = TRUE, auto_align = TRUE) %>%
fmt_number(
columns = "R2",
decimals = 2) %>%
fmt_number(
columns = setdiff(names(res_metricas), "R2"),
decimals = 0
) %>%
data_color(
columns = c("R2"),
fn = scales::col_numeric(
palette = c("red", "white", "green"),
domain = NULL
)
) %>%
data_color(
columns = setdiff(names(res_metricas), "R2"), # Todas las columnas excepto "R2"
fn = scales::col_numeric(
palette = c("green", "white", "red"),
domain = NULL
)
)| R2 | MAE | MAPE | MedAPE | |
|---|---|---|---|---|
| Random Forest | 0.85 | 149 | 93 | 15 |
| RF n_tree= 10 mtry= 5 | 0.84 | 158 | 95 | 16 |
| RF n_tree= 25 mtry= 5 | 0.85 | 149 | 94 | 15 |
| RF n_tree= 50 mtry= 6 | 0.86 | 144 | 90 | 14 |
| RF n_tree= 100 mtry= 5 | 0.86 | 143 | 94 | 14 |
Observamos, que hay dos modelos que sobresalen por encima del resto, el que tiene 50 árboles de decisión y el que tiene 100. Vamos elegir el modelo con 50 árboles por dos razones, en primer lugar que el modelo de 50 árboles es más sencillo, lo que lo vuelve más robusto a la sobreoptimización. En segundo lugar, aunque en 3 de las 4 métricas el resultado es mejor en el modelo de 100 árboles, esta diferencia es muy pequeña (inferior a la unidad), sin embargo la mejora en el MAPE es de 4 unidades (90% de error medio frente a un 94%).
Por lo tanto, este será el modelo que optimicemos en la siguiente, en la que analiceramos los resultados de las métricas para todos los tamaños mínimos de nodo entre 1 y 10. veamos los resultados de esta optimización:
#Eliminamos todos los modelos menos el escogido
res_metricas <- res_metricas %>%
subset((row.names(res_metricas) %in% "RF n_tree= 50 mtry= 6"))# Vector con los valores de nodesize
nodesize_values <- 1:10
# Bucle sobre los valores de ntree
for(nodesize in nodesize_values){
set.seed(seed)
# Ajuste del modelo randomForest con el mtry óptimo
modelo_rf_opt <- randomForest(formula_modelo_no_lin,
data = entrenamiento,
ntree = 50,
mtry = 6,
nodesize = nodesize,
importance = TRUE)
# Predicciones con el modelo optimizado
pred_rand_forest_opt <- predict(modelo_rf_opt, newdata = validacion)
# Agregar métricas al data frame de resultados
nombre_fila <- paste("RF n_tree= 50 mtry= 5 nodesize=", nodesize)
nuevo_registro <- as.data.frame(t(calcular_metricas(demanda_real, pred_rand_forest_opt)))
rownames(nuevo_registro) <- nombre_fila
res_metricas <- res_metricas %>%
rbind( nuevo_registro )
}res_metricas %>%
gt(rownames_to_stub = TRUE, auto_align = TRUE) %>%
fmt_number(
columns = "R2",
decimals = 2) %>%
fmt_number(
columns = setdiff(names(res_metricas), "R2"),
decimals = 0
) %>%
data_color(
columns = c("R2"),
fn = scales::col_numeric(
palette = c("red", "white", "green"),
domain = NULL
)
) %>%
data_color(
columns = setdiff(names(res_metricas), "R2"), # Todas las columnas excepto "R2"
fn = scales::col_numeric(
palette = c("green", "white", "red"),
domain = NULL
)
)| R2 | MAE | MAPE | MedAPE | |
|---|---|---|---|---|
| RF n_tree= 50 mtry= 6 | 0.86 | 144 | 90 | 14 |
| RF n_tree= 50 mtry= 5 nodesize= 1 | 0.86 | 145 | 93 | 14 |
| RF n_tree= 50 mtry= 5 nodesize= 2 | 0.86 | 146 | 95 | 14 |
| RF n_tree= 50 mtry= 5 nodesize= 3 | 0.86 | 145 | 90 | 14 |
| RF n_tree= 50 mtry= 5 nodesize= 4 | 0.86 | 144 | 91 | 14 |
| RF n_tree= 50 mtry= 5 nodesize= 5 | 0.86 | 144 | 90 | 14 |
| RF n_tree= 50 mtry= 5 nodesize= 6 | 0.86 | 148 | 93 | 14 |
| RF n_tree= 50 mtry= 5 nodesize= 7 | 0.86 | 146 | 95 | 14 |
| RF n_tree= 50 mtry= 5 nodesize= 8 | 0.86 | 148 | 95 | 14 |
| RF n_tree= 50 mtry= 5 nodesize= 9 | 0.85 | 149 | 96 | 14 |
| RF n_tree= 50 mtry= 5 nodesize= 10 | 0.85 | 149 | 98 | 14 |
En todas los tamaños mínimos del nodo entre 3 y 5 los valores son relativamente similares. Sin embargo, cuánto mayor sea el tamaño del nodesize más robusto será el modelo, por lo tanto elegiremos el modelo con tamaño mínimo de nodo igual a 5. En definitiva, el modelo que pondremos en producción será un Random Forest con 100 árboles, que considera 6 variables en cada división y necesita un número mínimo de 5 observaciones para realizar la división.
# Creamos el modelo con los parámetro óptimos
set.seed(seed)
modelo_rf_final <- randomForest(formula_modelo_no_lin,
data = entrenamiento,
ntree = 50,
mtry = 6,
nodesize = 5,
importance = TRUE)
#Insertamos la predicción en el dataframe de validación para poder evaluar los resultados
validacion <- validacion %>%
mutate(pred = predict(modelo_rf_final, newdata = validacion))Vamos a analizar el modelo resultante, en primer examinemos la importancia tiene cada variable en el Random Forest utilizando el Porcentaje de Incremento del Error Cuadrático Medio (% Inc. MSE). Esta métrica se basa en cuanto se deteriora el MSE cuando se permutan valores de una variable en particular. Veamos la importancia de las distintas variables en el modelo:
# Creamos el dataframe que recoge la importancia de variables
importancia_variables <- as.data.frame(importance(modelo_rf_final)) %>%
rename(Inc_MSE = !!names(.)[1]) %>%
arrange(Inc_MSE)
# Añadimos una columna con los nombres y los grupos para facilitar el trabajo
importancia_variables <- importancia_variables %>%
mutate(
Variable = factor(rownames(importancia_variables), levels = rownames(importancia_variables)),
Grupo = case_when(
Variable %in% Nombres_col_clima ~ "Clima",
Variable %in% Nombres_col_demanda ~ "Demanda",
Variable %in% Nombres_col_tiempo_no_lin ~ "Tiempo"
)
)
#importancia_variables$Variable <- factor(rownames(importancia_variables), levels = rownames(importancia_variables))
ggplot(data = importancia_variables, aes(x = Variable, y = Inc_MSE, fill = Grupo)) +
geom_bar(stat = "identity", color="black", alpha=0.7) +
coord_flip() +
theme_minimal() +
labs(title = "Importancias de las características", y = "% Inc. MSE", x = "Características") +
geom_text(aes(label=sprintf("%.0f%%", Inc_MSE)), hjust=-0.05) +
scale_y_continuous(labels = scales::percent_format(scale = 1)) +
scale_fill_manual(values = c("Clima" = "#00A9DA", "Demanda" = "#B920F7", "Tiempo" = "#D15301"))
Vemos que las características que más importancia tienen para el modelo
es la Hora y la Demanda a 24h ajustada, recordemos que esta variable era
la demanda que hubo a la misma hora para el mismo tipo de día
(laborable, sábado o domingo/festivo). Que estas variables sean las que
mas importancia tienen es algo lógico, si tenemos en cuenta que en el
análisis veíamos como la demanda variaba altamente dependiendo de la
hora y la alta correlación que había entre la demanda y la demanda a 24h
ajustada.
Por otra parte, las variables que menos importancia tienen son la cantidad de lluvia 24h antes y las relacionadas con la cantidad de nieve, todas con un 3% o menos de % Inc. MSE. Ante esta baja importancia podríamos pensar en eliminar estas variables del modelo, para asi tener un modelo más sencillo, sin embargo hay que tener en cuenta que la poca importancia de estas variables puede ser debido a la baja ocurrencia de días con lluvía o nieve. Además, en caso de eliminar estas variables estaríamos perdiendo información que, aunque poco frecuente, puede ser muy útil, especialmente en el caso de la nieve, que nos quedaríamos sin ninguna variable que reflejara si ha nevado.
Vamos a analizar el detalle, la distribución de los errores de predicción. Para estudiar el error, vamos a trabajar exclusivamente con la diferencia entre la predicción y la demanda real. Se ha decidido trabajar con este tipo de error porque es fácilmente entendible por el negocio y porque nos va a permitir ver si ciertos eventos hacen que el modelo sobrevalore o infravalore la demanda. En primer lugar vamos a analizar gráficamente como se distribuye el error dependiendo de las variables climáticas y de las horas. Para la representación gráfica se ha limitado el eje Y a valores entre -500 y 500 ya que en este rango ocurren la mayor parte de los errores y de esta manera los valores mayores a estos no nos dificultan la interpretación de los gráficos.
validacion <- validacion %>%
mutate(error = pred - Num_Alquileres,
error_abs = abs(pred - Num_Alquileres),
error_porc = pred/Num_Alquileres - 1,
error_porc_abs = abs(pred/Num_Alquileres - 1) )#Recogemos las variables climáticas en vector para facilitar el trabajo
variables_clima <- c("Cant_lluvia","Cant_nieve","Temperatura", "Humedad", "Vel_viento", "Visibilidad", "Rad_solar")
#Creamos el dataset
error_clima <- validacion %>%
select(c(error, variables_clima, Hora, Estacion))
#Representamos gráficamente los resultados
error_clima %>%
gather(-error, -Estacion, key = "var", value = "value") %>%
ggplot(aes(x = value, y = error)) +
facet_wrap(~ var, scales = "free") +
geom_point(aes(color = Estacion) ,alpha = 0.6, size = 1) +
stat_smooth(method = "loess", formula = "y ~ x") +
scale_color_manual(values = colores_estaciones) +
theme_minimal() +
theme(legend.position = "bottom", axis.title.x=element_blank()) +
ylim(-500,500) +
geom_hline(yintercept = 0, color = "black")
Vamos a analizar los resultados:
Cantidad de lluvia: vemos que en las horas que se produce lluvia el modelo tiende a sobrevalorar la demanda. Sabemos por el análisis previo que la lluvia reduce la demanda, sin embargo el modelo no ha sido entrenado con los datos del clima en esa hora, sino con datos de un mínimo de 24h de antelación. Aunque pueda haber variables con 24h de antelación que estén correlacionadas con la lluvia (Por ejemplo, si ha llovido 24h antes, la temperatura, etc.), el modelo difícilmente podrá anticiparse a la lluvia. Por lo tanto, en el caso de que llueva a esa hora, es esperable que la demanda sea sensiblemente inferior a la predicha por el modelo.
Cantidad de nieve: Ocurre un efecto similar a la lluvia aunque la magnitud es menor. Este se debe a que, como vemos en el gráfico las horas en las que nieva se acumulan en invierno, estación en la que la demanda es menor.
Hora: Ninguna hora tiene un sesgo positivo ni negativo. La amplitud de los errores es mayor en la horas de mayor demanda, algo lógico, ya que la posibilidad de error es mayor en una hora donde se alquilan miles de bicicletas que las 5 de la mañana cuando se alquilan decenas.
Humedad: El modelo sobrevalora la demanda a partir del 75% de humedad. Esto tiene relación con lo que veíamos en el análisis descriptivo, que la humedad a partir de este punto desincentivaba la demanda y, además, al igual que ocurre con la lluvia o la nieve, el modelo no tiene capacidad de predecir esta alta humedad.
Rad. Solar: No hay ninguna tendencia en el error, el signo del error es independiente de la radiación solar. Sin embargo, sin embargo, vemos que la amplitud del error es mayor cuando la radiación solar es este 0 y 1.
Temperatura: Hay cierta infravaloración de la demanda entre 20 y 25 ºC y sobrevaloración de la misma entre 35Cº y 40 Cº. En el análisis exploratorio veíamos como en los 30Cº se producía un punto de inflexión y la temperatura comenzaba a tener un efecto negativo, por lo que es razonable pensar que, en las temperaturas cercanas a esta se produzcan estas sobrevaloraciones e infravaloraciones de la demanda debido a cambios de la temperatura de un día a otro.
Velocidad del viento: No se observa ningún tipo de sesgo positivo o negativo ni una mayor amplitud de los errores en ninguna parte del eje X.
Visibilidad: vemos que hay cierta sobrevaloración cuando la visibilidad es baja, entre 0 y 1000. Además, cuanto mas cercana a 0 es la visibilidad mayor es el sesgo. De nuevo, un evento meteorológico difícilmente predecible por el modelo y que tiene una influencia significativa y negativa en la demanda hace que la predicción la sobrevalore.
Veamos ahora como si existe sobrevaloración o infravaloración de la demanda en las distintas estaciones y en los distintos tipos de día, para ello veamos como es las distribución de los errores por estación y por tipo de día, así como el promedio y mediana de los errores:
p1 <- ggplot(data = validacion, aes(x=error, color=Tipo_Dia)) +
geom_density(alpha=0.1, size=1) +
xlim(-500, 500) +
theme_minimal() +
ggtitle("Distribución errores por tipo de día") +
xlab("Dif. predicción y demanda real") +
geom_vline(aes(xintercept=0), color="black") +
theme(legend.position = "bottom")
p2 <- ggplot(data = validacion, aes(x=error, color=Estacion)) +
geom_density(alpha=0.05, size=1) +
xlim(-500, 500) +
theme_minimal() +
ggtitle("Distribución errores por Estación") +
ylab("Densidad") +
xlab("Dif. predicción y demanda real") +
scale_color_manual(values = colores_estaciones) +
geom_vline(aes(xintercept=0), color="black") +
theme(legend.position = "bottom")
print(p1 + p2)
#Calculamos los datasets que necesitamos para los graficos
calc_error_tipo_dia <- validacion %>%
group_by(Tipo_Dia) %>%
summarise(error_promedio = round(mean(error), 1),
error_median = round(median(error), 1 )
)
calc_error_estacion <- validacion %>%
group_by(Estacion) %>%
summarise(error_promedio = round(mean(error), 1),
error_median = round(median(error), 1)
)
# Convertir los dataframes de formato ancho a formato largo
calc_error_tipo_dia <- tidyr::gather(calc_error_tipo_dia, key = "metrica", value = "valor", -Tipo_Dia)
calc_error_tipo_estacion <- tidyr::gather(calc_error_estacion, key = "metrica", value = "valor", -Estacion)
#Grafico de erro por tipo de dia
p1 <- ggplot(data = calc_error_tipo_dia, aes(x = Tipo_Dia, y = valor, fill = metrica)) +
geom_bar(stat="identity", position="dodge", alpha=0.75, color="black") + # añade bordes negros a las barras
geom_hline(yintercept = 0, color = "black") + # añade una línea negra en y=0
geom_text(aes(label = sprintf("%.0f", valor),
y = ifelse(valor < 0, valor - 0.05, valor + 0.05),
color = ifelse(valor < 0, "red", "black")),
position = position_dodge(0.9), vjust = 1) +
labs(title="Errores por Tipo de Día", y="Error", x="Tipo de Día") +
scale_fill_manual(values=c("#10C6AD", "#D77901"),
name="Métrica",
breaks=c( "error_median", "error_promedio"),
labels=c( "Error Mediano","Error Promedio")) +
scale_color_identity(guide = "none") +
theme_minimal() +
theme(legend.position = "bottom")
#Grafico de error por estacion
p2<- ggplot(data = calc_error_tipo_estacion, aes(x = Estacion, y = valor, fill = metrica)) +
geom_bar(stat="identity", position="dodge", alpha=0.75, color="black") +
geom_hline(yintercept = 0, color = "black") +
geom_text(aes(label = sprintf("%.0f", valor),
y = ifelse(valor < 0, valor - 0.05, valor + 0.05),
color = ifelse(valor < 0, "red", "black")),
position = position_dodge(0.9), vjust = 1) +
labs(title="Errores por Estación", y="Error", x="Estación") +
scale_fill_manual(values=c("#10C6AD", "#D77901"),
name="Métrica",
breaks=c( "error_median", "error_promedio"),
labels=c( "Error Mediano","Error Promedio")) +
scale_color_identity(guide = "none") +
theme_minimal() +
theme(legend.position = "bottom")
#Imprimimos los gráficos
p1 + p2
Respecto a los errores por tipo de día, lo primero a destacar es que el modelo tiende a sobrevalorar la demanda para los domingos/festivos e infravalorar para los otros dos tipos de día (laborables y sábados). Además, en el caso de los sábados hay una gran diferencia entre la mediana (-5) y el promedio (-24) de los errores, lo que indica que a pesar de que en términos promedio a penas hay desviaciones particularmente grandes que están arrastrando el promedio hacia un valor más negativo.
Respecto a las desviaciones de la predicción por estaciones lo primero que debemos mencionar es que, la alta concentración de datos en torno al 0 se debe a que en invierno la demanda fue sensiblemente menor por lo que el error es normal que sea menor también. Por otra parte, en primavera la desviación de la demanda es mínima, lo que indica que no existe sesgo positivo ni negativo en las predicciones. No ocurre lo mismo en verano y en otoño, donde la predicción tienen a infravalorar la demanda, especialmente en verano.
Para terminar este apartado, me gustaría resaltar posibles ampliaciones del modelo que se podrían abordar en el futuro:
Incorporación de más datos históricos: actualmente el conjunto de datos del que disponemos solo tiene un año de históricos, lo cuál dificulta la generalización, ya que es posible que haya ciertos eventos de los que no disponemos de suficiente información (por ejemplo, días de lluvia o festivos).
Diferenciar el error de predicción: Se ha tratado de la misma manera un error de predicción que sobrestima o subestima la demanda. Estos errores pueden tener diferentes implicaciones para el servicio de alquiler de bicicletas. Por lo tanto, una posible ampliación de este trabajo puede consistir ponderar estos errores de forma diferente y entrenar un modelo que intente minimizar los errores ponderados.
Predicción geospacial: el modelo actual predice cuanta demanda de bicicletas va a haber, pero no donde. El servicio de alquiler puede tener decenas o cientos de puntos de servicio, por lo que una extensión valiosa puede ser calcular modelos para las distinas estaciones o zonas de las ciudad lo que permitiría una reposición más eficiente.
Incorporación de las predicciones meteorológicas: hemos señalado que una de las principales debilidades del modelo es que, cuando ocurren eventos climáticos que que desincentivan la demanda (por ejemplo lluvia o excesivo calor) tiende a sobrevalorar la demanda. Una posible de paliar estos es introducir las predicciones meteorológicas que se dieron 24h antes a la demanda que queramos predecir.
Inclusión de más variables retardadas: en el trabajo hemos estudiado un conjunto limitado de variables retardadas, orientadas, la mayoría de ellas, a un horizonte temporal entre 24-48h antes. Una posible via de estudio podría ser la inclusión de variables retardadas con otros horizontales temporales.
Informe Final
El objetivo del trabajo desarrollado en este documento ha sido proponer una solución al siguiente problema del equipo de logística de la ciudad de Seúl: determinar la cantidad adecuada de bicicletas en funcionamiento frente a una demanda fluctuante e incierta. Para abordar este problema se han abordado dos soluciones complementarias, por un lado, un análisis descriptivo de como es la relación de la demanda con el variables relacionadas con el tiempo (día, hora y mes), con el clima y con la propia demanda pasada. Por otra parte se ha desarrollado un modelo de machine learning para predecir la demanda con 24h de antelación.
Respecto al análisis descriptivo estas han sido las principales conclusiones:
- Patrón de la demanda muy diferente entre los días laborables y los fines de semana y festivos.
- Durante los días laborables se producen dos picos en la demanda, uno a las 8:00 y otro a las 18:00.
- Durante la semana la demanda está más repartida
- Existe una fuerte correlación entre la demanda actual y la demanda registrada a la misma hora durante el último día del mismo tipo (laborable, sábado o domingo/festivo). Por ejemplo, la demanda de un lunes a las 15:00 tiende a correlacionarse estrechamente con la demanda del viernes anterior a la misma hora, dado que ambos días son laborable.
- Hay una disminución significativa de la demanda los días que llueve o nieva.
- A partir del 75% la humedad tiene una influencia negativa en la demanda.
- Existe una relación positiva entre la radiación solar y la demanda.
- La temperatura incentiva la demanda hasta lo 25ºC. A partir de los 25ºC, comienza a desincentivar la demanda.
- La falta de visibilidad desincentiva la demanda.
En la construcción del modelo de Machine Learning se ha pasado por tres fases. En la primera se construyen tres modelos sencillos que no son de ML que nos servirán para comparar como funcionan los modelos de Machine Learning. En la segunda fase construyen tres modelos de Machine Learning, todos obtienen mejores métricas que los modelos base y, de estos modelos de ML el que mejor funciona es el de Random Forest. En la 3ª fase se optimizan los parámetros del modelo Random Forest, obtiendo el modelo final, que es el que puede utilizar el negocio para predecir la demanda.
El modelo de predicción ha obtenido las siguientes métricas:
- R2=0.86, esto implica que el modelo es capaz de explica el 86% de la variabilidad de la demanda.
- De media, comete un error de 144 unidades a la hora de predecir los datos en términos absolutos.
- En términos relativos el error medio es del 90%, pero este se debe a que los errores en valores bajos penalizan bastante este error en términos porcentuales. Por ejemplo, si el modelo predice una demanda de 110 unidades, pero debido a la lluvia se reduce a 10, el error es del 1000%. De hecho, el 50% de las predicciones tienen una desviación porcentual del 14% o menos.
Se ha analizado si la predicción tiende a sobrestimar o infraestimar la demanda dependiendo del clima, la hora, el tipo de día o la estación y estas han sido las principales conclusiones:
- Al no tener capacidad de predecir el clima, los eventos climático que tienen una alta incidencia en la demanda dificultan la capacidad de predicción del modelo, por ejemplo, como el modelo no tiene capacidad de predecir la lluvia, en caso de que llueva, el modelo tenderá a haber sobrestimado la demanda.
- El modelo tiende a sobrestimar la demanda en los domingos y festivos e infraestimarla en el resto de días.
- El modelo muestra una infraestimación de la demanda en los meses de verano y otoño.
En el Anexo “Productivización del modelo” se encuentra toda la información necesaria para poner el modelo en producción.
Anexos
Analisis exploratorio
En primer lugar vamos a cargar el conjunto de datos y comprobar los nombres de las columnas
SeoulBikeData <- read.csv("Datos Originales.csv", header=TRUE, fileEncoding="latin1")
print(colnames(SeoulBikeData))## [1] "Date" "Rented.Bike.Count"
## [3] "Hour" "Temperature..C."
## [5] "Humidity..." "Wind.speed..m.s."
## [7] "Visibility..10m." "Dew.point.temperature..C."
## [9] "Solar.Radiation..MJ.m2." "Rainfall.mm."
## [11] "Snowfall..cm." "Seasons"
## [13] "Holiday" "Functioning.Day"Vamos a acortar y traducir los nombres de las columnas para facilitar el trabajo que hagamos con el dataset.
Comprobemos que todas las variables tienen el tipo de variable que deberían tener:
str(SeoulBikeData)## 'data.frame': 8760 obs. of 14 variables:
## $ Fecha : chr "01/12/2017" "01/12/2017" "01/12/2017" "01/12/2017" ...
## $ Num_Alquileres: int 254 204 173 107 78 100 181 460 930 490 ...
## $ Hora : int 0 1 2 3 4 5 6 7 8 9 ...
## $ Temperatura : num -5.2 -5.5 -6 -6.2 -6 -6.4 -6.6 -7.4 -7.6 -6.5 ...
## $ Humedad : int 37 38 39 40 36 37 35 38 37 27 ...
## $ Vel_viento : num 2.2 0.8 1 0.9 2.3 1.5 1.3 0.9 1.1 0.5 ...
## $ Visibilidad : int 2000 2000 2000 2000 2000 2000 2000 2000 2000 1928 ...
## $ Punto_rocio : num -17.6 -17.6 -17.7 -17.6 -18.6 -18.7 -19.5 -19.3 -19.8 -22.4 ...
## $ Rad_solar : num 0 0 0 0 0 0 0 0 0.01 0.23 ...
## $ Cant_lluvia : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Cant_nieve : num 0 0 0 0 0 0 0 0 0 0 ...
## $ Estacion : chr "Winter" "Winter" "Winter" "Winter" ...
## $ Festivo : chr "No Holiday" "No Holiday" "No Holiday" "No Holiday" ...
## $ Dia_de_func : chr "Yes" "Yes" "Yes" "Yes" ...La variable “Fecha” no tiene asignado el tipo de variable correcto, por lo que corregiremos esto en el conjunto de datos.
# Convertimos la columna Fecha a tipo "Date"
SeoulBikeData$Fecha <- as.Date(SeoulBikeData$Fecha, format="%d/%m/%Y")En primer lugar vamos a comprobar el período temporal que abarca el dataset, para ello vamos a obtener la fecha mínima y máxima.
SeoulBikeData %>%
summarise(
fecha_min = min(Fecha, na.rm = TRUE),
fecha_max = max(Fecha, na.rm = TRUE)
)Tenemos un año completo de datos, desde el 1 de diciembre de 2017 hasta el 30 de noviembre de 2018. Esto puede ser una debilidad a la hora de construir el modelo, ya que solo tendremos información de los eventos que hayan ocurrido ese año. Por ejemplo, si ese año no ha nevado o las nevadas no fueron importantes no tendremos información de como se comportaría la demanda ante una nevada.
Vamos ahora a comprobar si el histórico es continuo o, si por el contrario, hay saltos entre las fechas. Para ello vamos a comprobar cual es la máxima diferencia entre fechas, en caso de que sea mayor a 1 significará que el histórico es discontinuo.
SeoulBikeData %>%
arrange(Fecha) %>%
mutate(diff = c(NA, diff(Fecha))) %>%
summarise(max_diff = max(diff, na.rm = TRUE))Vemos que el histórico es continuo, lo cual facilita bastante el trabajo ya que no tendremos que lidiar con saltos en el histórico.
La ultima comprobación que vamos a hacer relativa a las fechas es comprobar que todos los dias tienen 24 registros. Para ello vamos a calcular la fecha que tenga menos registros ya que, en caso de que este número fuera inferior a 24, significaría que hay por lo menos un día con datos incompletos.
SeoulBikeData %>%
group_by(Fecha) %>%
summarise(num_registros_hora = n()) %>%
summarise(minimo_registros_hora = min(num_registros_hora))La fecha que menos registros tiene es 24, por lo que podemos concluir que no hay fechas con datos incompletos.
Otra variable que merece la pena incluir en el análisis preliminar es, “Festivo”. Una pregunta que me gustaría responder es, ¿Qué está marcado como festivo en el dataset? Lo normal sería que solo las festividades, pero podría ser que los fines de semana o domingos estuvieran marcados también. Otra cuestión que podemos analizar es si existe algún mes con un número inusitadamente alto de festivos. Vamos a comprobar estas dos cuestiones creando un gráfico con el nº de festivos por mes.
#Creamos un dataset con los festivos que ha habido cada mes
festivos <- SeoulBikeData %>%
mutate(Year = year(Fecha), Month = month(Fecha)) %>%
group_by(Year, Month) %>%
summarise(num_festivos = sum(Festivo == "Holiday")/24, .groups = "drop") %>%
mutate(num_festivos = replace_na(num_festivos, 0))
#Creamos un gráfico para representar los resultados
ggplot(festivos, aes(x = interaction(Month, Year, sep = "-"), y = num_festivos)) +
geom_bar(stat = "identity", fill = "#8ddce5") +
geom_text(aes(label = num_festivos), vjust = -0.5) +
theme_minimal() +
labs(title = "Número de festivos por mes",
x = "Mes-Año",
y = "Número de festivos") +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
axis.text.y = element_blank())
# Eliminamos el dataset para ahorrar espacio en memoria
rm(festivos)Vemos que el mes con mas festivos tiene es diciembre con 4, un número mas alto que el resto pero no debería afectar a las predicciones. Además también vemos que no están incluidos los fines de semana, ya que si no habría muchos mas festivos.
LA última variable respecto a la fecha que voy a analizar es “Dia_de_func”. De acuerdo a la documentación esta variable recoge si ese día funcionaba el servicio o no.Comprobemos que, efectivamente, no hay un dato mayor que 0.
SeoulBikeData %>%
filter(Dia_de_func == "No") %>%
summarise(max_num_alquileres = max(Num_Alquileres))La conclusión es clara, tal como refleja la documentación, en estos días no existe servicio de alquiler de bicis, por lo tanto, a la hora de entrenar el modelo eliminaremos estos datos.
Veamos ahora como se distribuyen los Dias en los que no funciona el servicio por mes, con un doble motivo, por un lado, conocer si hay meses con un número muy alto de días sin funcionar.
#Creamos un dataset con los festivos que ha habido cada mes
no_func <- SeoulBikeData %>%
mutate(Year = year(Fecha), Month = month(Fecha)) %>%
group_by(Year, Month) %>%
summarise(dias_no_func = round(sum(Dia_de_func == "No")/24, 2 ), .groups = "drop") %>%
mutate(dias_no_func = replace_na(dias_no_func, 0))
#Creamos un gráfico para representar los resultados
ggplot(no_func, aes(x = interaction(Month, Year, sep = "-"), y = dias_no_func)) +
geom_bar(stat = "identity", fill = "#8ddce5") +
geom_text(aes(label = dias_no_func), vjust = -0.5) +
theme_minimal() +
labs(title = "Días sin funcionamiento por mes",
x = "Mes-Año",
y = "Nº Días") +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
axis.text.y = element_blank())
# Eliminamos el dataset para ahorrar espacio en memoria
rm(no_func)Vemos que no hay ningún mes con un número demasiado alto de días sin funcionar, además, también vemos que el número de días no es un número entero, lo que indica que hay días que tienen una horas etiquetadas como día no funcional y otras como funcional.
Vamos a comprobar como es la distribucción de los alquileres por hora:
psych:::histo.density(SeoulBikeData$Num_Alquileres, dcol = c('cadetblue','blue'), bcol='lightblue')
Sin necesidad de realizar ningún contraste, podemos concluir que la distribucción de los datos no es una distribucción normal- Deberemos tener esto en cuenta a la hora de realizar realizar los contrastes de hipótesis.
Por último, vamos a hacer un pequeño análisis de la correlación que tienen las variables climáticas entre sí, ya que es bastante probable que exista una alta correlación entre ellas.
# Selecciona las variables de interés
variables_clima <- c("Temperatura", "Humedad", "Vel_viento", "Visibilidad", "Punto_rocio", "Rad_solar","Cant_lluvia", "Cant_nieve")
correlaciones <- cor(SeoulBikeData[variables_clima], use="complete.obs")
# Creamos el Heatmap
correlaciones_melted <- melt(correlaciones)
ggplot(data = correlaciones_melted, aes(x=Var1, y=Var2)) +
geom_tile(aes(fill=value), colour = "white") +
geom_text(aes(label=round(value, 1)), size = 5) +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Correlation") +
geom_text(aes(label=ifelse(Var1 != Var2, round(value, 1), "")), size = 5) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, vjust = 1, size = 10, hjust = 1)) +
coord_fixed()
rm( correlaciones)Vemos que existe una alta correlación (0.9) entre el la temperatura y el punto de rocío, vamos a eliminar el punto de rocío para evitar introducir multicolinealidades en el modelo. La decisión de eliminar el punto de rocío y no la temperatura se basa en facilitar la compresión del modelo al negocio, ya que la temperatura es un parámetro más intuitivo y comprensible.
Productivización del modelo
Recordemos que uno de los objetivos del presente trabajo era que el servicio de alquiler de bicis pudiera tener un modelo que predijera la demanda de bicicletas. Por lo tanto, en esta sección nos encargaremos de realizar el código y funcionalidades necesarias para poder utilizar la el modelo para predecir nuevos datos.
El marco de trabajo del que partiremos es que el modelo recibirá como input un CSV con los datos de los días pasados y los nuevos datos que tiene que predecir y deberá producir como output un CSV con la predicción de la demanda.
En primer lugar vamos a exportar el modelo. De esta, el modelo se convierte en un objeto importable que cualquiera puede utilizar sin necesidad de entrenarlo.
save(modelo_rf_final, file="rf_model.RData")A continuación, vamos a reproducir todo el código necesario para poder importar un CSV y exportar un CSV con la prediccion. Para que el código funcione es necesario que en la carpeta donde el código estén los siguientes ficheros:
- CSV con los datos de los días anteriores y los que necesitamos predecir llamado “Anexo_Productivizacion.csv”
- Archivo RData con el modelo Random Forest llamado “Anexo_Modelo_RF.RData”.
Con el envío del presente documento se ha enviado el archivo “RData” y un ejemplo sintético de fichero CSV para la predicción.
library(dplyr)
library(tidyverse)
library(randomForest)
library(zoo)
#Cargamos el modelo y el CSV
load("Anexo_Modelo_RF.RData")
SeoulBikeData <- read.csv("Anexo_Productivizacion.csv", header=TRUE, sep = ";", fileEncoding="latin1")
#Cambiamos los nombres de las columnas
colnames(SeoulBikeData) = c("Fecha", "Num_Alquileres", "Hora", "Temperatura", "Humedad", "Vel_viento",
"Visibilidad", "Punto_rocio", "Rad_solar", "Cant_lluvia", "Cant_nieve",
"Estacion", "Festivo", "Dia_de_func")
#Creamos vectores con nombres de variables que nos facilitarán el código
variables_clima <- variables_clima <- c("Temperatura", "Humedad", "Vel_viento", "Visibilidad", "Rad_solar","Cant_lluvia", "Cant_nieve")
#Añadimos las variables con lag relacionadas con el clima
for (var in variables_clima) {
# Crear el nombre de la nueva columna para datos desplazados
new_col_name_24h_antes <- paste0(var, "_24h_antes")
# Desplazar los datos 24 filas hacia abajo
SeoulBikeData[[new_col_name_24h_antes]] <- lag(SeoulBikeData[[var]], 24)
# Crear el nombre de la nueva columna para promedio de las últimas 24 sesiones
new_col_name_prom_24_48h <- paste0(var, "_prom_24_48h")
# Desplazar datos 48 filas hacia abajo para comenzar desde el registro 48
shifted_data <- lag(SeoulBikeData[[var]], 48)
# Calcular el promedio de las sesiones entre 24 y 48 registros atrás usando rollmean de zoo
SeoulBikeData[[new_col_name_prom_24_48h]] <- rollmean(shifted_data, 24, align = "right", fill = NA)
}
#Transformamos los datos
SeoulBikeData <- SeoulBikeData %>%
replace_na(list(Num_Alquileres = 0)) %>%
mutate(Fecha = as.Date(Fecha, format="%d/%m/%Y")) %>%
select(-Punto_rocio) %>%
mutate(
Dia_Semana = weekdays(Fecha),
Mes = months(Fecha)
) %>%
mutate(Festivo = ifelse(Dia_Semana == "domingo", "Holiday", Festivo)) %>%
#Creamos variable de tipo de día
mutate(Tipo_Dia = ifelse(
Dia_de_func == "No", "No Funcional", ifelse(
Festivo == "Holiday", "domingo/festivo", ifelse(
Dia_Semana == "sábado", "sábado", "laborable"
)))) %>%
#Creamos la variable demanda 24h ajustada
group_by(Tipo_Dia) %>%
mutate(Demanda_24h_aj = ifelse(Tipo_Dia == lag(Tipo_Dia, 24), lag(Num_Alquileres, 24), NA)) %>%
ungroup() %>%
#Creamos la variable de demanda acumulada 24-53h
group_by(Tipo_Dia) %>%
mutate(Demanda_Acum_24_53h =
rollapply(Num_Alquileres, width = 53, FUN = sum, align = "right", fill = NA) -
rollapply(Num_Alquileres, width = 24, FUN = sum, align = "right", fill = NA) ) %>%
ungroup() %>%
#Aplicamos One Hot Encoding
mutate(
sabado = ifelse(Dia_Semana == "sábado", 1, 0),
laborable = ifelse(Tipo_Dia == "laborable", 1, 0)
) %>%
#Nos quedamos solo con los datos que queremos predecir
filter(Num_Alquileres == 0)
#Creamos un fichero con la predicción
prediccion <- SeoulBikeData %>%
mutate(prediccion = predict(modelo_rf_final, newdata = SeoulBikeData)) %>%
select(Fecha, Hora, prediccion)
#Guarda el fichero como CSV
write.csv(prediccion, file = "prediccion.csv", row.names = FALSE)Documento con código completo
Entre los archivos anexos a este documento se incluyen tanto el archivo “rmd” como un PDF con el documento con el código incluído.