Pour ce projet de prédiction de série temporelle, nous avons décidé de nous orienter vers une problématique qui nous intéressait particulièrement, celle des transports. Si le data-mining peut, dans ce domaine, s’appliquer à de nombreuses questions (de la prédiction de flux routier à l’étude du nombre d’usagers d’un réseau ferré), notre choix à été contraint par le peu de jeux de données libre d’accès. En particulier, la RATP, la SNCF, et la société JCDecaux ne mettant pas en ligne de données non agrégées, nous nous sommes tournés vers l’étude d’un jeu de données propre à un système de vélos en libre-service tel que le vélib’ à Paris, mais implanté aux États-Unis. Nous avons considéré plusieurs grandes villes américaines, et avons finalement retenu Washington D.C., qui présentait la plus riche base de données. Le système de vélos en libre-service de la capitale américaine, “Capital Bikeshare” (également appelé “CaBi”), a été mis en place en septembre 2010. Il est aujourd’hui considéré comme le système de vélos en libre-service s’étant le mieux développé. Il propose plus de 4 000 vélos sur plus de 400 stations réparties sur la ville et son agglomération (comtés d’Arlington,de Fairfax, de Montgomery et ville d’Alexandria).
Les données disponibles en ligne consistent en des tableaux regroupant par trimestre l’ensemble des trajets effectués par les vélos de la société CapitalBikeShare. Chaque ligne, correspondant à un trajet, contient la durée du trajet, la date de départ (jour et heure), la date d’arrivée (jour et heure), la station de départ, la station d’arrivée, le numéro du vélo et le type d’utilisateur (abonné ou occasionnel). La première étape de traitement des données a donc été de rassembler ces différents fichiers dans un même data frame, en convertissant les dates et les durées (parfois notées avec des conventions différentes suivant les trimestres), regroupant des colonnes dont les noms changeaient et en supprimant les colonnes redondantes qui apparaissaient pour certains trimestres seulement.
Une fois ce nettoyage effectué, pour prendre en main le jeu de données, nous avons cherché à le visualiser. Nous nous sommes tout d’abord intéressés à l’évolution du nombre de vélos empruntés au cours de la semaine.
L’influence du jour de la semaine sur le nombre de trajets semble faible : en particulier on ne détecte par de baisse du nombre de vélos empruntés en semaine ou le week-end.
Nous nous sommes ensuite intéressés aux heures d’emprunt des vélos.
On observe deux pics correspondant aux heures d’arrivée au et de départ du lieu de travail. On émet l’hypothèse que le service proposé par CapitalBikeShare est principalement utilisé pour faire les trajets entre le domicile et le lieu de travail.
Les données liées aux vélos étant intrinsèquement spatiales, nous souhaitions visualiser les trajets empruntés; cependant le grand nombre de trajets et de stations rendait le résultat difficilement interprétable. Nous avons donc changé de stratégie. Nous avons calculé le nombre moyen de vélos empruntés par heure et par station, puis avons, pour chaque heure, classifié ces stations selon le nombre de vélos qui y étaient déposés (plus précisément, dans les plans si dessous, une station est étiquetée -1 si le nombre de vélos déposés ou pris est inférieur de plus d’un écart type du nombre moyen de vélos déposés ou pris par heure).
En affichant cette classification sur une carte pour les différentes heures, nous avons mis en évidence l’importance des trajets de la banlieue vers le centre ville vers 8 et 9h, et en sens inverse vers 17, 18 et 19h, ce qui confirme notre hypothèse. En général, peu de vélos sont empruntés la nuit.
Ce faisant, nous avons mis en évidence le fait que CapitalBikeShare se chargeait de déplacer les vélos des stations encombrées vers les stations plus vides. Ci-dessous nous affichons la carte des stations, chaque station étant colorée suivant le nombre de trajets partant de la station moins le nombre de trajets y arrivant. On observe notamment que le nombre de vélos redéposés par CaptialBikeShare en centre ville est plus faible que le nombre dé vélos partant du centre ville (les utilisateurs font plus de trajet vers le centre ville que vers la banlieue, d’altitude plus haute).
Les données brutes sont trop riches pour que nous puissions les exploiter telles quelles. Dans un premier temps, nous avons décidé des nous intéresser au problème de la prédiction du nombre de trajet par jour. Ce problème, simple, nous semblait plus abordable en première approche que celui de la prédiction du nombre de trajets par heure (en particulier, nous n’avions pas accès aux données météo heure par heure, ce qui compliquait considérablement la phase de pré-processing des données), ou station par station. Nous nous intéressons aux données postérieures au 1er janvier 2011, le dernier trimestre de 2010 correspondant à la mise en place du service.
Pour plus de lisibilité, nous avons aussi étudié la série lissée en effectuant une moyenne mobile de largeur un mois.
On observe que le nombre moyen par an de trajets augmente au cours du temps. Ceci s’explique par le développement du service de location qui a multiplié ses stations et le nombre de vélos en circulation. Pour tenir compte de cette augmentation, nous avons souhaité ajouter comme variable explicative le nombre de vélos en circulation et le nombre de stations. Comme ces données n’étaient pas disponibles, nous avons dû les estimer. Nous l’avons fait de la manière suivante. Chaque vélo étant repéré par un identifiant unique, il suffisait de repérer la première date d’apparition de chaque vélo, puis le nombre de vélos apparus à chaque date. Nous avons corrigés ce chiffre par celui des vélos retirés de la circulation (obtenu de la même façon), puis nous avons tenté d’effacer le bruit en modélisant le nombre de vélos réellement en service comme un fonction continue par morceaux, à l’aide du package changepoint. Nous avons procédé de même pour estimer le nombre de stations.
Pour prédire le trafic journalier, nous avons également utilisé comme variables explicatives des données météorologiques (fournies par le National Climate Data Center). Pour des informations plus précises, on peut trouver la documentation relative à ces données météorologiques sur le même site ou directement ici. Nous avons utilisées les données récoltées par la station de l’aéroport Reagan, en agrégeant certains facteurs qui apparaissaient très rarement (tels que différents niveaux de brume ou de pluie).
D’autres données auraient pu être intéressantes mais sont inacessibles. Par exemple, nous aurions voulu pouvoir utiliser des données concernant le traffic automobile de Washington (kms de bouchons, ralentissements en terme de durée, accidents .), le réseau de transport en commun (lignes fermées/ouvertes, en maintenance, grèves .) ou les différents évènemens dans Washington et sa banlieu, car ces dernières pourraient contenir de l’information utile à la résolution du problème.
Une fois les données brutes transformées et agrégées aux données météorologiques, nous avons pu appliquer les différentes méthodes vues en cours pour prédire le traffic journalier.
Nous avons tout d’abord tenté de prédire le traffic journalier à l’aide d’une simple régression linéraire. Les paramètres explicatifs sont les suivant :
Le modèle entraîné sur les données de 2011 à 2016 prenant en compte ces paramètres présentait de nombreux défauts. Tous d’abord, un certain nombre de ces variables étaient non significatives. Ensuite, il apparaissait que les résidus n’étaient pas identiquement distribués suivant les valeurs des différents paramètres : par exemple, on observait toujours une tendance dans les résidus par rapport à la température, et une oscillation dans les résidus affichés en fonction du temps. Pour corriger ces défauts, nous avons testé différents modèles comportant des termes d’ordre 2 et 3 en la température maximale et en le nombre de vélos, et avons ajouté une base de Fourier de 10 sinus et 10 cosinus, de période 52.2 semaines. Nous avons sélectionné pour chaque nombre de variables le meilleur sous modèle grâce au package leaps, puis nous avons sélectionné le meilleur sous modèle parmi cette liste de modèle par validation croisée. Cependant le modèle obtenu n’était pas satisfaisant pour plusieurs raisons qui nous semblaient liées au même problème : par exemple, les résidus n’étaient pas homoscédastiques. L’explication tient au fait que le nombre de vélos en service augmente au cours du temps, et donc que logiquement la variance du nombre de trajets augmente également. Pour prendre en compte ce phénomène, il faudrait donc faire un modèle exprimant une dépendance linéaire en le nombre de trajet renormalisé par le nombre de vélos en service et le reste des variables explicatives. Plus exactement, en observant l’augmentation du nombre de trajets moyens par rapport à celle du nombre de vélos en service, nous avons choisit d’utiliser un modèle exprimant une dépendance linéaire du nombre de trajets divisé par le nombre de vélos en service en le reste des paramètres. Nous avons, là aussi, augmenté les variables explicatives de termes quadratiques et cubiques en la température maximale et par une base de Fourier de 10 sinus et 10 cosinus. Nous avons ensuite sélectionné le meilleur sous-modèle à l’aide d’une méthode Gauss-Lasso, en choisissant la pénalisation proposée par E. Candès : \(2*\sigma*\sqrt((2 \log (p)/n)\) où \(n\) est le nombre d’observations, \(p\) le nombre de variables explicatives, et \(\sigma\) une estimation de la variance, calculée sur le modèle complet. Le modèle final prend en compte les variables suivantes : Wday, PRCP, SNWD, TMAX, WSF2, FOG, RAIN, cos1, cos2, sin2, sin4. Pour obtenir une régression minimisant l’erreur sur le nombre de trajets non renormalisé par le nombre de vélos en service, nous avons pondéré notre régression linéaire par le nombre de vélos disponibles. Lors d’une précédente tentative, nous avions voulu prendre en compte le fait que le nombre de trajets était nécessairement positif, et construire un modèle reliant linéairement les paramètres et le logarithme du nombre de trajets. Ce modèle n’était pas concluent (les erreurs se propagaient considérablement lors du passage à l’exponentielle pour obtenir le nombre de trajets prédit). Cette fois, nous avons simplement choisi de proposer comme prédiction la partie positive du nombre de trajets préduit par le modèle linéaire.
##
## Call:
## lm(formula = expr.gauss.lasso, data = TrainData, weights = TrainData$velos.dispo)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -5449.3 -840.1 25.2 791.0 6450.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 92.32921 2.92269 31.590 < 2e-16 ***
## Wday 1.87364 0.27010 6.937 5.26e-12 ***
## PRCP -1.64182 0.07504 -21.879 < 2e-16 ***
## SNWD -0.19896 0.02440 -8.156 5.80e-16 ***
## TMAX 3.11769 0.09380 33.239 < 2e-16 ***
## WSF2 -1.29588 0.19996 -6.481 1.13e-10 ***
## FOG -9.21931 1.27980 -7.204 8.03e-13 ***
## RAIN -18.48693 2.01311 -9.183 < 2e-16 ***
## cos1 -15.16285 1.30876 -11.586 < 2e-16 ***
## cos2 -12.71404 0.76739 -16.568 < 2e-16 ***
## sin2 -9.37471 0.77417 -12.109 < 2e-16 ***
## sin4 4.72980 0.76592 6.175 7.85e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1285 on 2179 degrees of freedom
## Multiple R-squared: 0.7782, Adjusted R-squared: 0.7771
## F-statistic: 694.9 on 11 and 2179 DF, p-value: < 2.2e-16
On remarque que les variables explicatives retenues ont bien l’effet qu’on attendrait intuitivement : la pluie, la neige, le vent, la brume entraînent une diminution du nombre de trajets, tandis que plus il fait chaud, plus le nombre de trajet est important.
Nous testons la qualité du modèle entraîné sur les données du 01/01/2011 au 31/12/2016 en calculant la RMSE des prédictions sur la période du 01/01/2017 au 31/12/2017. La RMSE obtenue est 1725 (pour un nombre moyen de trajets par jour sur la période de 10 298 trajets, et un écart type de 3728), tandis que la RMSE obtenue pour le premier modèle linéaire “naïf” était de 2511. On observe de plus que les résidus du modèles sont bien hétéroscédastiques, et ne présentent pas de tendance évidente quand on les représente en fonction des variables explicatives. On constate de plus que le modèle n’est pas très performant pour expliquer le nombre de trajets effectués la première année (il le surestime). On suppose que ce phénomène est dû à un régime différent lié à la nouveauté du service proposé (les clients ne connaissent pas très bien ce service, ne savent pas où sont situées les stations, etc).
Dans la suite du projet, nous avons voulu comparer ce modèle à d’autres modèles. Le premier modèle alternatif que nous avons testé est un modèle CART. Nous avons pour cela utilisé le package rpart, avons estimé l’arbre maximal sur les données du 01/01/2011 au 31/12/2016 et l’avons ensuite élagué.
## [1] "RMSE"
## [1] 1928.77
On remarque que le noeud à la racine de l’arbre divise les données en deux selon la température (inférieure ou non à 18 degrés Celsius), tandis que la variable explicative sélectionnée ensuite pour les deux noeuds fils est le nombre de vélos.L’arbre obtenu comporte 1825 noeuds. Nous avons ensuite évalué ce modèle en calculant la RMSE sur l’année 2017. La RMSE obtenue est de 1907. Celle-ci est plus élevée que la RMSE du meilleur modèle linéraire : on peut supposer que ceci est dû au fait que le service a continué à se développer en 2017, le nombre de vélos en service a augmenté, et que notre modèle CART peine à généraliser ses prédictions. En particulier, on observe que le nombre de trajet prédit est souvent inférieur aux observations, ce qui est particulièrement flagrant en été (lorsque le nombre de trajets maximal est atteint).
On essaye d’améliorer la prédiction à l’aide d’un modèle de forêt aléatoire. On utilise pour cela le package randomForest. On optimise successivement les paramètres ntree (nombre d’arbres), mtry (nombres de variables tirées au hasard à chaque noeud), sampsize (nombre d’observations pour le bagging), et nodesize (taille minimale d’une feuille), qui sont respectivement fixés à 100, 4, 400 et 1. Pour ce faire, on a divisé les données en un échantillon d’entraînement (du 01/01/2011 au 31/12/2015) et de validation (du 01/01/2016 au 31/12/2016). Une fois les paramètres choisis, on entraîne à nouveau l’arbre sur l’union des deux échantillons.
La RMSE obtenue pour les données de l’année 2017 est de 2043. Là encore, on observe que la forêt aléatoire peine à généraliser les résultats des années précédentes pour prendre en compte la croissance du service : le modèle sous-estime toujours le nombre de vélos empruntés, ce qui est particulièrement visible en été.
Parmi les modèles mis en oeuvre, le modèle le plus performant est le modèle par régression linéaire. Néanmoins, ce modèle est celui qui nécessite la plus longue optimisation, alors que le modèle CART, qui performe légèrement moins bien, a l’avantage de se calibrer très rapidement. De plus le modèle le plus performant est celui dans lequel nous avons incorporé notre conviction selon laquelle le nombre de trajet devait être renormalisé par la racine carrée du nombre de vélos en service pour être linéairement dépendant des autres paramètres. En ce sens, les modèles d’arbre CART et de forêt aléatoire ont l’avantage d’être agnostiques.
Une des idées que nous avons eu pour répondre à un autre problème sur ce même jeu de données qui ne s’éloigne pas trop du problème de départ a été d’essayer de prédire le pic de consommation (par tranche de dix minutes ~ durée médiane d’un trajet) sur chaque journée. En pratique, cela permettrait d’anticiper différents problèmes rencontrés par les utilisateurs tels que des stations pleines (impossibilité de reposer le vélo emprunté) ou vides (impossibilité d’emprunter un vélo) et qui le restent pour une certaine durée. Par manque de temps, nous n’avons pas pue exploiter au mieux tous les modèles. Parmi ceux que nous avons essayé, nous avons retenu un modèle par GAM, que nous avons ensuite amélioré en appliquant un modèle ARIMA sur les résidus.
La partie d’optimisation des paramètres et de sélection de modèle a étée réalisée de manière un peu hâtive mais le résultat final n’est pas si mauvais, bien qu’il puisse probablement être amélioré. Nous avons effectué une sélection de variables en s’inspirant du modèle linèaire car la série du traffic journalier et de la hauteur du pic journalier semblent assez corrélées (ci-dessous en noir la série des pics et en bleu la série des trafics journalier multipliée par un coefficient (0.25)).
Nous avons finalement retenu le modèle gam suivant: “PEAK ~ s(Wday) + te(X,Yday) + s(PRCP) + s(BikesDB) + s(TMAX) + s(SNWD)” Nous avons simplement choisi comme base de splines ‘cr’ pour les variables “X”,“PRCP”,“BikesDB”,“TMAX” et “SNWD” et ‘cc’ pour les autres variables qui sont cycliques (“Wday”,“Yday”).## [1] "Training RMSE"
## [1] 28.92
## [1] "Testing RMSE"
## [1] 51.35
Le modèle est plutôt satisfaisant mais en regardant l’autocorrélogramme et l’autocorrélogramme partiel des résidus, on se dit qu’on peut sûrement faire des améliorations.
On peut notamment voir sur ces graphiques que la périodicité hebdomadaire est toujours présente dans les résidus. Pour tenter d’améliorer notre modèle, on entraîne un modèle arima sur les résidus (sélection de ce modèle par la fonction auto.arima).
## [1] "Training RMSE"
## [1] 27.47
## [1] "Testing RMSE"
## [1] 48.77
On observe que le problème est à peu près équivalent à celui de la prédiction du trafic journalier, que ce soit en terme de difficulté ou en terme de limite de performance des méthodes classiques appliquées dans le cours. De plus, la prédiction de la hauteur exacte du pic journalier est peut-être moins intéressant que la prédiction d’un dépassement de seuil pour les pics (qui correspondrait grossièrment à une situation de forte disparité dans la répartition spatiale des vélos parmi les stations), accompagnée d’une prédiction des heures auquelles le seuil est dépassé. Néanmoins la résolution du problème pratique, anticipation des pénuries et surabondances de vélos à chaque station, semblent trop complexe pour n’être résolu qu’avec les méthodes de prévision de séries temporelles qu’on utilise ici. Il semble se prêter beaucoup plus à une résolution via des méthodes d’apprentissage en ligne avec des données également plus complexes (nombres vélos empuntés/déposés à chaque station par tranche de dix minutes). Cela nécessiterait sûrement de prendre en compte la dimension spatiale des données et pas seulement l’aspect temporel. S’intéresser à ce problème nous a permis de saisir un peu plus la richesse mais aussi la complexité des données brutes. Différentes pistes de réflexion s’ouvrent à celui qui veut améliorer notre approche. On peourrait par exmple reprendre le même principe que celui que nous avons adopté en diminuant le pas de discrétisation (ici 24 heures). Le cours s’intéresse davantage à des séries temporelles, mais comme nous venons de le dire, introduire de la spatialité dans le jeu de données, via un clustering par exemple, pourrait se révéler pertinent. Bien que très intéressante, cette partie à amener beaucoup plus de questions que de réponses, ouvrant notamment sur des problèmes dont la résolution semble nécessiter beaucoup de temps et de connaissances pratiques, que nous n’avons pas.
Les méthodes de prédictions mises en oeuvre dans ce projet nous ont permis d’obtenir des modèles satisfaisant. Cependant on peut regretter que les séries temporelles à prédire n’aient pas présenté de réels challenges (une fois la phase d’agrégation des données passée, phase particulièrement pénible en raison des différentes conventions employées). Il aurait sans doute été plus intéressant de prédire le nombre de vélos empruntés dans la prochaine heure pour les différentes stations, cependant ceci posait plusieurs problèmes. Tout d’abbord, nous n’avions pas accès à des données météos heures par heures. Ensuite, nous n’avions pas directement accès au nombre des stations à chaque date, ni à leur capacité. De plus, le nom des stations changeant au cours du temps, il aurait sans doute été assez pénible de reconstruire une série temporelle par station. Enfin, le jeu de données correspondant aux trajets par heure et par station aurait été plus lourd et les calculs plus longs, nous avons donc préféré appliquer les techniques vues en cours à un jeu de données “jouet” sur lesquels la durée des calculs n’était pas un facteur limitant. Cependant, Il aurait sans doute été intéressant d’appliquer un modèle linéaire à effet mixtes à ce problème, et de comparer l’influence des différentes variables explicatives. Une autre question que nous nous sommes posée était liée à la classification des stations entre “lieu de travail” et “lieu de résidence”. On a vu en effet que l’histogramme des vélos empruntés varie d’une station à l’autre; il aurait été intéressant de parvenir à classifier ces stations suivant leur mode d’utilisation. Dans cet esprit, nous avons fait une ACP sur la matrice ayant comme ligne les différentes stations, et comme colonne les nombres de vélos empruntés et déposés par heure. Il en est ressorti que la variance du dataset est bien expliquée par les trois premières directions de l’ACP, qui comptent pour 80% de la variance. Cependant ceci ne nous permettait pas de classifier facilement les stations. Pour procéder de manière plus agnostique, nous pourrions utiliser des méthodes de classification non supervisée comme des méthode de k plus proches voisins, avec une distance mesurant la similarité entre deux histogrammes décrivant la station (par exemple la distance de Wasserstein).