BiciMAD

Análisis y predicción de la demanda del servicio de público de bicicletas eléctricas de la ciudad de Madrid

Ezequiel Paura
2019-01-14

Tabla de contenidos

Este documento es un reporte técnico sobre la demanda de bicicletas eléctricas en la ciudad de Madrid.

Objetivos

  1. Análisis

Explorar las variables de interés para entender la demanda del servicio público BiciMAD.

  1. Predicción

Predecir esta demanda a nivel día y hora.

Datos

Se parte de dos conjuntos de datos, uno con información sobre las estaciones de bicicletas eléctricas y otro con la descripción de los trayectos realizados por usuarios de este servicio en 28 días del mes de septiembre de 2018.

De acuerdo con literatura reciente sobre el tema (e.g. El-assi et al,20171; Garcia-Palomares et-al, 20122), el factor climático juega un papel fundamental en la demanda de bike sharing. Los días de lluvia, por ejemplo, es razonable esperar una menor demanda del servicio. Por esta razón se complementan los datos iniciales con información climatológica.

Por un lado los de la Agencia Estatal de Meteorología OpenData disponibles al público a través de una API REST. Si bien los datos que consumiría un modelo predictivo deplegado en producción serán los pronósticos del clima, el registro histórico podría funcionar como un buen proxy. Por otro, se extraen datos climáticos de mayor granularidad del sitio Windguru para complementar los datos de la Agencia de Meteorología.

Luego de una limpieza inicial, se combinan todas las fuentes de datos.

Análisis exploratorio

Estaciones

El dataset con información sobre las estaciones consta de 172 filas y 15 variables. Cada fila contiene información sobre una estación (dirección, ubicación geográfica, cantidad de plazas, etc.).

Para entender la distribución geográfica del servicio, el mapa a continuación muestra la localización de cada estación, junto a la cantidad de plazas disponibles y los trayectos totales realizados desde la misma.

Figure 1: Estaciones de BiciMAD, número de plazas y trayectos totales realizados desde la misma (click sobre los íconos)

El siguiente mapa de calor muestra como se distribuyen la totalidad de los trayectos en cada estación. A mayor intensidad, más trayectos.

Las estaciones más activas. Distribución de densidad de la solicitud o devolución de bicicletas por estación

Figure 2: Las estaciones más activas. Distribución de densidad de la solicitud o devolución de bicicletas por estación

El máximo de trayectos realizados desde una estación es 6291 y 594 el mínimo. En promedio, se realizan 2369 trayectos al mes.

Figure 3: Ranking de las estaciones de acuerdo a la demanda en el mes de septiembre. Al pasar el ratón por cada punto se muestra el ID de la estación y su total de trayectos

Demanda de bicicletas

La variable de mayor interés de este análisis es la demanda de bicicletas por día y hora. El histograma a continuación presenta una distribución de frecuencia de ésta variable. En general, la cantida de trayectos es menor los fines de semana (con una excepción en la zona de la Embajada de Estados Unidos).

Distribución de la demanda de bicicletas por día y hora. La línea vertical indica la media.

Figure 4: Distribución de la demanda de bicicletas por día y hora. La línea vertical indica la media.

Se pueden ver dos picos pronunciados alrededor de 560 y de 1000 bicicletas por día/hora, y uno menos pronunciado alrededor de 1400 bicicletas por día/hora.

Es de esperar que la demanda varíe según el día de la semana. El siguiente gráfico presenta la demanda de cada día, con una marcada caída de la actividad los días sábado y domingo. Aunque menor, el viernes también parece un día de más baja actividad.

Figure 5: Demanda por día de la semana. Pasando el ratón por encima, se presentan la mediana, el primer y el tercer cuartil, máximos y mínimos

Estas variaciones se ven también a lo largo del día. A las horas pico (de 7 a 9, y de 17 a 21 horas) la demanda se incrementa notablemente. Un incremento menor se da de 14 a 16 horas.

Figure 6: Demanda por hora del día. Pasando el ratón por encima, se presentan la mediana, el primer y el tercer cuartil, máximos y mínimos

En los datos se identifican 3 tipos de usuarios:

Esta variable, como cualquier otra relacionada al usuario (e.g. edad), no podrá ser usada en la predicción ya que incurriríamos en el error de utilizar información obtenida ex-post al fenómeno a predecir. Éste problema se conoce como data leakage.

De todas formas puede ser interesante, en esta etapa exploratoria, saber si existen diferencias en la demanda de bicicletas entre los usuarios, en especial entre los usuarios con pase anual y los ocasionales. Como vemos en el gráfico a continuación, los patrones de uso a lo largo del día son consistentes para todos los tipos de usuarios.

