Predicción del uso de la bicicleta pública en Donostia

Ander Fernández Jauregui.

Ander Fernández Jauregui

Data Scientist y Business Intelligence.

Problema Inicial

Suponemos que el área de movilidad del Ayuntamiento de Donostia/San Sebastián nos ha contactado para crear un modelo que prediga el uso de la bicicleta pública a un mes vista, para que puedan saber si necesitarán, o no, más bicicletas.

Asimismo, aunque no es prioritario, recalcan que les puede resultar muy interesante conocer la previsión de bicicletas a más corto-medio plazo para que puedan planificar las reparaciones de las mismas.

Cuestiones Involucradas

  • Machine Learning
    • ARIMA Univariante
    • ARIMA Multivariante
    • Random Forest
  • Integración de Datos
    • GeoJSON
    • API

Solución

Se ha creado un modelo de Machine Learning para cada uno de los enfoques que requería el Ayuntamiento. Los modelos elegidos han sido ARIMA Univariante (un mes vista), ARIMA Multivariante (una semana vista), Random Forest (a un día vista).

Para mejorar el rendimiento de los modelos, los datos se han enriquecido con información de la API de AEMET.

Desarrollo del Proyecto

Share on linkedin
Share on twitter
Share on email
Share on linkedin
Share on twitter
Share on email

Aviso. El caso expuesto a continuación es un caso ficticio, aunque se ha intentado simular lo más realista posible.

Business Understanding

Determinar el Objetivo de Negocio

Supongamos que el área de movilidad del Ayuntamiento de Donostia/San Sebastián nos ha contactado. Para poder conocer si el Ayuntamiento tendría que comprar más bicicletas de cara al mes que viene.

Fijar la situación

  • Datos: Nos comentan que disponen de sensores en todas las estaciones de bicicleta, que está accesibles mediante una API pública (la cual no facilitan). El sensor funciona desde hace ya varios años.
  • Enfoque: Aunque se podría fijar como un modelo de optimización que maximice la “disponibilidad” de las bicicletas basándose en el número de bicicletas y de aparcamientos en cada estación de bicis, el Ayuntamiento insiste en estimar la demanda futura.

Objetivo de Minería

  • Granularidad: necesitaremos los el número de viajes realizados en un día a nivel agregado.
  • Cadencia: necesitaremos los datos a nivel de día.

Data Understanding

Obtención de los datos

Como he comentado anteriormente, no disponemos de acceso a la API. Sin embargo, Donosti dispone de un panel abierto en el que cual se muestra, a tiempo real, el uso de las bicicletas.

Este dashboard seguramente esté “alimentado” por la API a tiempo real. Si inspeccionamos la página, vemos como uno de los ficheros Javascript que corren por detrás, hace una llamada a la API para obtener un JSON.

Si analizamos el resto del javascript, veremos como lanza una query con la siguiente estructura: http://www.donostia.eus/info/ciudadano/Videowall.nsf/movimientos.xsp? + fecha_inicial= + fecha en formato %Y-%M-%d%20HH:MM + fecha_final= + fecha en formato %Y-%M-%d%20HH:MM

Por ejemplo, la URL de los datos para el uso de bicicletas el 30 de junio del 2019 sería el siguiente http://www.donostia.eus/info/ciudadano/Videowall.nsf/movimientos.xsp?fecha_inicial=2019-06-30%2000:00&fecha_final=2019-06-30%2023:59. En él se encuentran los datos en modo de JSON.

Por tanto, para poder analizarlo, requereimos leer datos JSON desde una url, para lo cual usaremos la función fromJSON del paquete jsonlite. Más información aquí.

Veremos qué obtenemos al usar esa función sobre los datos del 30 de junio del 2019.

library(jsonlite)

url <- "http://www.donostia.eus/info/ciudadano/Videowall.nsf/movimientos.xsp?fecha_inicial=2019-06-30%2000:00&fecha_final=2019-06-30%2023:59"

datos <- fromJSON(url)
str(datos)
## 'data.frame':    422 obs. of  9 variables:
##  $ estacion_desenganche: chr  "16" "14" "8" "1" ...
##  $ fecha_desenganche   : chr  "2019-06-30 22:53:06" "2019-06-30 22:55:03" "2019-06-30 22:57:59" "2019-06-30 22:49:29" ...
##  $ estacion_enganche   : chr  "3" "4" "4" "12" ...
##  $ fecha_enganche      : chr  "2019-06-30 23:57:26" "2019-06-30 23:07:45" "2019-06-30 23:06:55" "2019-06-30 23:05:59" ...
##  $ idbici              : chr  "00000025" "00000164" "00000048" "00000028" ...
##  $ idrfid              : chr  "f428ec7d" "92dfe9b5" "d439ebdd" "d5be1b5e" ...
##  $ tipo_usuario        : chr  "ANUAL" "ANUAL" "ANUAL" "ANUAL" ...
##  $ edad                : chr  "29" "19" "22" "47" ...
##  $ codigo_postal       : chr  "20005" "20010" "20010" "20018" ...

Como vemos, la granularidad del origen de cada dato es cada uso, mientras que nosotros necesitamos una granularidad por día. Por tanto, tendremos que transformar todo lo que obtenemos para poder usarlo de forma más sencilla.

Por un lado, el número de trayectos es igual al número de observaciones de ese día. Siendo cada día un df, el número de trayectos es igual a dim(datos)[1].

Por tanto, simplemente podemos recopilar todos los datos y después agruparlos por día. En nuestro caso manejaremos datos de 5 años, lo que equivale a 5x365x400 = 730K registros, que puede ser asumible.

Realizamos un for para poder obtener y almacenar los datos. Tras varios intentos, puedes ver como:

  1. Si pides demasiadas fechas de golpe (de inicio a fin,por ejemplo), la API no funciona.
  2. Si pides una granularidad muy pequeña (a nivel de día, por ejemplo), si quieres analizar muchos días, tienes que hacer tantas peticiones como días, y el proceso se hace muy lento.

Por tanto, hemos decidido hacer peticiones de 100 días, para así agilizar al máximo el proceso de extracción de datos.

library(lubridate)

fecha_desde <- ymd("2014-06-30")
fecha_hasta <-  ymd("2019-06-30")
fechas <- seq(fecha_desde, fecha_hasta, 50)

Ahora, para cada una de esas fechas, obtenemos la url y apartir de esa url obtenemos los datos, que los iremos guardando en “datos”.

for(i in 1:length(fechas)){
  if(i == length(fechas)){
    j <- i
  }else{
    j <- i + 1 
  }
  
  url <- paste("https://www.donostia.eus/info/ciudadano/Videowall.nsf/movimientos.xsp?fecha_inicial=",fechas[i],"%2000:00&fecha_final=",fechas[j],"%2023:59",sep="" )
  if(i == 1){
    datos <- fromJSON(url, flatten = TRUE, simplifyDataFrame = TRUE)
    i <- i + 1
  } else{
    x <- fromJSON(url)
    datos <- rbind(datos,x)
    i <- i + 1
  }
   # descomentar si quieres ver el progreso al reproducirlo
   # print(i) 
}

Descripcción de los datos

Como podemos ver, tenemos un df con 1072127 observaciones y 9 variables. Estas variables son: estacion_desenganche, fecha_desenganche, estacion_enganche, fecha_enganche, idbici, idrfid, tipo_usuario, edad, codigo_postal

Estas variables se refieren a:

  • Estación desenganche: el número de la estación de la cual la bici ha sido desenganchada.
  • Fecha desenganche: cuándo esa bicicleta ha sido desenganchada.
  • Estación enganche: el número de la estación en la cual la bici ha sido enganchada.
  • Fecha enganche: cuándo la bicicleta ha sido enganchada.
  • idbici: el identificador de la bici que se ha usado.
  • idrfid: supongo
  • tipo_usuario: el tipo de suscripción que el usuario tiene contratado.
  • edad: edad del usuario.
  • código postal: código postal en el que reside el usuario.

Vamos a codificar correctamente las variables.

datos$estacion_desenganche <- as.factor(datos$estacion_desenganche)
datos$estacion_desenganche <- factor(datos$estacion_desenganche, levels = c(1:16))
datos$fecha_desenganche <- as.Date(datos$fecha_desenganche)
datos$estacion_enganche <- as.factor(datos$estacion_enganche)
datos$estacion_enganche <- factor(datos$estacion_enganche, levels = c(1:16))
datos$fecha_enganche <- as.Date(datos$fecha_enganche)
datos$idbici <- as.factor(datos$idbici)
datos$idrfid <- as.factor(datos$idrfid)
datos$tipo_usuario <- as.factor(datos$tipo_usuario)
datos$edad <- as.integer(datos$edad)
datos$codigo_postal <- as.factor(datos$codigo_postal)

Exploración de los datos

Analizaremos cada una de las variables una por una para tener una mejor comprensión de los datos con los que estamos tratanto.

Estación de enganche-desenganche

Cada estación estará distribuida en una zona geográfica diferente. Por tanto, es posible que no todos tengan la misma cantidad de enganches o desenganches.

library(ggplot2)
library(gridExtra)
library(dplyr)

datos %>%
  ggplot(aes(estacion_desenganche))+geom_bar(fill = "#00263e") + theme_minimal() + 
    labs(title = "Número de Desenganches por Estación", y = "Nº Casos", x = "Estación de Desenganche") +
    theme(panel.grid.major.x = element_blank()) 

 

datos %>%
  ggplot(aes(estacion_enganche))+geom_bar(fill = "#00263e") + theme_minimal() + 
    labs(title = "Número de Enganches por Estación", y = "Nº Casos", x = "Estación de Enganche") +
    theme(panel.grid.major.x = element_blank())





Como podemos ver, el número de enganches y desenganches cambia para cada estación. Esto indica que no todas las personas hacen trayectos de ida y vuelta o, al menos, no desde la misma estación donde engancharon la bicicleta.

Fechas de enganche/desenganche

Cada enganche se refiere a un servicio que se ha prestado y es de suponer que, en la mayoría de los casos, un enganche es derivado de un desenganche en el mismo día. Por tanto, analizando la evolución del número de enganches o desenganches (cuál analizar debería ser prácticamente indiferente), estamos analizando la evolución del uso del servicio de bicicletas.

datos %>%
  group_by(fecha_desenganche) %>%
  summarize(n_servicios = n()) %>%
  ggplot(aes(fecha_desenganche, n_servicios)) + geom_line(col =  "#8dc8e8") + geom_smooth(method ="lm", se=FALSE, col = "#00263e") +
    labs(title = "Evolución y Tendencia del uso de las bicicletas", x= "Fecha", y="Número de Desenganches") 

En ese sentido, podemos ver dos cuestiones claras:

  1. El uso de la bicicleta en Donostia tiene una clara tendencia positiva, tal como muestra la línea azul
  2. El uso de la bicicleta es estacional, es decir, se usa más en ciertas fechas que en otras. Esto es algo bastante lógico atentiendo al clima de la región, puesto que al tratarse de una región lluviosa, en aquellas épocas cuando llueve más, la bicicleta se usará menos.

Tipo de Usuario

Existen dos modalidades sobre las cuales puedes usar la bicicleta. Por un lado, la suscripción anual del servicio, que van desde los 40euros el abono anual, hasta los 8euros por el abono del día. Existen, además, diferentes abonos denominados “ocasionales”. Fuente.

Es de esperar que los residentes sean los usuarios del abono anual, mientras que los abonos ocasionales sean utilizados por turistas.

datos %>%
  ggplot(aes(tipo_usuario, fill = tipo_usuario)) + geom_bar() +
    labs(title = "Distribución según tipo de usuario", x = "Tipo de Usuario", y = "Nº de casos", fill = "") +
    theme(panel.grid.major.x = element_blank(), legend.position = "bottom") +
    scale_fill_manual(values=c("#00263e", "#8dc8e8"))

Como podemos ver, la mayoría de los usuarios de las bicicletas en Donostia son residentes en la ciudad. Esto es un aspecto importante, puesto tiene una gran influencia a la hora de hacer predicciones, ya que nos indica qué variables no deberíamos incorporar al modelo.