Demanda por hora del día, por tipo de usuario.

Figure 7: Demanda por hora del día, por tipo de usuario.

Modelos predictivos

Para predecir la demanda de bicicletas por día y hora, partiremos los datos en dos partes. El 75% de los datos seleccionados de manera aleatoria, se utilizarán para entrenar un modelo, el 25% restante se utilizará para testearlo. Así y con sucesivas particiones de los datos durante el entrenamiento, evitamos un sobreajuste (overfitting) del modelo a los datos.

Métricas de evaluación

Para predecir el problema de regresión planteado, utilizaremos dos métricas que darán cuenta de cuán alejados estan los valores predichos de los “verdaderos”. La primera, y la más comunmente usada, es el Error Absoluto Promedio (o MAE, por sus siglas en inglés). El MAE mide la magnitud promedio del error, sin tener en cuenta la dirección del mismo (i.e. si es una sobre o una infra estimación). Y se calcula como el promedio de las diferencias (absolutas) entre cada valor “verdadero” y el valor predicho. O en notación matemática:

\[MAE = \frac{1}{n} \sum_{j = 1}^{n} \mid y_j - \hat y_j \mid \] Otras de las medidas utilizadas habitualmente para medir el error de una regresión, es la raíz del error cuadrático medio (o RMSE por sus siglas en inglés). El RMSE es una regla de puntuación cuadrática que mide la magnitud promedio del error. Es la raíz cuadrada del promedio de las diferencias cuadradas entre la predicción y la observación real.

\[RMSE = \sqrt{\frac{1}{n} \sum_{j = 1}^{n} (y_j - \hat y_j)^2} \]

Tanto el MAE como el RMSE expresan el error promedio de la predicción en la misma unidad que la variable de interés. Ambas pueden oscilar entre 0 y infinito, y son indiferentes a la dirección del error. En ambos casos, valores más pequeños indican mejor performance.

La diferencia entre ambas es que en el RMSE, al elevar al cuadrado los errores antes de calcular la media, se da relativamente más peso a los errores más grandes. Cuando los errores de gran magnitud son particularmente dañiños, el RMSE será más útil que el MAE. La principal ventaja del MAE, es su interpretabilidad, ya que permite ser leído como “el error medio”.

Para el entrenamiento de los modelos, la medida a optimizar será el RMSE.

No obsante, se incluirá en el reporte de performance otra medida de control menos utilizada pero común en este tipo de problemas (competiciones Kaggle, por ejemplo), la raíz del logaritmo del error cuadrático medio (o RMSLE).

\[RMSE = \sqrt{\frac{1}{n} \sum_{j = 1}^{n} (log(p_i +1) - log(a_i + 1))^2} \]

Por último, se tendrá el cuenta el tiempo de entrenamiento de los modelos (luego de la puesta a punto de hiperparámetros) a la hora de evaluar su performance.

Predictores

A partir de los datos disponibles, se generaron variables nuevas para predecir la demanda. Un grupo de variables se relaciona a la fecha y la hora en que se hizo el trayecto, por ejemplo: día de la semana, si es laborable o no, hora del trayecto, etc. Otro grupo de predictores está relacionado al clima al momento del trayecto, a nivel diario: temperatura media diaria, temperatura máxima y mínima diaria, velocidad media del viento, presión, etc. Un tercer grupo está relacionado con el clima pero a nivel horario: viento, temperatura, nubosidad y lluvias. Éste último grupo provee información para intervalos de 3 horas, por lo se imputa el mismo valor para cada hora de cada intervalo.

Reducción de dimensionalidad

Como presenta el siguiente mapa de calor, las variables de los grupos relacionados al clima están fuertemente asociadas. Esto puede ser problemático para un modelo predictivo ya que agrega ruido a la predicción y la hace menos estable al predecir nuevos de datos.

Figure 8: Mapa de calor de la matriz de correlaciones. Pasando el ratón por cada celda se presenta su coeficiente. Haciendo zoom-in aparecen el resto de las variables.

Dependiendo del tipo de modelo, se verá más o menos afectado por este efecto (i.e. multicolinearidad). Como medida de limpieza , se elimina una las variables de cada par con más de un 0.75 de asociación (coeficiente Pearson). Son 6, todas climáticas. Además se realiza un Análisis de Componentes Principales (PCA) y la siguiente transformación de los datos para retener los 6 componentes que explican el mayor porcentaje de varianza. Dada la reducción de interpretabilidad que implica éste tipo de transformación, se evalúan también opciones de modelos que no hagan uso de ésta transformación.