Edad

datos %>%
  filter(!is.na(edad)) %>%
  group_by(edad, idrfid) %>%
  summarize(n_personas = n()) %>%
  ggplot(aes(x =edad)) + geom_bar() +
    labs(title = "Distribución de la edad de los usuarios", x = "Edad", y = "Nº de usuarios")

Como podemos ver, la edad parece tener ciertos valores extraños, puesto que es imposible tener una edad superior a los 1.500 años (y menos seguir andando en bicicleta con esa edad ;). En el siguiente apartado “Calidad del Datos” nos encargaremos de analizar estos casos. Por el momento, filtraremos a los usuarios para ver las edades inferiores a 100 años.

datos %>%
  filter(!is.na(edad), edad<100) %>%
  group_by(edad, idrfid) %>%
  summarize(n_personas = n()) %>%
  ggplot(aes(x = edad)) + geom_bar() +
    labs(title = "Distribución de la edad de los usuarios (filtrado)", x = "Edad", y = "Nº de usuarios") +
    scale_x_continuous(breaks = seq(0,100,5))

Vemos como la mayoría de los usuarios son jóvenes de entre 15 y 25 años o personas adultras de entre 50 y 55. Asimismo, podemos observar como hay casos de niños muy jóvenes (menos de 5 años), e incluso niños con edad negativa, lo cual muestra un error en los datos.

Código Postal

Puede ser interesante analizar cómo es geográficamente el uso de la bicicleta de una manera visual, puesto que nos puede dar ciertas características interesantes sobre la misma. Para ello, nos hemos descargado los datos geoespaciales en formato GEOJSon de GitHub, a través del cual podremos visualizar esta información.

library(rgdal)
cp <- readOGR("Recursos\\GUIPUZCOA.geojson")
## OGR data source with driver: GeoJSON 
## Source: "C:\Users\Ander\Documents\Big Data\Proyectos\ML\Bicis Donosti\bicis\Recursos\GUIPUZCOA.geojson", layer: "GUIPUZCOA"
## with 171 features
## It has 4 fields
cp <- cp[(cp$COD_POSTAL %in% c(20001:20018)),] #Filtramos por Donostia
plot(cp, main = "Datos geoespaciales de Donostia/San Sebastián")

Ahora separamos esta información y la añadimos a la BBDD que ya disponemos.

library(spbabel)

geo <- sptable(cp)[,c(1,5,6)]

means <- aggregate(geo[,2:3], geo[1], mean) #Obtenemos el centroide de los poligonos
colnames(means) <- c("ID","lat","lon")
means$ID <- cp$COD_POSTAL


cp@data$lat <- means$lat
cp@data$lon <- means$lon

Ahora añadimos la información agrupa por cada código postal al geoJSON. Como hay CP que aparecen en más de una geometría, haremos aparecer esos valores más de una vez.

x <- datos %>% group_by(codigo_postal) %>% summarize(n_viajes = n()) %>% filter(codigo_postal %in% c(20001:20018))
x <- rbind(x, c(20017,103),c(20017,103),c(20017,103),c(20017,103),c(20017,103),c(20014,734),c(20003,669),c(20006,677))
cp@data$numero_viajes <- as.vector(x$n_viajes)

Ahora podemos hacer un mapa interactivo con la librería leaflet.

library(leaflet)

pal <- colorNumeric(
 palette = "Blues",
 domain = cp$numero_viajes)

leaflet() %>%
  addProviderTiles(providers$Esri.NatGeoWorldMap) %>%
  addPolygons(data=cp, color = ~ pal(numero_viajes), fillOpacity = 1, stroke=FALSE) %>%
  addMarkers(cp$lat,cp$lon, label = as.character(cp$numero_viajes),labelOptions = labelOptions(noHide = F))

Puedes ver el mapa interactivo aquí.

Analizando el mapa vemos como hay una zona donde el uso de la bicicleta es casi nulo, concentrándose, principalmente en dos códigos postales. Sin conocer la horografía de la ciudad, seguramente, se trate de la zona neurálgica de la ciudad.

De todos modos, el uso está bastante expandido, por lo que el código postal no parece que vaya a resultar un buen indicador para predecir el uso de la bicicleta.

Calidad del Dato

Valores Perdidos

Vamos a analizar si hay algún caso en el que el sensor haya fallado y haya resultado en un valor perdido (NA).

colSums(is.na(datos))[colSums(is.na(datos))>0]
## estacion_desenganche    estacion_enganche                 edad 
##                   31                  106                33846 
##        codigo_postal 
##                38922

Como podemos ver, tenemos valores perdidos principalmente para el código postal y la edad y algunos para la estación de enganche o desenganche.

Rspecto al primero, seguramente se trate de casos en los que el billete es ocasional y, por tanto, no se dispone de esos datos. Lo comprobamos.

datos %>%
  filter(is.na(edad) | is.na(codigo_postal)) %>%
  distinct(idrfid, tipo_usuario) %>%
  group_by(tipo_usuario) %>%
  summarize(casos = n())
#tipo_usuariocasos
1ANUAL216
2OCASIONAL7867

Vemos como no es así y hay bastantes casos en los que ususarios con abono anual no tienen o bien la edad o el código postal. Veámos qué porcentaje sobre el total suponen estos números.

datos %>%
  distinct(idrfid, tipo_usuario) %>%
  group_by(tipo_usuario) %>%
  summarize(n_ususarios = n())
#tipo_usuarion_usuarios
1ANUAL6445
2OCASIONAL7884

Vemos como en el caso del ocasional la falta de datos es del 99,78%, ya que en el proceso de obtención del abono no se requieren esos datos. En el caso de los usuarios anuales la completitud del dato es del 96.64%. Se trata por tanto de un dato que, en caso de ser significativo, podríamos usar.

Respecto a los datos de enganche o desenganche, para ambos casos son muy pocos datos los que faltan y, en cualquier caso, como requerimos de datos agregados a nivel de día, no es indispensable contar con ambos. Siplemente nos vamos a asegurar que no haya observaciones donde ambos son nulos.

datos %>%
  filter(is.na(estacion_desenganche),is.na(estacion_enganche)) %>%
  summarize(casos = n())

#casos
131
Coherencia interna
  1. Edad del usuario

Como hemos visto en el análisis, hay ususarios con edades poco normales: gente con más de 100 años y con menos de 0 años. Analicemos cuántos casos son:

datos %>%
  filter(edad >100 | edad < 0) %>%
  distinct(edad,idrfid, edad) %>%
  summarize(cantidad = n())

#cantidad
59

En 59 usuarios la edad está incorrectamente añadida. Esto implica una completitud de la edad del 99.70%, por lo que se trata de una variable que nos podría servir.

  1. Duración del trayecto

Tendría poco sentido que para un mismo trayecto dejes la bicicleta antes de cogerla, es decir, que la fecha de enganche sea anterior a la de desenganche. Analicemos si se da el caso.

datos %>%
  mutate(duracion_trayecto = fecha_enganche - fecha_desenganche) %>%
  filter(duracion_trayecto<0) %>%
  summarize(cantidad = n())

#cantidad
1982

Como podemos ver, en 982 ocasiones el sensor funcionó mal. Aunque parezca mucho, al dispones de muchos datos la importancia es relativa. De hecho, completitud de las fechas de enganche y desenganche del 99.90%.


Conclusión: para todas las variables disponemos de una completitud suficiente con la que poder realizar un análisis. Ahora, toca elegir aquellas variables que de verdad tengan una relación con el número de viajes que se realizan.

Data Preparation

Construcción de datos

Las variables más sencillas que podemos construir son aquellas referentes a la fecha. Mediante el paquete lubridate, podremos obtener cuestiones tales como el día de la semana, el día del mes o el mes para posteriormente analizar si son variables que influyen en el uso de la bicicleta.

datos$mes <- as.factor(month(datos$fecha_desenganche))
datos$dia_mes <- as.factor(day(datos$fecha_desenganche))
datos$dia_semana <- as.factor(wday(datos$fecha_desenganche))
datos$semana <- as.factor(week(datos$fecha_desenganche))
datos$fecha <- ymd(datos$fecha_desenganche)

Integración y limpieza de datos

Hay dos datos principales que podemos integrar:

1. Clima: seguramente dispongamos de datos climatolóicos históricos (precipitaciones, viento, etc.) una granularidad suficientemente buena como para poder utilizarla.