Modelos

Se entrenan 4 modelos:

  1. XGBoost con PCA (xgb_PCA): Extreme Gradient Boosting utilzando los 6 componentes del Principal Components Analysis.
  2. XGBoost sin PCA (xgb): El mismo modelo que en el punto 1 pero utilizando todas las variables menos las de correlación mayor a 0.75 (ver el apartado Reducción de dimensionalidad)
  3. Regresión Lineal con PCA (glm_PCA): Modelo de regresión lineal generalizado (GLM), entrenado con los 6 componentes del PCA.
  4. Regresión Lineal con regularización (glm_NET): Modelo de regresión lineal generalizado con mecanismo de regularización ElasticNet.

Para todos los modelos, menos el número 3 se realiza la puesta a punto de hiperparámetros.

Performance

Figure 9: Métricas de performance por modelo. Pasando el ratón por encima, se presentan el valor de cada barra.

El Extreme Gradient Boosting supera con amplias diferencias la performance de los modelos lineales, con o sin penalización y con o sin transformación en componentes principales. Focalizando en el RMSE, el modelo lineal con PCA tiene una desviación standard de la varianza no explicada de ~318 bicicletas por día y hora. Mientras que los modelos de árbol tienen un RMSE de ~130 para el modelo con todos los predictores y ~126 para el reducido (vía PCA).

Una de las ventajas del Extreme Gradient Boosting es la posibilidad de utilizar una gran cantidad de variables predictoras, sin verse muy afectado por la multicolinearidad. Una adecuada selección de hiperparámetros hace que las diferencias de precisión entre los dos algoritmos de boosting sean mínimas.

En términos de tiempo, en cambio, la versión con PCA es muy superior. El costo es una menor interpretabilidad del modelo.

Para saber qué variables fueron las de mayor aporte al modelo podemos ver las 20 primeras del modelo sin PCA en un ranking.

Importancia de las variables en el modelo Extreme Gradient Boosting.

Figure 10: Importancia de las variables en el modelo Extreme Gradient Boosting.

Las variables de temperatura aparecen en los primeros puestos seguidas, como cabe esperar, por dos variables de horas pico (8AM y 7PM).

Conclusión

Los modelos presentados ofrecen algunas opciones con un desempeño razonable en términos de precisión y tiempo. Dependiendo del uso que se le fuese a dar, diferentes camminos serían recomendables. Algunos de ellos se discuten en la siguiente sección.

A considerar

Límites actuales y posibilidades de mejora a futuro.

Base de entrenamiento

En este caso aumentar el tamaño de la base de datos puede ser una vía poco costosa para obtener mayor fiabilidad 3. Series temporales más pronlongadas permitirían observar efectos estacionales. Es probable, por ejemplo, que los usuarios ocasionales aumenten en los meses de primavera mientras se mantienen estables los usuarios anuales. Esto perjudicaría el desempeño del modelo actual pero sería controlable en un período más prolongado de datos.

Profundizar la puesta a punto de hiperparámetros

Durante el entrenamiento de los modelos se evaluaron rangos no muy amplios de hiperparámetros (debido a tiempo y costo computacional). Estos rangos podrían expandirse más para asegurar que no existe una configuración mejor cada modelo.

Datos climáticos hora a hora

En línea con la literatura, las variables climatológicas demostraron ser de gran importancia a la hora de predecir la demanda agregada de bicicletas. Ambas fuentes de datos, de la Agencia Estatal de Meteorología y de Windguru, tienen inconvenientes. Esto se debe a que los archivos de datos no tienen tanta o tan buena información como las consultas en períodos de tiempo más cercanos al actual.

El archivo de la Agencia Estatal de Meteorología sólo provee datos a nivel agregado por día y no cuenta con datos de nubosidad. Windgurú supera estos inconvenientes, pero agrega otros:

En caso fuese a utlizarse este modelo diariamente como parte de un proceso de producción, sería posible extraer datos del pronóstico en tiempo real y seguramente mejorar así la predicción.

Apéndice

  1. El-Assi, W., Mahmoud, M. S., & Habib, K. N. (2017). Effects of built environment and weather on bike sharing demand: a station level analysis of commercial bike sharing in Toronto. Transportation, 44(3), 589-613.↩

  2. García-Palomares, J. C., Gutiérrez, J., & Latorre, M. (2012). Optimizing the location of stations in bike-sharing programs: A GIS approach. Applied Geography, 35(1-2), 235-246.↩

  3. El portal de Datos Abiertos de Madrid ofrece al menos 1 año y medio de datos (desde abril de 2017)↩