Para poder obtener los datos climáticos, utilizaremos la (API de AEMET)[https://opendata.aemet.es/centrodedescargas/inicio]. Los datos que obtendrás son un fichero sh. jsonlite no puede leer estos archivos directamente, así que simplemente los he descargado y convertido a formato .JSON con Notepad.

  • Nota: AEMET no deja descargar datos de más de 5 años, por lo que tendrás que hacer dos peticiones, juntar los dos JSON e importarlos.
clima <- fromJSON("Recursos\\datos_clima.json")
str(clima)
## 'data.frame':    1825 obs. of  20 variables:
##  $ fecha      : chr  "2014-06-30" "2014-07-01" "2014-07-02" "2014-07-03" ...
##  $ indicativo : chr  "1024E" "1024E" "1024E" "1024E" ...
##  $ nombre     : chr  "DONOSTIA/SAN SEBASTIÁN, IGUELDO" "DONOSTIA/SAN SEBASTIÁN, IGUELDO" "DONOSTIA/SAN SEBASTIÁN, IGUELDO" "DONOSTIA/SAN SEBASTIÁN, IGUELDO" ...
##  $ provincia  : chr  "GIPUZKOA" "GIPUZKOA" "GIPUZKOA" "GIPUZKOA" ...
##  $ altitud    : chr  "251" "251" "251" "251" ...
##  $ tmed       : chr  "16,1" "16,4" "15,6" "17,2" ...
##  $ prec       : chr  "0,0" "11,8" "9,5" "18,7" ...
##  $ tmin       : chr  "12,0" "13,9" "13,7" "16,0" ...
##  $ horatmin   : chr  "04:00" "19:10" "04:20" "03:30" ...
##  $ tmax       : chr  "20,2" "19,0" "17,4" "18,5" ...
##  $ horatmax   : chr  "13:30" "08:40" "18:40" "16:05" ...
##  $ dir        : chr  "02" "29" "30" "32" ...
##  $ velmedia   : chr  "2,8" "3,1" "3,1" "1,7" ...
##  $ racha      : chr  "6,9" "10,8" "7,8" "5,6" ...
##  $ horaracha  : chr  "17:00" "15:03" "09:25" "00:05" ...
##  $ sol        : chr  "12,2" "1,4" "0,0" "0,0" ...
##  $ presMax    : chr  "991,5" "985,7" "992,0" "991,7" ...
##  $ horaPresMax: chr  "00" "23" "23" "00" ...
##  $ presMin    : chr  "983,7" "981,5" "984,7" "986,7" ...
##  $ horaPresMin: chr  "24" "06" "04" "24" ...

Nos deshacemos de todas las variables que no nos interesan (datos geográficos y temas relacionadas con las horas).

borrar <- c("indicativo","nombre","provincia","altitud","horatmin","horatmax","horaracha","horaPresMax","horaPresMin")
clima <- clima[,!colnames(clima) %in% borrar]

Por último, agregamos los datos del uso de la bicicleta a nivel de día e incorporamos la información meteorológica.

datos <- datos %>%
  group_by(fecha) %>%
  summarize(n_viajes = n())

Volvemos a introducir los datos que ya habíamos generado.

datos$mes <- as.factor(month(datos$fecha))
datos$dia_mes <- as.factor(day(datos$fecha))
datos$dia_semana <- as.factor(wday(datos$fecha))
datos$semana <- as.factor(week(datos$fecha))

Finalmente juntamos los datos que hemos obtenido.

clima$fecha <- as.Date(clima$fecha)
datos <- inner_join(datos, clima, by="fecha")

2. Festivos: puede ser que el uso de la bicicleta cambie en función de si un día es o no festivo o, por ejemplo, si hay una festividad para incentivar el uso de la bicicleta.

De cara a los festivos, según se indica en el calendario laboral son: Año Nuevo, Jueves Santo, Viernes Santo, Pascua y lunes de Pascua, el Día Internacional de los Trabajadores, Santiago Apóstol, Asunción de la Virgen, la Fiesta Nacional de España, Día de Todos los Santos, el Día de la Constitución y la Navidad.

De estos días, solo aquellos relacionados con la Semana Santa ocurren en fechas diferentes cada año, puesto que esta festividad ocurre con el primer plenilunio de primavera.

Hacemos pues, festivo todos aquellos días que siempre caen en las mismas fechas:

datos$festivo <- NA
datos$temporal <- paste(month(datos$fecha),"-",mday(datos$fecha), sep="")
festivos <- c("1-1","1-5","25-7","15-8","12-10","1-11","6-12","25-12")
datos[datos$temporal %in% festivos,"festivo"] <- 1
datos[!datos$temporal %in% festivos,"festivo"] <- 0

datos$temporal <- NULL

Por último, vamos a ver si los valores de la serie temporal pueden tener sentido y es que puede que el número de bicicletas no haya sido el miso a lo largo del tiempo y, por tanto, haya habido cambios sustanciales en el uso de la bicicleta. Para ello, mostramos el gráfico de la evolución del uso de la bicicleta.

ggplot(datos, aes(fecha,n_viajes)) + geom_line(col = "#00263e") +
  labs(title = "Evolución del uso de las bicis", x= "Tiempo",y="Nº diario de usos") 

Como podemos ver, a finales de 2015 el uso de la bicicleta descendió de manera significativa. Veamos qué ocurrió en ese periodo.

datos %>%
  filter(fecha>as.Date("2015-03-30"), fecha<as.Date("2016-01-01")) %>%
  ggplot(aes(fecha,n_viajes)) + geom_line(col = "#00263e") +
  labs(title = "Evolución del uso de las bicis (filtrado)", x= "Tiempo",y="Nº diario de usos")

Como no conocemos las razones por las que ocurre este fenómeno, vamos a eliminar estos datos de la serie.

datos <- datos %>%
  filter(fecha>as.Date("2015-12-31")) 

Formato de Datos

Finalmente, asignaremos a cada dato su formato correspondiente. Para ello, primero cambiaremos las comas de las variables numéricas de clima por puntos, puesto que es el formato en el que están el resto de datos e interpreta R.

datos$tmed <- gsub(",",".", datos$tmed, fixed = TRUE)
datos$prec <- gsub(",",".", datos$prec, fixed = TRUE)
datos$tmin <- gsub(",",".", datos$tmin, fixed = TRUE)
datos$tmax <- gsub(",",".", datos$tmax, fixed = TRUE)
datos$velmedia <- gsub(",",".", datos$velmedia, fixed = TRUE)
datos$racha <- gsub(",",".", datos$racha, fixed = TRUE)
datos$sol <- gsub(",",".", datos$sol, fixed = TRUE)
datos$presMax <- gsub(",",".", datos$presMax, fixed = TRUE)
datos$presMin <- gsub(",",".", datos$presMin, fixed = TRUE)

Ahora, procedemos a asignar el formato correcto a cada una de las variables.

datos$tmed <- as.numeric(datos$tmed)
datos$prec <- as.numeric(datos$prec)
## Warning: NAs introducidos por coerción
datos$tmin <- as.numeric(datos$tmin)
datos$tmax <- as.numeric(datos$tmax)
datos$dir <- as.integer(datos$dir)
datos$velmedia <- as.numeric(datos$velmedia)
datos$racha <- as.numeric(datos$racha)
datos$sol <- as.numeric(datos$sol)
datos$presMax <- as.numeric(datos$presMax)
datos$presMin <- as.numeric(datos$presMin)

Siguiendo el mensaje del código anterior, revisamos las variables que contengan NAs. Estas se habrán creado debido a la falta de datos de una observación para esa variable.

colSums(is.na(datos))
##      fecha   n_viajes        mes    dia_mes dia_semana     semana 
##          0          0          0          0          0          0 
##       tmed       prec       tmin       tmax        dir   velmedia 
##          0         49          0          0          0          0 
##      racha        sol    presMax    presMin    festivo 
##          0          1          0          0          0

Así pues, convertimos los NAs en 0.

datos$prec[is.na(datos$prec)] <- 0
datos$sol[is.na(datos$sol)] <- 0

Selección de los datos

Ahora tendremos que elegir aquellos datos que sean más sifnificativos a la hora de explicar el número de viajes en bicicleta.

Para ello,realizaremos dos análisis:

  • Un Random Forest para conocer la importancia de cada una de las variables, ya sean numéricas o categóricas.
  • Una tabla de correlaciones entre las variables numéricas para poder elegir variables sin que haya multicolinearidad.

Importancia de Variables con Random Forest

library(randomForest)

#Quitamos el día para poder hacer el análisis
temporal <- datos[,!colnames(datos) %in% "fecha"]

modelo <- randomForest(y =  temporal$n_viajes, x = temporal[,!colnames(temporal)%in% "n_viajes"], , ntree = 500, importance = TRUE)

importancia <- importance(modelo)
importancia <- data.frame(Variables = row.names(importancia), MSE = importancia[,1])
importancia <- importancia[order(importancia$MSE, decreasing = TRUE),]

ggplot(importancia, aes(x=reorder(Variables, MSE), y = MSE, fill=MSE)) + geom_col() + coord_flip() +
  labs(title = "Importancia de las Variables", x="variable",y="% incremento del MSE si la variable se permuta al azar") +
  theme(legend.position =  "none", panel.grid.major.x = element_blank())

Como podemos ver, la variable más importante es el día de la semana, seguido de si hay o no precipitaciones o si hace sol o no, la temperatura y el viento que hace. Como veoms, varias de esas variables se refieren a lo mismo (temperatura, viento, presión…), por lo que seguramente estén correlacionados.

Por tanto,** antes de asegurar con qué variables nos quedamos, vamos a realizar un análisis de multicolinearidad** entre las variables numéricas.

Análisis de Multicolinearidad
library(psych)

factor <- select_if(datos, is.factor)
numerico <- datos[, !colnames(datos) %in% c(colnames(factor),"fecha")]

correlaciones <- cor(numerico)

cor.plot(correlaciones, numbers=TRUE, xlas = 2, upper= FALSE, main="Correlación entre variables", zlim=c(abs(0.65),abs(1)))

rm(correlaciones, numerico, factor, clima)

Como podemos ver, la temperatura media está linearmente muy corelacionada tanto con la temperatura mínima (0.96) y con la temperatura máxima (0.97), y a la vez, estas dos también están muy correlacioneadas (0.87). Lo mismo ocurre con las rachas de viento y la velocidad media del viento (0.82) y la presión máxima y mínima (0.91).

Por tanto, teniendo esto en consideración, las variables numéricas con las que nos quedamos son: día de la semana, precipitación, sol, temperatura media, dir, presión máxima, racha y semana.

Modeling

A la hora de predecir el número de viajes, existen dos enfoques posibles:

  • Serie temporal: nos permite predecir un periodo de tiempo de más allá de un día. La serie temporal puede adoptar dos formas:
  • Univariante: únicamente disponemos de una variable para predecir el uso de la bicicleta.
  • Multivariante: disponemos de más de una variable para predecir el uso e la bicicleta.
  • Modelo de regresión: se trata de estimar el uso que tendrá un día en base a un modelo de regresión al uso (RandomForest, XGBoost, etc.). Estos modelos serán validos en tanto y cuanto las predicciones meteorológicas que tengamos sean válidas. Lógicamente, su precisión depende, a su vez, de la precisión de las predicciones meteorológicas.

En nuestro caso siguiendo las indicaciones del cliente (Ver el apartado de Business Understanding), realizaremos los dos modelos de serie temporal. Además, probaremos un modelo Random Forest para ver cómo funcionan a la hora de predecir a corto plazo, por si puediera llegar a ser de interés.

ARIMA Univariante

Creamos un dataframe únicamente con los datos que nos interesan.

datos1 <- datos %>%
  select(n_viajes, fecha)

ggplot(datos1, aes(fecha,n_viajes)) + geom_line(col = "#8dc8e8") +
  labs(title = "Evolución del número de viajes", x="fecha", y="Nº Viajes") +
  geom_smooth(method ="lm", se=FALSE, col = "#00263e")

Ahora, dividimos los datos entre train y test, dejando todos los datos para el train y el último mes para el test.

train <- datos1[(1:1221),]
test <- datos1[(1222:1251),]

Los modelos ARIMA se deben usar en series temporales estacionarias. Por tanto, lo primero que quebemos hacer es convertir nuestra serie temporal en estacionaria.

Comprobamos primero que la serie no es estacionaria con el test aumentado de Dickey-Fuller.

library(tseries)
adf.test(train$n_viajes, k = 1)
## Warning in adf.test(train$n_viajes, k = 1): p-value smaller than printed p-
## value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train$n_viajes
## Dickey-Fuller = -18.739, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationary

Como vemos, rechazamos la hipótesis nula (h0 = no es estacionaria), por lo que la serie es estacionaria y ya podemos aplicar los modelos ARIMA.

Para elegir los mejores valores de p y q para la parte Autoregresiva (AR) y de Medias Móviles (MA) respectivamente, vamos a usar el modelo autoarima. Sin embargo, para poder entender mejor la serie y asegurarnos de los resultados del modelo autoarima, vamos a realizar las autocorrelaciones (ACF) y autocorrelaciones parciales (PACF).

par(mfrow=c(1,2))

Acf(train$n_viajes, lag.max = 50, main="")
Pacf(train$n_viajes, lag.max = 50, main = "")

Vemos como a las autocorrelaciones les “cuesta” entrar en los valores fiajdos en el ACF, por lo que suponemos el valor de p será alto. Respecto a las autocorrelaciones parciales en lags pequeños ya “encajan” bien, por lo que suponemos que el vaor de q será bajo. Con estos datos, realizamos ajustamos el mejor modelo ARIMA.

modelo <- auto.arima(train$n_viajes, max.p=20,max.q=20)
modelo
## Series: train$n_viajes 
## ARIMA(6,1,0) 
## 
## Coefficients:
##           ar1      ar2      ar3      ar4      ar5      ar6
##       -0.5947  -0.6144  -0.5839  -0.5457  -0.5567  -0.3543
## s.e.   0.0268   0.0274   0.0286   0.0285   0.0274   0.0268
## 
## sigma^2 estimated as 60301:  log likelihood=-8443.52
## AIC=16901.04   AICc=16901.13   BIC=16936.79

Como podemos ver, el sistema considera que el ARIMA que mejor se ajusta (para un p y q menores de 20) tiene una p de 10 y una q de 3, lo cual concuerda con lo que habíamos analizado en las funciones ACF y PACF. Por tanto, damos por bueno el modelo.

Con este modelo, realizaremos la predicción y compararemos cuánto acierta o se equivoca.

prediccion <- forecast(modelo, h = 30)
y_ARIMA <- as.vector(prediccion$mean)
plot(prediccion, main = "Ajuste de la Predicción a 30 periodos vista",
     col = "#00263e", xlim = c(1200, 1251))
lines(ts(datos1$n_viajes), col = "red", lty = 2, lwd = 1.5)
legend("topleft", lty=1, col=c("#00263e","red"),
  legend=c("Predicción","Valor Real"))

Como podemos ver, las predicciones están bastante suavizadas. Sin embargo, casi todos los valores caen dentro del intervalo de confianza del 95%, lo cual es bueno. En definitiva, no acertaremos el número exacto para ningún día, pero podremos predecir la demanda de un mes sin fallar demasiado. Veamos cúal es exactamente el accuracy del modelo.

accuracy(prediccion,test$n_viajes)
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  1.224991 244.8571 180.5669 -17.98829 38.05493 0.7802769
## Test set     50.833724 279.2218 221.3275 -45.35971 73.30516 0.9564141
##                    ACF1
## Training set 0.01249985
## Test set             NA

Resultado: El modelo tiene un RMSE de 279.2218.

Modelo Multivariante

En este caso, intentaremos predecir el número de viajes, no solo usando la fecha que hemos obtenido, sino también el tiempo, en este caso, las precipitaciones. Para ello, vamos a realizar un modelo ARIMA multivariante.

Lo primero que debemos hacer es asegurarnos de que la serie es estacionaria. Esto ya lo hemos comprobado para la fecha en el modelo anterior, por lo que únicamente deberemos comprobarlo con la variable.

adf.test(train$prec, k = 1)
## Warning in adf.test(train$prec, k = 1): p-value smaller than printed p-
## value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train$prec
## Dickey-Fuller = -19.566, Lag order = 1, p-value = 0.01
## alternative hypothesis: stationary

Como vemos, rechazamos la hipótesis nula, por lo que se trata de una serie estacionaria. Lo comprobamos gráficamente.

ggplot(train, aes(x = fecha, y = prec)) + geom_line(col = "#8dc8e8") + geom_smooth(method = "lm", se=FALSE, col = "#00263e") + 
  labs(title = "Evolución y tendencia de las precipitaciones", x="", y = "Cantidad de precipitación")

Una vez nos hemos asegurado de tratarse de una serie estacionaria, deberemos calcular el Vector de Auto Regresión (VAR), según el cuál, cada variable es una función linear de sus propios valores pasados y el resto de valores pasados de otras variables. Y es que, el VAR es capaz de entender y usar la relación entre las diferentes variables. Para ello, usaremos la librería vars.

library(vars)
x <- train[, colnames(train) %in% c("n_viajes","prec")]

VARselect(x, lag.max = 10, type = "both")
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      9      9      7      9 
## 
## $criteria
##                   1            2            3            4            5
## AIC(n) 1.535426e+01 1.535286e+01 1.535615e+01 1.535879e+01 1.535438e+01
## HQ(n)  1.536695e+01 1.537188e+01 1.538152e+01 1.539050e+01 1.539243e+01
## SC(n)  1.538795e+01 1.540339e+01 1.542353e+01 1.544300e+01 1.545543e+01
## FPE(n) 4.658777e+06 4.652240e+06 4.667591e+06 4.679909e+06 4.659313e+06
##                   6            7            8            9           10
## AIC(n) 1.522452e+01 1.505737e+01 1.504663e+01 1.503673e+01 1.504219e+01
## HQ(n)  1.526892e+01 1.510811e+01 1.510370e+01 1.510014e+01 1.511195e+01
## SC(n)  1.534242e+01 1.519212e+01 1.519821e+01 1.520516e+01 1.522746e+01
## FPE(n) 4.091922e+06 3.462062e+06 3.425064e+06 3.391326e+06 3.409901e+06

En este caso, la función VarSelect nos indica el orden de lag para cada uno de los cuatro citerios: Akaike Information Criteria, Hannan-Quinn, Schwartz y Akaike’s Final Predicton Error.

Seguiremos el criterio AIC, puesto que la varianza de los errores que produce suelen ser menores que en el resto de los casos. (Fuente) Así pues, elegiremos un VAR con lag de 9.

var = VAR(x, p=9)
summary(var)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: n_viajes, prec 
## Deterministic variables: const 
## Sample size: 1212 
## Log Likelihood: -12512.882 
## Roots of the characteristic polynomial:
## 0.955 0.9484 0.9484 0.9017 0.9017 0.7607 0.7607 0.7158 0.7158 0.6884 0.6884 0.6845 0.6845 0.647 0.647 0.6468 0.6468 0.2563
## Call:
## VAR(y = x, p = 9)
## 
## 
## Estimation results for equation n_viajes: 
## ========================================= 
## n_viajes = n_viajes.l1 + prec.l1 + n_viajes.l2 + prec.l2 + n_viajes.l3 + prec.l3 + n_viajes.l4 + prec.l4 + n_viajes.l5 + prec.l5 + n_viajes.l6 + prec.l6 + n_viajes.l7 + prec.l7 + n_viajes.l8 + prec.l8 + n_viajes.l9 + prec.l9 + const 
## 
##              Estimate Std. Error t value Pr(>|t|)    
## n_viajes.l1   0.27101    0.03093   8.761  < 2e-16 ***
## prec.l1      -6.38460    0.85071  -7.505 1.20e-13 ***
## n_viajes.l2   0.02840    0.03182   0.893 0.372268    
## prec.l2       1.26542    0.88344   1.432 0.152300    
## n_viajes.l3   0.01662    0.02939   0.566 0.571698    
## prec.l3      -0.84284    0.87750  -0.961 0.336998    
## n_viajes.l4   0.01351    0.02877   0.469 0.638815    
## prec.l4       0.08679    0.87218   0.100 0.920755    
## n_viajes.l5  -0.04732    0.02882  -1.641 0.100959    
## prec.l5      -1.46841    0.87151  -1.685 0.092270 .  
## n_viajes.l6   0.20809    0.02890   7.201 1.05e-12 ***
## prec.l6       1.42956    0.87203   1.639 0.101404    
## n_viajes.l7   0.41626    0.02950  14.110  < 2e-16 ***
## prec.l7       5.25857    0.87179   6.032 2.16e-09 ***
## n_viajes.l8   0.04326    0.03216   1.345 0.178866    
## prec.l8       3.22423    0.88400   3.647 0.000276 ***
## n_viajes.l9  -0.12284    0.02990  -4.109 4.24e-05 ***
## prec.l9      -0.43805    0.87469  -0.501 0.616605    
## const       112.85835   32.31180   3.493 0.000495 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 230.8 on 1193 degrees of freedom
## Multiple R-Squared: 0.4912,  Adjusted R-squared: 0.4835 
## F-statistic: 63.99 on 18 and 1193 DF,  p-value: < 2.2e-16 
## 
## 
## Covariance matrix of residuals:
##          n_viajes   prec
## n_viajes  53249.0 -720.5
## prec       -720.5   71.4

Al tener un lag elevado, no haría falta analizar la cointegración mediante el test de Johansen.En cualquier caso, realizaremos el test de causalidad de Granger para comprobar si la precipitación es buena a la hora de estimar el número de viajes.

grangertest(train$n_viajes ~ train$prec, order = 9)

## Granger causality test
## Model 1: train$n_viajes ~ Lags(train$n_viajes, 1:9) + Lags(train$prec, 1:9)
## Model 2: train$n_viajes ~ Lags(train$n_viajes, 1:9)
##     Res.Df Df      F    Pr(>F)
## 1   1193                        
## 2   1202 -9 14.216 < 2.2e-16 ***
## --
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Rechazamos la hipótesis nula, por lo que se puede afirmar que los lags de prec explican las variaciones en el número de viajes. Por tanto, hacemos finalmente la predicción.

prediccion <- predict(var, n.ahead = 30, ci = 0.95, dumvar = NULL)

Comparamos los valores predichos con los reales de una forma gráfica, para ver cómo se ajusta la predicción.

y_MARIMA <- prediccion$fcst$n_viajes[,1]
y <- cbind(test, prediccion = y_MARIMA)

cols <- c("Prediccion"="#00263e","Real"="#8dc8e8")
ggplot(y, aes( x = fecha)) + 
  geom_line(aes(y=prediccion, col = "Prediccion"), size=2) + 
  geom_line(aes(y=n_viajes, col = "Real"), size=2) +
  labs(title = "Viajes reales vs Viajes Predichos", x = "") +
  scale_colour_manual(name="",values=cols) +
  theme(legend.position = "bottom")

Ahora, para poder comparar cuál de los dos modelos es mejor, estimaremos el RMSE. Al tratarse de un modelo que no se ha elaborado con la librería forecast, no podemos obtener el RMSE de la misma manera, por lo que lo calculamos de forma manual.

y %>%
  mutate(ER = (n_viajes-prediccion)^2) %>%
  summarize(RMSE = sqrt(sum(ER)/30))
## RMSE     
## 262.5904 

Random Forest

Vamos a estimar otro modelo de predicción, más allá de las dos series temporales. Este modelo, aunque no sirvan para poder predecir la demanda de todo el mes, puede que nos sirva para predecir mejor la demanda del día en función de variables meteorológicas que disponemos a priori.

datos1 <- datos[, !colnames(datos) %in% "fecha"]
train <- datos1[(1:1221),]
test <- datos1[(1222:1251),]

tuneRF(x = train[,-which(names(train) == "n_viajes")], y = train$n_viajes, tepFactor = 1.5, improve = 1e-5, ntreeTry = 500)
## mtry = 5  OOB error = 48203.06 
## Searching left ...
## mtry = 3     OOB error = 48946.85 
## -0.01543031 1e-05 
## Searching right ...
## mtry = 10    OOB error = 49693.63 
## -0.03092274 1e-05

##    mtry OOBError
## 3     3 48946.85
## 5     5 48203.06
## 10   10 49693.63

Vemos cómo el Out-Of-Bag Error (OOB) es menor cuando mtry = 4. Por tanto, utilizamos este parámetro para realizar el Random Forest.

modelo <- randomForest(n_viajes ~. , data = train, mtry = 4, ntree=1000)
test_predictors <- test[,-which(names(test) == "n_viajes")]
prediccion <- predict(modelo, test_predictors)
y_rf <- prediccion

Tras haber entrenado el modelo y predecir el valor para los próximos días, analizaremos cómo de bien se ha comportado el modelo, mediante un RMSE.

test <- cbind(test, prediccion)
test %>%
  mutate(ER = (n_viajes-prediccion)^2) %>%
  summarize(RMSE = sqrt(sum(ER)/30)) 
## RMSE     
## 199.8421

Elección del Modelo

Tras haber entrenado varios modelos para los mismos datasets de train y test hemos obtenido los siguientes resultados:

ModeloRMSE
ARIMA Univariante279.22
ARIMA Multivariante262.59
Random Forest199.92
library(tidyr)
resultados <- as.data.frame(cbind(ARIMA = y_ARIMA, MARIMA = y_MARIMA,  Random_Forest = y_rf, Valor = test$n_viajes))
resultados$fecha <- datos$fecha[1222:1251]

cols <- c("ARIMA"="#80ceff","MARIMA"="#40a3d8","Random Forest"="#1e6b95", "Valores" = "#00263e")
ggplot(resultados, aes(x = fecha)) +
  geom_line(aes(y = ARIMA, col = "ARIMA"), size = 1) +
  geom_line(aes(y = MARIMA, col = "MARIMA"), size = 1) + 
  geom_line(aes(y = Random_Forest, col = "Random Forest"), size = 1) +
  geom_line(aes(y = Valor, col = "Valores"), size =2, alpha = 0.8) +
  labs(title = "Comparación del Ajuste de los Modelos", x = "", y="Nº Desplzamaientos")+
  scale_colour_manual(name="",values=cols) +
  theme(legend.position = "bottom") 

En este caso, el modelo Random Forest bien tuneado estima mejor la demanda de bicicletas que los modelos ARIMA Multivariante y Univariante. Sin embargo, Random Forest requiere de diversos datos meteorológicos para poder hacer sus estimaciones, por lo que no podría ser utilizado a un mes vista.

Por tanto, la solución final pasaría por:

  1. Usar el modelo de Random Forest para predecir la demanda a corto plazo (un día para otro, por ejemplo).
  2. Usar el modelo ARIMA Multivariante para predecir la demanda a corto-medio plazo (una semana o 14 días en función de hasta dónde dispongamos de previsión meteorológica de lluvia)
  3. Usar el modelo ARIMA Univariante para predecir la demanda a largo plazo (a un mes vista).

Conclusión

Son dos las principales conclusiones que podemos obtener de este análisis:

  1. El enriquecimiento de datos resulta indispensable para poder mejorar la estimación de los modelos. Precisamente, los dos modelos con mejor precisión han sido aquellos que han sido entrenados con datos externos. En este sentido, el Análisis Descriptivo (EDA) resulta imprescindible.
  2. A la hora de elegir un modelo es indispensable atender a las restricciones de los modelos entrenados. Más allá de la capacidad explicativa del modelo, hay que ser conscientes del propósito que el modelo cumple para la organización y, por tanto, no basar la elección únicamente en el RMSE. Un modelo “perfecto” a la hora de predecir, si no cumple con las necesidades de la organización, no sirve para nada.