Prédire la circulation routière à Paris

Projet ML, M2 Statistiques et Machine Learning

1 Introduction : la circulation à Paris

La ville de Paris rend publiquement accès à de nombreux jeux de données sur la circulation routière de la ville. La découverte du site web Paris Data (“Paris Data” visité le 26.12.2021) nous a motivé à utiliser cette ressource dans notre projet de prédiction. En voyant les données, nous avons développé l’idée de relier les différents capteurs dans une structure de graphe. D’une part, cela motive des études du réseau entier : peut-on construire un algorithme global, qui apprend les interactions dans le réseau entier ? D’autre part, l’étude de rues individuelles mène à des questions sur la corrélation locale de la circulation. Un grand boulevard se comporte-t-il comme ses avenues voisines ?

Dans une première partie, nous introduisons le jeu de données étudié dans notre projet. Nous expliquons notre démarche de création du jeu de données en détaillant les variables récoltées, leurs sources et la complétion de leurs valeurs manquantes. Une brève analyse descriptive sera également réalisée. Dans une seconde partie, nous réalisons la tâche de prédiction de notre projet. Après avoir introduit les problématiques, la métrique utilisée et les détails de notre stratégie de prédiction, nous rapportons les différentes méthodes utilisées et leurs résultats. Les modèles classiques sont passés en revue puis agrégés par une agrégation d’experts.

Matériel accompagnant ce document

Ce document est le rapport de notre projet de Machine Learning, il contient peu de code et encore moins de données. Pour accéder à ces ressources, nous recommandons de cloner notre repository sur github : https://github.com/jakobmaier-eu/paris-traffic-prediction.git . Dans le rapport, nous faisons parfois des références à ce code.

Par ailleurs, nous avons préparé un package qui contient les fonctions les plus utilisées dans nos calculs ainsi que les données que nous avons utilisées pour les prédictions. Il peut être installé à partir du fichier .tar.gz qui peut être téléchargé suivant ce lien One Drive : https://1drv.ms/u/s!AorCvx6dgh2IhIo6VoTCkoM5A-Cw5A?e=sdiFvV .

2 Préparer le jeu de données

2.1 Les données et leurs attribus

Sur le site Paris Data, on trouve énormément de jeux de données issus des activités de la ville, notamment sur le comptage routier (“Données Du Comptage Routier” visité le 26.12.2021). Ces données sont collectées en continu grâce à des boucles électromagnétiques implantées dans la chaussée sur plus de 3000 tronçons de voies. Ces boucles et bornes de comptage mesurent principalement deux choses: le nombre de voitures qui passent et l’occupation du tronçon de route. Disponibles de 2014 jusqu’à aujourd’hui, ces deux données sont mesurées au pas horaire, résultant en environ \[ 8 \text{ ans } \times 365 \text{ jours } \times 24 \text{ heures } \times 3000 \text{ capteurs } = 2.1\times10^8 \text{ observations dans le data set.} \]

Chacune de ces lignes contient les informations suivantes, qui serons nos variables explicatives \(X_i\). Les noms des variables que nous utilisons dans notre analyse sont écrits en gras tout au long de cet exposé.

• Un timestamp t_1h indiquant date et heure de la mesure.
• Le nombre de voitures ayant passé le point de comptage pendant une heure, noté q et renommé par nbCar.
• Le taux d’occupation, correspondant au temps de présence de véhicules sur la boucle en pourcentage d’une heure, noté k et renommé par rateCar.
• La variable etat_barre indique si la route où se trouve la borne de comptage est “barée,” etat_barre=2, ou “ouverte,” 1. Nous la renommons state. Malheureusement, cette variable a souvent la valeur “inconnu,” 0, ou “invalide,” 4. La différence entre 3 et 4 est un mystère et par conséquence nous nous servons peu de cette variable.
• L’identifiant iu_ac de chaque borne. Chaque ligne contient de informations sur les compteurs adjacents.
• Le libellé (libelle) du tronçon de rue où est placé le capteur, par exemple “Champs Elysées partie B.”

Les données historiques de circulation sont enregistrées sur le site web Paris Data sous le format .txt. Chaque fichier couvre toutes les bornes pour une semaine. Pour commencer, il faut donc télécharger les fichiers de toutes les semaines et les convertir en dataframes. Ce faisant, nous effectuons également un tri par identifiant de borne iu_ac.

2.2 Simplification et aggrégation

Au vue de la quantité énorme de 200 millions d’observations, nous sommes obligés de réduire la taille du dataset. La première démarche consiste à garder les variables intéressantes: t_1h, nbCar, rateCar et iu_ac.

Ayant réduit la dimension de chaque observation, la deuxième étape est la réduction du nombre de capteurs de comptage. Pour ce faire, nous utilisons un modèle simplifié des rues de Paris. Etant donné que les capteurs couvrent les grandes axes routiers, dont le périphérique, notre objectif était d’agréger les données selon ces grands axes. Afin d’identifier où passent le plus de voitures, nous calculons les moyennes de nbCar pour chaque libellé, c’est-à-dire les moyennes à travers le temps et pour toutes les bornes associées à chaque libellé. En prenant les 200 libellés dont le nombre de voitures est le plus grand (hors périphérique), nous obtenons le graphique 2.1 qui représente très bien les axes routiers auxquels on s’attendait.
Ensuite, nous ignorons les libellés et regardons exclusivement ce graphique. De manière arbitraire, nous en déduisons le graphe de la figure 2.2, un schéma très simplifié de la circulation à Paris. Grâce à notre démarche d’utiliser les rues les plus fréquentées, cette abstraction nous permet d’avoir un modèle représentatif de la circulation parisienne.

Figure 2.1: Représentation en rouge des bornes d’observation associés aux 200 libellés les plus fréquentés (hors périphérique)

Figure 2.2: Graphe simplifié. Les points noirs indiquent les intersections entre nos arêtes.

Chaque arête dans le graphe 2.2 regroupe plusieurs capteurs dont nous agrégeons les données : après une collecte minutieuse et laborieuse des identifiants de toutes les bornes sur chacune des arêtes, nous associons pour chaque heure la moyenne de nbCar et rateCar des bornes de l’arête dont ils font partie. Cela a trois avantages:

1. Prendre la moyenne évite d’avoir des données manquantes, issues entre autres de capteurs dysfonctionnels sur une certaine période. Ce n’est pas toujours le cas, comme on le verra dans la section suivante.
2. S’il y a des capteurs où les données fluctuent beaucoup en raison de mauvais positionnement ou des perturbations indésirées, la moyenne réduit cette variance.
3. Le problème est réduit à 69 arêtes représentées par 69 dataframes.

Cette réduction des données résulte en la construction des jeux de données que l’on va manipuler et en la concrétisation du problème que l’on souhaite résoudre. Elle permet également d’obtenir des objets que nos ordinateurs sont capables de stocker efficacement. Il reste pourtant un problème: la présence de données manquantes.

2.3 Complétion de valeurs manquantes

En examinant des figures basiques des taux d’occupation, nous observons régulièrement une image comme dans le graphique 2.3 Malgré l’agrégation le long des arêtes, il reste un nombre considérable de valeurs manquantes, NA, dans le dataset. Afin d’exécuter des tâches de prédiction, il nous faut cependant des données complètes. Une partie conséquente de notre travail a donc été de compléter les valeurs manquantes.

Figure 2.3: Taux d’occupation pour une journée ouvrière avec une valeur manquante

Figure 2.4: Visualisation de valeurs manquantes pour 30 arêtes (gris = NA)

En regardant le graphique 2.4, on considère qu’il y a deux types de “trous” dans le data set:

1. Des valeurs manquantes ponctuelles comme on a pu observer dans la figure 2.3 ou celles marquées en bleu dans le graphique ci-dessus. On peut deviner qu’il s’agit de simples erreurs dans la collection des données ou des travaux d maintenance par exemple.
2. Des trous plus grands, couvrant parfois plusieurs mois, comme ceux marqués en rouge dans le graphique 2.4. Il peut s’agir d’une absence étendue de plusieurs bornes, une clôture de la route en raison d’une construction, etc. Savoir si une route est effectivement inaccessible pendant de longues période serait une tâche de recherche en soi : comme mentionné plus tôt, la variable state ne contient presque aucune information. Nous remplissons alors les grands trous sans respecter d’éventuelles fermetures, ajoutant une autre simplification.

Le fait que les arêtes ont de différentes structures de valeurs manquantes peut aussi être observé en regardant la répartition du pourcentage de NA’s à travers les dataframes dans la figure 2.5.

Figure 2.5: Distribution des valeurs manquantes

Nous observons que les données de taux d’occupation sont généralement plus complètes : ceci est dû au fait qu’il y a plus de capteurs de ce type. Par ailleurs, nous avons vérifié qu’une valeur manquante de nbCar implique toujours un NA dans rateCar.

Pour remplir les trous, la première étape consiste à ajouter des variables explicatives supplémentaires. En fonction du type de trou, nous verrons qu’ils auront un impact important sur la complétion.

2.3.1 Variables supplémentaires : Proximité temporelle et locale

Pour l’instant, nous disposons des variables nbCar et rateCar pour chaque arête et chaque heure entre 2014 et 2020. Nous exploitons deux propriétés structurelles des données pour obtenir de l’information supplémentaire:

D’abord, il s’agit de séries temporelles, donc nous rajoutons les valeurs historiques de nbCar et rateCar décalées d’une heure, d’un jour et d’une semaine (nbCarLaggedHour, nbCarLaggedDay, et nbCarLaggedWeek) comme variables explicatives. L’espoir étant que la circulation reste identique à travers ces cycles. Avec l’objectif spécifique de compléter les données manquantes, il est aussi pertinent d’ajouter ces valeurs du futur (rateCarFuturHour etc.), étant donné que le modèle d’imputation a pour but d’interpoler au lieu d’extrapoler. On fait l’hypothèse que toutes ces variables décalées seront efficaces pour remplir les trous locaux : une heure ou une journée manquante devrait facilement être inférée en utilisant des données qui sont proches dans le temps.

Pour les valeurs manquantes au milieu des grands trous, nous n’avons pas accès aux heures et jours d’avant car elles sont aussi manquantes. Cependant, grâce à l’interconnexion de nos données à travers le graphe présenté plus haut, les mesures de circulation d’une arête peuvent bien être expliquées en fonction de celle de ces voisines. Dans un effort manuel, nous avons pour chaque arête \(A\) collectionné des voisins \(V_i\), dans le sens que les voitures dans \(A\) passent ensuite par l’un des \(V_i\) ou bien l’inverse. Nous faisons ici encore des choix arbitraires de qui est voisin de qui pour deux raisons. Premièrement, inférer statistiquement quelle arête influence quelle autre nécessiterait des données déjà complètes pour une régression par exemple. Deuxièmement, utiliser toutes les autres 68 arêtes au lieu de juste 2 à 6 voisins augmenterait trop le temps d’exécution de la complétion.
Donnons un exemple illustrant le choix de voisins ainsi que la dénomination des nouvelles variables. Parmi les voitures qui entrent dans Paris par l’arête “pont amont - pont austerlitz,” il y a une partie considérable qui continue tout droit sur l’arête “pont austerlitz - chatelet.” Par conséquent, nous ajoutons la variable rateCar_pontausterlitzTOchatelet au dataframe de “pont amont - pont austerlitz.”

Jusqu’ici, nous n’avons utilisé que des informations déjà présentes dans le jeu de données. Dans la suite, nous ajoutons des variables extérieures qui pourraient également expliquer la circulation routière.

2.3.2 Variables supplémentaires : Facteurs externes

Variables temporelles

Le traffic routier étant relié à l’activité humaine, nous avons ajouté de nombreuses variables temporelles (notamment grâce à la librairie lubridate).

• year, month, day, hour résultent d’une simple décomposition des timestamps t_1h,
• time est un simple compteur de journées, commençant à \(1\) le 1er janvier 2014,
• toy (time of year), il s’agit d’un numéro entre 0 et 1 indiquant la position de l’observation dans l’année en cours,
• weekdays le numéro du jour de la semaine et weekendsIndicator l’indicatrice si le jour est un jour du weekend,
• winterHolidaysIndicator et summerHolidaysIndicator les indicatrices des vacances d’hiver et d’été définies à partir de (“Les Archives Du Calendrier Scolaire” visité le 26.12.2021),
• bankHolidaysIndicator l’indicatrice des jours fériés définie à partir de (“Jours Fériés” visité le 26.12.2021).

Index de la situation sanitaire en rapport avec le Covid-19

Nous avons également ajouté, covidIndex, un index allant de 0 à 100 représentant les restrictions du gouvernement sur la situation sanitaire en rapport avec le Covid-19. Il est calculé à partir de nombreux indicateurs et est fourni par l’université d’Oxford (“Covid Index” visité le 26.12.2021). Cette variable permettra éventuellement une étude de l’année 2020 qui voit une circulation perturbé. Pour cela, il faudra utiliser des modèles qui s’adaptent très vite à l’influence de nouvelles variables.

Météo

A partir de données de l’Organisation Météorologique Mondiale (“Observations Météorologiques Historiques En France” visité le 26.12.2021), nous avons extrait 2 variables météorologiques : temperature la température en Kelvin et precipitation les précipitations dans les 3 dernières heures en mm. Ces données ont été mesurées à Athis-Mons en Essonne.

On dispose d’un relevé tous les 3 heures environ donc on procède à une interpolation pour compléter les données. Etant donné que les mesures sont uniformément réparties, on utilise une interpolation linéaire basique à l’aide de la fonction na.approx de la librairie zoo.

2.4 Complétion via Random Forest

Avec un dataframe enrichi à notre disposition pour chacune des arêtes du graphe de circulation, nous pouvons passer au remplissage des valeurs manquantes. Nous sommes pourtant restreint dans notre recherche d’un algorithme d’imputation car plusieurs variables explicatives contiennent également des NA. D’une part, nous voulons remplir deux variables en même temps (nbCar et rateCar). D’autre part, les voisins n’ont pas moins de valeurs manquantes, laissant de grands trous dans le dataframe.
Heureusement, il y a un modèle qui peut être entraîné en dépit de valeurs manquantes parmi les variables : les forêts aléatoires constituées d’arbres CART. Une règle de décision dans un tel arbre peut être ignorée si la valeur de la variable nécessaire n’est pas renseignée. Le package miceRanger profite de ce fait en utilisant un algorithme dit d’imputation multiple : en commençant par la variable \(V_1\), une forêt aléatoire est entraînée sur les lignes où \(V_1\) n’est pas manquant. Puis, les autres lignes sont “prédites” par la forêt. Avec le nouveau dataset moins vide, le processus est reproduit et ainsi de suite (Van Buuren 2018). Nous nous arêtons pourtant à la deuxième itération vu qu’il n’y a pas d’intérêt à compléter les données des voisins.

Comme nous répétons le processus de complétion 69 fois, le temps d’exécution est de quelques heures sur nos ordinateurs. Nous n’avons donc pas le luxe d’optimiser des hyperparamètres comme nous ferons plus tard pour la tâche de prédiction. Le choix de 100 arbres par forêt et 7 variables considérées à chaque split semble un bon compromis. (voir le fichier preprocessing/missing_data.R sur git pour le code exact).

Au délà des capacités de complétion, les forêts aléatoires permettent de calculer des scores d’importance pour chaque variable explicatives. Ces scores peuvent aider à identifier lesquelles des variables ont le plus contribué à l’estimation de la cible (le score est calculé comme la réduction de variance à chaque split). Pour une arête du graphe, jussieu - saint michel, nous illustrons les variables qui ont le plus contribué à l’imputation :

Figure 2.6: Variables avec le score d’importance le plus haut pour rateCar de l’arête ‘’jussieu - saint michel’’

Figure 2.7: Variables avec le score d’importance le plus haut pour nbCar de l’arête ‘’jussieu - saint michel’’

Sur l’arête dont traitent les deux plots, on observe que rateCar et nbCar exhibent deux comportement assez différents. Pour le taux d’occupation, ce sont presque exclusivement les valeurs temporellement décalées qui contribuent à l’imputaion. Comme il y seulement 0.8% de valeurs à compléter, ceci confirme notre hypothèse que les petits trous d’une ou plusieurs heures sont très bien complétés par les valeurs des heures juste avant ou après. Pour le nombre de voitures, où la proportion de NA était plus importante, les deux variables ayant le plus réduit les variances aux splits sont des variables nbCar issues d’arêtes voisines. Les trous importants en nbCar entre Jussieu et Saint Michel semblent être le mieux expliqué par ce qui se passe autour. Cela confirme notre choix de variables pour combler les 2 types de trous. Un phénomène similaire peut s’observer dans les autres arêtes.

Après l’imputation arête par arête, nous pouvons enfin compléter toutes les données grâce au fait que les variables reliées aux voisins ont été imputées dans les dataframes qui correspondent à ces voisins. Il est important de noter ici que nous allons dans la suite couper le jeu de données en deux (train - test) et que cette imputation sera effectuée individuellement sur chacune des deux parties. Faire ceci deux fois séparément est nécessaire pour préserver l’indépendance et permet ensuite d’estimer les erreurs de nos modèles.
Dans la suite, nous supposons que les données sont complètes. Il faut pourtant remarquer que ces forêts aléatoires font indirectement partie des modèles de prédiction que nous allons utiliser. De plus, nous gardons à l’esprit que nous avons introduit un biais dans nos données d’apprentissage comme de test.

2.5 Analyse descriptive des données

Maintenant que nous avons collecté et complété toutes les variables que l’on souhaite utiliser, nous pouvons réaliser une brève analyse descriptive de celles-ci. Premièrement, la variable rateCar étant une série temporelle, on étudie ses saisonnalités et tendance. On peut voir qu’elle possède des saisonnalités journalière, hebdomadaire et annuelle. Cela est typique des activités humaines. Représentons une semaine de mesure quelconque d’une certaine arête:

Figure 2.8: Taux d’occupation de l’arête ‘’saint-michel châtelet’’ la semaine du lundi 6 janvier 2014

On constate bien une saisonnalité journalière. De plus, on voit un schéma similaire pour les 5 jours de la semaine et un schéma particulier les jours de week-end. C’est une motivation de l’ajout de la variable weekendsIndicator. On peut interpréter ces mesures de l’activité routière. En effet, il semble logique qu’il n’y ait peu de bouchons la nuit et qu’il y en ait aux heures de pointe. Ces dernières sont caractérisées par des ‘’doubles pics’’ le matin et en fin d’après-midi. Enfin, on voit que le traffic est plus dense les jours de week end que les jours de semaine. Ces saisonnalités et interprétations se retrouvent sur toutes les arêtes, bien que sur certaines ce sont les jours de semaine où la circulation est plus dense.

On étudie maintenant la corrélation de nos variables. On procède en 2 étapes: premièrement, on l’évalue sur les variables décalées et quelques variables descriptives, ensuite on l’évalue sur les voisins.

Figure 2.9: Représentation de la corrélation des variables décalées et de variables descriptives

On voit directement que les variables décalées sont toutes positivement corrélées les unes avec les autres. Ce n’est pas le cas de day qui est non corrélée aux autres variables. La variable temperature est peu corrélée aux autres variables également et c’est la même chose pour precipitation.

Figure 2.10: Représentation de la corrélation des variables rateCar des voisins

On remarque sur ce graphique que rateCar et les variables rateCar des voisins sont très corrélées. Cela nous paraît normal car ce sont les mesures de la même quantité sur des points géographiquement proches.

3 Problèmes et modèles de prévision

Avec un jeu de données complet, nous pouvons enfin utiliser des algorithmes de Machine Learning pour faire de la prévision. Vu la complexité du jeu de données, il n’est pas évident de définir les axes de travail. Nous allons alors commencer par la définition de différents problèmes et de métriques de comparaison avant d’entraîner des modèles.

3.1 Définition des problèmes

Il y a deux dimensions principales dans nos données:

1. La première dimension du problème est temporelle : à quel horizon voulons-nous prédire ?
   1. Le cas une heure correspond à un problème du type “Google Maps” : si nous prenons la voiture maintenant, quelle sera la circulation dans les prochaines 60 minutes ?
   2. Etant donnée la circulation du jour précédent, la tâche devient intéressante pour les régulateurs des feux de circulation ou pour la police qui pourront ainsi éviter des bouchons.
   3. Si nous voulons prédire une semaine voire un mois en avance, le modèle devra comprendre davantage le comportement habituelle de la circulation routière. Ce sera plutôt un pari basé sur des fluctuations habituelles.
2. La deuxième dimension est reliée au réseau routier : pouvons-nous profiter de la structure du graphe pour améliorer la prédiction ?
   1. La question de base sera si on réussit à battre des benchmarks sur une seule arête.
   2. Puis nous pourrons réutiliser les “voisins” introduits pendant l’imputation.
   3. Finalement, une question est si nous réusissons à trouver un seul modèle décrivant toutes les arêtes en même temps, une sorte de prédiction multivariée.

Avant de nous lancer, nous précisons comment nous évaluons nos modèles et quels sont les benchmarks, les modèles naïfs, que nous voulons battre.

3.2 Evaluation des modèles

Train - test - split

On découpe notre jeu de données en 2 parties de proportion 2/3 et 1/3 : une partie apprentissage de 2014 à 2017 et une partie test de 2018 à 2019. Faire le découpage de cette manière assure qu’on est évalué sur les prédictions du futur plutôt que du passé. L’année 2020 est omise à cause du Covid, mais elle pourrait faire partie d’un autre projet sur le jeu de données en relation avec la pandémie !

Mesurer la performance

Pour évaluer les performances des modèles, nous considérerons le Root Mean Square Error (RMSE) qui représente la racine carrée du deuxième moment d’échantillonnage des différences entre les valeurs prédites et les valeurs observées. Une autre métrique souvent utilisée est le Mean Absolute Percentage Error (MAPE) qui n’est pas applicable à notre situation car un nombre important des valeurs de rateCar sont nulles et on ne peut pas diviser par elles.

Modèles naïfs à battre

Afin de pouvoir comparer nos modèles à une référence, une sorte de benchmark, on construit 3 modèles témoins. Ces modèles sont dit naïfs car extrêmement simpliste :

• naiveModel1 prévoit le taux d’occupation d’une heure \(t \in \{ 0, \dots, 23 \}\) d’un jour \(j \in \{ 1, \dots, 7\}\) d’un mois \(m \in \{ 1, \dots, 12 \}\) d’une année \(a \in \{ 2018, 2019 \}\) en moyennant les taux d’occupation de l’heure \(t\) du jour \(j\) du mois \(m\) pour toutes les années \(a \in \{ 2014, \dots, 2017 \}\).

• naiveModel2 fait de même en calculant la médiane et non la moyenne.

• naiveModel3 est simplement le taux d’occupation à l’heure précédente rateCar_LaggedHour.

Pour ces trois modèles simples, on calcule leurs erreurs RMSE pour chacune des 69 arêtes pour les moyenner après (voir le fichier naive_model.r pour le code). Ceci nous donne :

Table 3.1: Performances des modèles naïfs témoins

naiveModel1	naiveModel2	naiveModel3
5.844	5.899	3.955

Les 2 premiers modèles naïfs ont des scores similaires. Le troisième modèle naïf, qui consiste à prendre la variable rateCar de manière décalée, est meilleur. On interprète cette bonne performance par la très forte corrélation du taux de circulation d’une heure avec la suivante. Dans la suite, il sera intéressant de comparer ce modèle naïf avec des modèles plus complexes pour justifier leur utilisation.

Validation croisée

Pour la séléction de modèles et plus particulièrement pour la recherche de paramètres (d’un arbre CART par exemple), nous nous servirons beaucoup de la validation croisée. Ceci n’est pas évident pour les séries temporelles comme les données ne sont pas indépendant d’un ligne à l’autre. En plus, il vaut mieux éviter l’utilisation de données futurs pour prédire le passé. Entre ces problèmes et l’efficacité algorithmique, nous utilisons trois approches différentes:

• Validation croisée sur blocs regroupant des mois successifs: on divise les 4 années de données d’apprentissage en blocs de 3 mois ou 6 mois (correspondant à 16 et 8 “CV-folds” respectivement), ce qui englobe les différentes saisonnalités de nos données. C’est la méthode qui nous semble le plus pertinent pour les séries temporelles car les blocs sont à peu près indépendants entre eux et on respecte les saisonnalités dans le découpage.

• Validation croisée progressive: on fixe une fenêtre initiale d’apprentissage de 2 ans, que l’on incrémente dans l’ordre chronologique de 2 mois à chaque itération, pour entraîner le modèle et mesurer sa performance à l’horizon 1 (c’est-à-dire prédire le taux d’occupation de l’heure suivante). Pour utiliser cette méthode, on utilise le package caret avec le paramètre timeSlice.

Figure 3.1: Validation croisée progressive (“Progressive Cross-Validation” visité le 22.02.2022)

• Validation croisée habituelle : on mélange toutes les lignes et les découpe en K blocs. Cette approche habituelle pour des données i. i. d. est moins pertinente pour les séries temporelles, on utilise particulièrement le futur pour prédire le passé.

Si nous ignorons la première idée numériquement plus difficile, il semble plus naturel d’utiliser l’approche progressive pour notre type de données. Il y a pourtant le souci de largement surévaluer les données des premières deux années dans le score de CV, vu qu’elles sont représentées dans chaque couche (fold). En plus, les différentes couches seront alors fortement corrélées comme elles ont plus de la moitié de leur contenu en commun. Les deux idées ont alors leurs failles.

Finalement, lors de nos recherches de paramètres, nous avons observé que les deux méthodes sélectionnent des paramètres assez similaires. La plupart du temps, nous nous contentons alors de simplement utiliser la validation croisée habituelle avec 16 couches. Ce chiffre est inspiré par l’approche de blocs sur mois successifs. L’avantage de ceci est que le package caret permet une implémentation efficace. Nous notons ici que nous fixons toujours une seed du générateur aléatoire pour le mélange des lignes effectué par caret. Ainsi le lecteur pourra relancer le code lui-même.

3.3 Google Maps : Prédire la prochaine heure

Commençons par la tâche la plus simple : prédire la circulation dans les prochaines 60 minutes.

3.3.1 GAM pour une arête

Tout d’abord, nous nous restreignons à une arête isolée dans le sens que nous ne prenons pas en compte la circulation sur ces voisins. Notre objectif initial était de faire une simple régression linéaire, mais il y a quelques relations non linéaires. Si on essaie, par exemple, d’écrire le taux d’occupation d’une avenue en fonction de l’heure, le graphique suivant montre bien que les heures de pointe seraient mieux décrites par un polynôme que pour une droite :

Figure 3.2: Une journée typique sur l’arête

Par conséquent, nous travaillons plutôt sur des GAM que sur de simples régressions. Constituer un modèle additif est une tâche qui peut être longue. C’est pourquoi nous adapterons une approche greedy pour la sélection de modèles où nous rajoutons une variable explicative à la fois. Idéalement ce sera celle qui apporte le plus de puissance prédictive supplémentaire (forward search).

Pour prédire une heure dans le futur, il est incontestable que l’heure précédente contienne le plus d’informations possible. La variable rateCar_LaggedHour sera donc la première à être considérée. L’effet de l’occupation d’une heure à l’autre n’est pourtant pas homogène à travers les heures de la journée. En plus, comme nous voyons dans la figure inférieure gauche, il y a peut-être une relation linéaire entre les variables, qui a pourtant une grande variance:

Les 3 autres figures ci-dessus représentent la même comparaison rateCar_LaggedHours vs. rateCar, mais sur différentes heures. On a bien l’impression que des régressions linéaires pourront nous servir. Comparer le nuage de points à 8 heures avec ceux de 12 et de 13 heures nous donne envie de regrouper certaines heures. Ce phénomène était déjà observable dans le graphique 3.2 : dans l’heure de pointe matinale (morning rushhour), la circulation augmente beaucoup entre 7 et 8 heures tandis qu’entre 12 et 13 heures elle reste plutôt constante.

Afin de regrouper les heures de la journées en périodes semblables, nous procédons de manière heuristique : pour chaque heure de la journée, nous effectuons une régression linéaire et comparons les coefficients (intercept et slope). Ceci donne le tableau suivant.

(#tab:linReg_heureclustering)LinReg coefficients pour chaque heure

hour	0.000	1.000	2.000	3.000	4.000	5.000	6.000	7.000	8.000	9.000	10.000	11.000	12.000	13.000	14.000	15.000	16.000	17.000	18.000	19.000	20.000	21.000	22.000	23.000
intercept	-0.148	0.038	-0.218	0.039	0.007	0.163	0.414	1.155	-0.630	-1.185	1.691	2.253	1.688	2.278	0.756	1.778	2.442	1.527	1.230	1.529	2.020	0.709	2.252	1.675
slope	1.050	0.749	0.736	0.651	0.645	0.578	0.575	0.568	2.993	2.298	0.846	0.888	0.949	0.879	0.808	0.842	1.046	1.013	0.916	0.810	0.795	0.685	0.274	0.438

Les résultats sont congruents avec le graphique 3.2, par exemple l’explosion de slope à 8 heures où commence l’heure de pointe matinale. En nous basant sur ce tableau, nous regroupons de manière heuristique les heures en différents groupes: “night” de 0 à 6 heures, des groupes individuels de 7 à 9 heures (“7heures” par exemple), “noon” de 10 à 15 heures, “afternoon” de 16 à 17 heures, “evening_rush” de 18 à 20 heures et “evening” de 21 à 23 heures. Ainsi, chaque ligne du dataframe reçoit le nouveau attribut GRPS pour “groups.” Avec ces groupes, nous pouvons établir notre premier modèle mod1, une régression linéaire avec un “mixed effect” GRPS * rateCar_LaggedHour :

cluster_hours = function(h){
  if (h %in% c(0,1,2,3,4,5,6,7)){return("night")}
  if (h %in% c(7)){return("7heures")}
  if (h %in% c(8)){return("8heures")}
  if (h %in% c(9)){return("9heures")}
  if (h %in% c(10,11,12,13,14,15)){return("noon")}
  if (h %in% c(16,17)){return("afternoon")}
  if (h %in% c(18,19,20)){return("evening_rush")}
  if (h %in% c(21,22,23)){return("evening")}
  else{return("asfd")}
}
df_train$GRPS = as.factor(sapply(as.factor(df_train$hour), cluster_hours))
df_test$GRPS = as.factor(sapply(as.factor(df_test$hour), cluster_hours))
# for (group in unique(df_train$GRPS)){
#   indices = df_train$GRPS == group
#   plot(df_train[indices,]$rateCar_LaggedHour, df_train[indices,]$rateCar,
#        main = group, xlab = "rateCar_laggedHour", ylab = "rateCar")
# }
form1 = rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour
mod1 = lm(data = df_train, formula = form1)
summary(mod1)

## 
## Call:
## lm(formula = form1, data = df_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.480  -0.836  -0.301   0.749  51.016 
## 
## Coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         0.846441   0.024047   35.20   <2e-16 ***
## GRPS8heures:rateCar_LaggedHour      2.201884   0.038810   56.73   <2e-16 ***
## GRPS9heures:rateCar_LaggedHour      1.935035   0.013804  140.18   <2e-16 ***
## GRPSafternoon:rateCar_LaggedHour    1.097172   0.003656  300.10   <2e-16 ***
## GRPSevening:rateCar_LaggedHour      0.572480   0.003810  150.26   <2e-16 ***
## GRPSevening_rush:rateCar_LaggedHour 0.881992   0.002559  344.71   <2e-16 ***
## GRPSnight:rateCar_LaggedHour        0.563815   0.011286   49.96   <2e-16 ***
## GRPSnoon:rateCar_LaggedHour         0.928103   0.002810  330.23   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.421 on 35052 degrees of freedom
## Multiple R-squared:  0.8753, Adjusted R-squared:  0.8753 
## F-statistic: 3.516e+04 on 7 and 35052 DF,  p-value: < 2.2e-16

Nous n’avons gardé que les effets croisés dans notre formule : notons que toutes les variables ont une contribution significative d’après les tests de student. Nous calculons le score de validation croisée pour ce modèle en utilisant la stratégie sur blocs regroupant des mois successifs, ce qui est facile à implémenter ici.

CV_gam = function(formula, df_train, no_folds){
  K = no_folds 
  data = df_train
  folds <- cut(seq(1,nrow(data)),breaks=K,labels=FALSE)
  fold_scores = c()
  start_time = Sys.time()
  for(i in 1:K){
    testIndexes <- which(folds==i,arr.ind=TRUE)
    testData <- data[testIndexes, ]
    trainData <- data[-testIndexes, ]
    g = gam(data=trainData, formula=formula)
    Y_test = testData$rateCar
    Y_predict = predict(g, newdata=subset(testData, select = -c(rateCar)))
    fold_scores = c(fold_scores, rmse(Y_test, Y_predict))
  }
  exec_time = Sys.time()-start_time
  score = mean(fold_scores)
  ret = data.frame(score = score, no_folds = no_folds, exec_time = exec_time, 
                   formula = paste(as.character(formula)[c(2,1,3)], collapse = " "))
  return(ret)
}
knitr::kable(CV_gam(form1, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")

(#tab:CV_mod1)Model score from CV

score	no_folds	exec_time	formula
2.42125	8	3.364563 secs	rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour

Ce premier modèle simple bat alors largement le meilleur modèle naïf.

Afin de choisir la prochaine variable à rajouter pour expliquer rateCar, nous les explorons une par une :

1. hour : Comme vu dans le graphique 3.2, la relation entre l’heure et le taux d’occupation n’est pas du tout linéaire. Pour trouver une fonction lisse adéquate, nous utilisons une base de splines cycliques car les heures 0 et 24 coïncident. Pour les nœuds, nous pouvons augmenter le paramètre k jusqu’au plus la valeur 24. En augmentant ce k et en regardant edf, les degrés de liberté estimés, nous nous rendons compte que k = 24 donne toujours un edf de presque 22, ce qui nous fait choisir cette valeur de k :

form1_hour = rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour + s(hour, bs="cc", k = 24, by=weekendsIndicator)
mod1_hour = gam(form1_hour, data=df_train)
summary(mod1_hour)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour + 
##     s(hour, bs = "cc", k = 24, by = weekendsIndicator)
## 
## Parametric coefficients:
##                                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                         1.102448   0.027395   40.24   <2e-16 ***
## GRPS8heures:rateCar_LaggedHour      2.302777   0.108018   21.32   <2e-16 ***
## GRPS9heures:rateCar_LaggedHour      2.344126   0.033255   70.49   <2e-16 ***
## GRPSafternoon:rateCar_LaggedHour    1.026940   0.007463  137.60   <2e-16 ***
## GRPSevening:rateCar_LaggedHour      0.525643   0.006972   75.39   <2e-16 ***
## GRPSevening_rush:rateCar_LaggedHour 0.819604   0.004874  168.17   <2e-16 ***
## GRPSnight:rateCar_LaggedHour        0.656110   0.016709   39.27   <2e-16 ***
## GRPSnoon:rateCar_LaggedHour         0.867199   0.004684  185.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                              edf Ref.df      F p-value    
## s(hour):weekendsIndicator0 21.94     22 206.40  <2e-16 ***
## s(hour):weekendsIndicator1 21.48     22  75.86  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.894   Deviance explained = 89.4%
## GCV = 5.0039  Scale est. = 4.9965    n = 35060

knitr::kable(CV_gam(form1_hour, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")

(#tab:CV_mod1_hour)Model score from CV

score	no_folds	exec_time	formula
2.24125	8	14.60161 secs	rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour + s(hour, bs = “cc,” k = 24, by = weekendsIndicator)

D’après le summary de ce modèle GAM, le score de validation croisée a clairement été amélioré. N’oublions pas de mentionner que nous avons aussi rajouté un effet croisé de hour avec weekendsIndicator. Ceci est dû au simple fait que la circulation à des courbes différentes le week-end, ce qui est confirmé par les p-valeurs des deux variables. Dans le graphique ci-dessous, nous visualisons les fonctions lisses obtenus par GAM, à gauche pendant la semaine et à droite le week-end :

draw(mod1_hour)

D’après ce graphique, la variable hour agit comme terme correctif pour le modèle initial des regréssions par GRPS.

Pour les variables suivantes, nous donnerons moins de commentaires et n’afficherons plus les summary ou des plots de fonctions. L’approche restant la même.

2. weekdays : numérotées de 1 à 7, nous utilisons encore une fois une base de splines cycliques. Le choix de k s’effectue comme avant et nous essayons d’ajouter un effet croisé avec l’une des variables HolidaysIndicator, i.e. summerHolidaysIndicator etc. Les indicateurs de vacances n’ont jamais des effets significatifs, raison pour laquelle nous rejetons cette approche. Au final, nous obtenons le score de CV suivant, qui n’améliore qu’à peine le premier modèle.

form1_weekdays = rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour + s(weekdays, bs="cc", k = 7) 
#mod1_weekdays = gam(form1_weekdays, data=df_train)
#summary(mod1_weekdays)
#plot(mod1_weekdays, residuals=T, rug=T, se=F, pch=20)
knitr::kable(CV_gam(form1_weekdays, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")

(#tab:CV_mod1_wd)Model score from CV

score	no_folds	exec_time	formula
2.41625	8	3.95377 secs	rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour + s(weekdays, bs = “cc,” k = 7)

3. state : Cette variable catégorielle indiquant l’état de la rue a été partiellement détruite pendant la simplification du dataset où nous avons pris des moyennes, résultant en state \(\not \in \{1,2,3,4\}\) dans notre data set. Comme attendu, le score de validation croisée ne s’améliore pas et nous allons ignorer cette variable dans la suite.

form1_state = rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour + s(state, bs="cr", k = 5)
#mod1_state = gam(form1_state, data=df_train)
#summary(mod1_state)
knitr::kable(CV_gam(form1_state, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")

(#tab:CV_mod1_state)Model score from CV

score	no_folds	exec_time	formula
2.42125	8	2.131716 secs	rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour + s(state, bs = “cr,” k = 5)

4. toy (time of year): Y a-t-il des tendances dans la circulation à travers l’année ? toy est une autre variable exigeant des splines cycliques, mais comme nous voyons après tuning, la rajouter n’améliore pas le score de validation croisée :

form1_toy = rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour + s(toy, bs="cc", k = 20)
#mod1_toy = gam(form1_toy, data=df_train)
#summary(mod1_toy)
knitr::kable(CV_gam(form1_toy, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")

(#tab:CV_mod1_toy)Model score from CV

score	no_folds	exec_time	formula
2.42	8	4.610196 secs	rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour + s(toy, bs = “cc,” k = 20)

4. nbCar_LaggedHour: cette variable pose la question intéressante de savoir si le nombre de voitures peut nous aider à prédire le taux d’occupation. Comme avec rateCar_LaggedHour, nous pouvons observer un effet hétérogène de nbCar_LaggedHour sur rateCar à travers la journée. Voici, par exemple, les nuages de points à 8 et à 14 heures :

g_tuple = list(); i=1 
for (h in c(8,14)){
  g = ggplot(data=df_train[(df_train$hour == h),], aes(y = rateCar, x = nbCar_LaggedHour))+ 
  geom_point() + xlab("nbCar_laggedHour") +
  ylab("rateCar (présent)") + 
  labs(title = paste0("Heure ", h))
  g_tuple[[i]] = g
  i = 2
}
plot_grid(g_tuple[[1]], g_tuple[[2]], nrow = 1, ncol = 2)

Contrairement à l’effet linéaire que nous avons observé avec rateCar_LaggedHour, il semble y avoir des non-linéarités dans les nuages de points. Par conséquent, nous utiliserons une fonction à base de splines cubiques. Nous avons testé par CV que des thin plate splines ou des P-splines ont la même performance. Pour considérer un effet différent à travers les heures de la journée, nous réutilisons les groupes d’heures définies avant en ajoutant le paramètre by = GRPS dans gam. La validation croisée nous donne le score

# Have tried bs \in c('cr', 'tp', 'ps', 'ad')
form1_nbCarHour = rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour + s(nbCar_LaggedHour, bs="cr", k = 15, by=GRPS)
#mod1_nbCarHour = gam(form1_nbCarHour, data=df_train)
#summary(mod1_nbCarHour)
#plot(mod1_nbCarHour, residuals=T, rug=T, se=T, pch=20)
#draw(mod1_nbCarHour, residuals=T)
knitr::kable(CV_gam(form1_nbCarHour, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")

(#tab:CV_mod1_nbCarHour)Model score from CV

score	no_folds	exec_time	formula
2.3475	8	35.76525 secs	rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour + s(nbCar_LaggedHour, bs = “cr,” k = 15, by = GRPS)

C’est une amélioration significative du premier modèle, même si heure a amélioré encore plus la performance.

Nous ne détaillons pas ici la suite de nos reflexions. Les variables heure et nbCar_LaggedHour continuent d’avoir la contribution la plus forte. Dans le fichier 1_heure/GAM_listEdges_1hour.R, le lecteur peut avoir un aperçu du feature engineering sur precipitation, les scores de chaque modèle, etc. Tout ceci nous mène finalement au modèle final dont nous affichons la formule et le score de validation croisée :

form_final = rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour + s(hour, bs="cc", k = 24, by=weekendsIndicator) + s(nbCar_LaggedHour, bs="cr", k = 15, by=GRPS)
knitr::kable(CV_gam(form_final, df_train, no_folds = 8), align = "llll", caption = "Model score from CV")

(#tab:CV_modfinal)Model score from CV

score	no_folds	exec_time	formula
2.21125	8	1.522581 mins	rateCar ~ GRPS * rateCar_LaggedHour - GRPS - rateCar_LaggedHour + s(hour, bs = “cc,” k = 24, by = weekendsIndicator) + s(nbCar_LaggedHour, bs = “cr,” k = 15, by = GRPS)

Malgré le fait d’avoir travaillé sur une seule arête, nous sommes curieux si l’efficacité de ce modèle se généralise aux autres rues de Paris. Par ailleurs, nous avons jusqu’à maintenant seulement calculé les erreurs de validation croisée. Comme la sélection de modèles est finie, nous voulons bien tester la capacité de généralisation du modèle final sur le test set.

Pour l’arête arbitraire que nous avons choisie pour cette analyse, concorde - saint michel, nous obtenons un RMSE de 2.24, pas loin du score de validation croisée. Sur cette arête, le test set ne semble alors pas se comporter très différemment du train set. Pourtant, la capacité de généralisation sur les autres arêtes est limitée comme montre l’histogramme de la répartition des erreurs quadratiques sur le test set à travers les 69 arêtes :

scores_edges_final = readRDS("../1_heure/scores_edges_GAM_final.rds")
hist(scores_edges_final, breaks = 30, main = "Histogram of test scores across 69 edges")

Il existe des arêtes où le score est meilleur que sur concorde - saint michel, mais plusieurs arêtes font exploser la moyenne de ces erreurs à 3.89, à peine meilleure que les modèles naïfs.

Avant de conclure que GAM généralise mal à travers le graphe, nous revenons brièvement sur le premier modèle qui n’utilisait que rateCar_LaggedHour comme variable explicative. Si nous faisons pareil en calculant le RMSE sur le test set, l’erreur à travers des arêtes est 3.35, donc significativement meilleure. Le groupement des heures, le rôle de nbCar, semble causer de l’overfitting. Un bon GAM doit alors être bricolé pour chaque arête individuellement.

3.3.2 Arbres aléatoires

Après cet exercise de bricolage avec GAM, nous passons à d’autres modèles dont le tuning pourra se faire de manière plus automatique. Pour commencer, nous nous intéressons à des arbres CART et nous nous servons du package rpart pour les implémenter. Au lieu d’examiner la contribution de chaque variable, le choix des \(k\) etc., il suffit de tuner quelques paramètres de l’algorithme. Vu que nous ferons un tuning étendu pour les forêts aléatoires après, nous nous focalisons ici sur un seul paramètre qui est de plus haute importance : le “complexity parameter” cp. Dans la construction de l’arbre, si le gain de performance après une découpe n’est pas meilleur d’un facteur de cp, alors la découpe n’est pas réalisée. Plus cp est grand, moins l’arbre sera alors complexe, d’où le nom.

Pour optimiser cp, nous utilisons de la validation croisée et une recherche sur deux grilles : la première de 0.0 à 0.1 et la seconde, plus fine, de 0.0 à 0.01. Nous affichons les scores de CV obtenues (voir le fichier 1_heure/tree.R pour le code) :

Figure 3.3: Recherche du paramètre cp optimal

La valeur optimale d’après cette CV est alors cp=0, c’est-à-dire que l’on accepte toutes les découpes et que notre arbre est un arbre profond. On maintient à leur valeur par défaut les autres paramètres de CART et on entraîne le modèle sur notre train dataset. Nous faisons ceci dans deux cas différents qui représentent deux dimensions locales différentes : rappelons-nous que pour la complétion, nous avions rajouté des informations sur les voisins pour enrichir les données. Pour voir si ces mêmes informations nous servent aussi maintenant, nous entraînons un arbre avec et un CART sans données voisines. Voici les scores sur le test set dans les deux cas :

Table 3.2: Performances des arbres aléatoires

Arbre avec voisins	Arbre sans voisins
3.783	3.582

Premièrement, on note dans les deux cas une performance légèrement meilleure que le modèle naïf qui recopie l’heure d’avant. Ce qui est frappant est que l’ajout d’information sur les voisins baisse le score, indiquant qu’il y a soit de l’overfitting, soit un problème de grande dimension qui perturbe le bon fonctionnement des arbres aléatoires. En tout cas, ce résultat sème le doute sur l’utilité de la composante spatiale de nos données. Mais peut-être, il faut juste utiliser des techniques plus puissantes, comme les forêts aléatoires qui rassemblent plusieurs arbres en utilisant du bagging et de l’aléatoire, c’est ce que nous allons étudier.

3.3.3 Forêts aléatoires

Pour cette méthode d’ensemble basé sur les arbres de Breimann, le package ranger nous servira à construire les forêts et le package vip aidera à visualiser l’importance des variables. Comme pour les arbres, il y a des paramètres à optimiser pour obtenir le meilleur prédicteur possible. Nous nous focalisons sur deux paramètres que nous optimisons l’un après l’autre. C’est une approche moins prometteuse qu’une recherche sur grille, mais elle diminue le temps de calcul et nous permet de visualiser la recherche de paramètres :

• Le premier paramètre à optimiser est min.node.size (minimal node size). L’algorithme ne découpe pas une feuille dont le nombre d’éléments est inférieur ou égal à sa valeur. Comme cp précédemment, ce paramètre contrôle la taille des arbres.

Figure 3.4: Recherche du paramètre min.node.size optimal

On remarque que l’erreur décroît jusque la valeur 1 avec un pic négatif à 9. Choisir min.node.size=1 donnerait des arbres de profondeur maximale, rendant le modèle susceptible à sur-apprendre. Par conséquent nous nous contentons de choisir 9, malgré les fluctuations importantes des CV-scores dans le graphique à droite. Ces flucuations pourraient être dues au fait qu’ici on a fait une validation croisée habituelle (voir nos explications sur CV en haut) qui mélange les lignes corrélées avant de les séparer en blocs. Dans un cadre iid, nous nous attendrions à la forme typique d’un “coude” par exemple.

• Le deuxième paramètre est mtry : le nombre de variables, tirées au hasard, qui sont considérées pour la construction de chaque arbre. Rendre ce paramètre plus petit réduit alors la corrélation entre les arbres de la forêt. Pour l’optimiser, nous utilisons deux de nos approches CV :

Figure 3.5: Recherche du paramètre mtry optimal

Celle-ci est l’unique fois où nous affichons la méthode de validation croisée par timeslice. On remarque que l’allure globale des courbes sont identiques mais que le coude est plus marquée avec timeslice. On choisit de poser \(\textbf{mtry} = 10\), qui semble être un bon choix d’après les deux graphiques.

Le dernier paramètre à déterminer est ntree, le nombre d’arbres dans la forêts. Ce n’est pas un paramètre à optimiser comme les 2 précédents car le sur-apprentissage dans les forêts aléatoires est peu lié au nombre d’arbres. On optimise ntree afin de déterminer à partir de quel nombre d’arbre on gaspille du temps d’exécution pour un faible gain de performance.

Figure 3.6: Recherche du paramètre ntree ‘’optimal’’

On voit que le temps d’exécution est linéaire selon le nombre d’arbres et que la courbe de l’erreur marque un coude à partir de \(150\), valeur que l’on choisit pour ntree.

Avec ces trois paramètres choisis, nous entraînons la forêt aléatoire et l’évaluons sur le test set. Comme avec les arbres, nous faisons ceci deux fois, une fois avec voisins et une fois sans.

Table 3.3: Performances des forêts aléatoires

Forêt aléatoire avec voisins	Forêt aléatoire sans voisins
3.092	3.129

Comme attendu, les performances des forêts aléatoires sont généralement meilleures que celles des arbres aléatoires. Pourtant, contrairement à ces derniers, les arêtes voisins semblent avoir un effet plutôt positif. Si nous nous intéressons davantage au rôle des voisins, nous pouvons exploiter le fait que les forêts aléatoires livrent des scores d’importance pour chaque variable. Sur l’arête pont amont - pont austerlitz, par exemple, le ranking d’importance des variables est le suivant :

Figure 3.7: Scores d’importance dans l’ordre décroissante (15 premières arêtes) pour l’arête ‘’pont amont - pont austerlitz’’

Sur cette arête, les voisins dominent les scores d’importance. Ceci nous donne une motivation pour davantage explorer le rôle des voisins et ainsi le rôle de la dimension spatiale du problème de prédiction de la circulation :

Une question immédiate est quelles arêtes profitent le plus des données de leurs voisins. Nous essayons de visualiser ceci sur le graphique suivant où pour chaque arête du graphe routier, nous visualisons la différence des erreurs (sur test set) avec et sans voisins. Le trait vertical noir pointillé marque la séparation entre les rues à l’intérieur Paris et les parties du périphérique. Ces dernières se trouvent à droite de la ligne. Le trait horizontal rouge pointillé est un seuil arbitraire à partir duquel nous interprétons que l’ajout des voisins a un impact significativement positif sur la performance de la forêt. Les arêtes qui ont une valeur sous ce seuil (et qui profitent donc de l’inclusion de leurs voisins) sont représentées par une croix rouge, les autres par une croix noire.

Seules 12 arêtes sur 69 ont une performance “significativement” meilleure lorsqu’on ajoute les voisins. De plus, pour les arêtes du périphérique, on remarque qu’inclure les voisins a, à 2 exceptions près, un effet néfaste sur la performance. Même en général, la performance sur le périphérique est moins bonne que dans Paris. En effet, les scores pour Paris sont 2.747 (avec voisins) et 2.898 (sans voisins), et les scores du périphérique sont 4.069 (avec voisins) et 3.787 (sans voisins).

3.3.4 Analyse des résidus

Nous venons de voir les arbres CART et les forêts aléatoires, tous les deux avec une performance peu meilleure qu’un simple GAM. On se pose alors la question s’il reste peut-être de l’information dans les résidus : comme nous travaillons sur des séries temporelles, SARIMA ou un lissage exponentiel pourrait découvrir de l’information qui est cachée pour un arbre CART.

Supposons alors qu’il reste de l’information dans les résidus des forêts aléatoires. Nous choisissons ce modèle avec l’espoir qu’on pourra ainsi atteindre une même performance qu’avec GAM. L’information obtenue grâce à l’analyse des résidus sera ajoutée comme terme additif au modèle de base. Vu la taille de nos données et les temps d’exécution reliés, nous nous restreignons dans cette partie à une arête du graphe.

Afin d’intégrer la saisonnalité des résidus dans le modèle, nous nous servons de SARIMA, méthode pour laquelle il faut déterminer plusieurs paramètres. Pour ce faire, nous aimerions nous servir de la fonction auto.arima qui trouve des paramètres optimaux. Pour ce faire, cette fonction exige des bornes supérieures sur les paramètres p, q, P et Q. Nous déterminons ces bornes avec l’approche Box - Jenkins :

Parmi les saisonnalités (journalières, hebdomadaires, annuelles), on choisit de traiter la saisonnalité journalière qui semble être la plus importante. On pose donc s=24, D=1 et on différentie la série temporelle avec lag=24. Cela fait, nous affichons les auto-corrélogrammes et les auto-corrélogrammes partiels. Les deux diagrammes du haut se basent sur la première saison, de 1 à 24 observations, et nous servent à déterminer q et p. Les deux du bas présentent la situation à une échelle temporelle plus grande, 20 et 17 saisons respectivement, et servent à déterminer P et Q.

Figure 3.8: Analyse de la série temporelle de la première arête

En haut à gauche, nous observons tout d’abord une décroissance relativement forte, ce qui nous fait choisir d=0. Ensuite, on déduit les valeurs p.max=9 et q.max=15 des 2 figures du haut et P.max=15 et Q.max=1 des 2 figures du bas (en comptant les pics). La fonction auto.arima utilise ces paramètres pour choisir le modèle \(SARIMA(9,0,0)(4,1,0)[24]\). Ce nouveau modèle, qui a pour but d’expliquer les résidus de random forest, a lui-même des résidus. Ce sont ces résidus que nous analysons ensuite avec un test de Box-Pierce:

Table 3.4: Test de Box-Pierce

lags	statistics	pvalue
1	0.0129069	0.9095483
24	781.0112511	0.0000000

Avec ces p-valeurs, nous jugeons que l’auto-corrélation sur la prochaine heure a bien été éliminée mais il reste clairement de l’auto-corrélation journalière. Ce résultat est confirmé par un plot des résidus qui montre d’ailleurs que la saisonnalité hebdomadaire n’a pas été éliminée non plus.

Nous avons essayé de modifier les paramètres à la main et en utilisant les critères AIC/BIC pour trouver un meilleur modèle SARIMA. Malheureusement, augmenter les valeurs de Q, p etc. est impossible car cela a un coût en mémoire trop important (nous parlons de vecteurs 13GB pour P=16 par exemple). Il y a évidemment d’autres manières d’optimiser ces paramètres, mais nous n’avions pas pris le temps de creuser davantage. En particulier, parce que chaque arête du graphe nécessiterait une optimisation individuelle, si nous voulons obtenir un modèle complet.

Si nous essayons, malgré l’auto-corrélation qui reste, de combiner le modèle SARIMA avec la forêt aléatoire initiale, nous obtenons un RMSE de 8.7. La technique ici était de faire des prédictions séquentielles, comme vu dans le cours avec la fonction rollArima. Pour les détails de nos calculs, nous renvoyons au fichier 1_heure/rfArima. Comme c’est la pire performance parmi nos modèles, nous ne voulons pas insister dessus.

En recherche d’une alternative à SARIMA, nous avons aussi essayé des méthodes de lissage exponentiel, les fonctions ets et HoltWinters. Mais comme précédemment, nous étions soit bloqués par des restrictions de performance, soit les scores de prédiction étaient terribles.

Pour conclure sur l’idée d’analyser les résidus des forêts aléatoires, on ne parvient pas à un résultat satisfaisant. Avec une hardware plus puissante et surtout plus d’investissement de temps, nous avons l’espoir d’obtenir une amélioration de la prédiction.

3.3.5 Boosting

Après les forêts aléatoires, nous voulons essayer un dernier algorithme de Machine Learning pour le scénario de prévision “Google Maps” : le boosting. Plus précisément, nous nous servons du package XGBoost qui est réputé pour ses très bonnes performances.

On construit un modèle à base d’arbres (booster = bgtree) et nous fixons arbitrairement le taux d’apprentissage eta = 0.1 et le paramètre d’avancement gamma = 0.2. Faire ces choix nous évite d’optimiser deux paramètres de plus, mais ce temps de calcul pourrait évidemment coûter en termes de performance.

Il y a tout un ensemble de paramètres que nous optimisons avec la package caret. Une validation croisée réalisée sur 5 blocs sert à évaluer les performances. Un nombre de blocs plus élevé serait désirable, mais le temps d’exécution avec 5 blocs est déjà très important (environ 3 heures). La recherche de paramètres est effectuée sur la grille suivante, qui montre également aussi quelles valeurs nous optimisons.

Table 3.5: Grille de recherche des paramètres optimaux

nrounds	max_depth	colsample_bytree	eta	gamma	min_child_weight	subsample
100	3, 5, 7, 10	0.70, 0.75, 0.80, 0.85, 0.90	0.1	0.2	1,3,5,7,9	0.70, 0.75, 0.80, 0.85, 0.90

La combinaison optimale trouvée par validation croisée est

Table 3.6: Paramètres optimaux

nrounds	max_depth	colsample_bytree	eta	gamma	min_child_weight	subsample
100	7	0.90	0.1	0.2	1	0.90

On optimise ensuite le nombre d’itérations nrounds à l’aide de la validation croisée sur 5 blocs (fonction xgb.cv) et d’un critère d’arrêt à partir de 5 itérations sans amélioration significative de la performance du modèle. On obtient la valeur de 133.

Après prédiction des 69 arêtes, on obtient les scores suivants

Table 3.7: Performances des méthodes de boosting

Boosting avec voisins	Boosting sans voisins
3.295	3.284

Les 2 scores sont très proches et on ne peut déterminer si l’ajout des voisins est pertinent dans cette méthode de boosting. Par rapport aux autres méthodes, ces scores sont meilleurs que ceux des arbres aléatoires mais moins bons que ceux des forêts aléatoires. C’est surprenant car incongruent avec nos attentes que le boosting ait une meilleure performance.

Une raison pour cela peut être notre approche d’optimisation qui n’est pas parfaite. On aurait aimé optimiser les paramètres gamma, eta et aussi les paramètres de régularisation lambda et alpha. Cela aurait été utile, surtout parce que notre boosting semble sur-apprendre. Une idée serait d’utiliser une méthode d’optimisation moins coûteuse comme de la recherche aléatoire de paramètre dans la grille ou une méthode bayésienne : idées que nous ne poursuivons pas davantage.

3.4 Agrégation d’experts

On conclut notre travail de prédiction par l’agrégation de tous les modèles que l’on a construits. Pour ce faire, nous utilisons dans le package opera le modèle MLpol avec la fonction de perte square. La nature séquentielle de nos données se prête bien à l’utilisation d’une agrégation en ligne des prédicteurs.

On présente rapidemment les modèles retenus, with indique qu’on a ajouté les voisins et without que l’on ne l’a pas fait.

• naive est le modèle naïf consistant à prendre le taux d’occupation de l’heure précédente

• cart_with et cart_without sont les arbres aléatoires

• rf_with et rf_without sont les forêts aléatoires

• xgb_with et xgb_without sont les modèles de boosting

• gam_simple est le modèle GAM qui n’utilise que rateCar_LaggedHour tandis que gam_full est le GAM avec plus de termes.

Le modèle agrégé a de meilleurs performances que le meilleur des experts: 2.806. Etudions le à l’aide des 2 figures suivantes:

Figure 3.9: Poids des experts

Ce graphique des poids selon les prédictions montre l’importance des experts au fur et à mesure des prédictions. On voit que certains experts contribuent beaucoup à l’agrégation, rf_with et xgb_without notamment, alors que d’autres contribuent très peu, cart_with par exemple. On voit que ces contributions ne sont pas uniformes dans le temps. Celle de xgb_without devient importante uniquement à partir de la 4093-ème prédictions. Dans l’ensemble, tous les prédicteurs ont des poids significatifs sur une certaine plage de prédictions. On en conclut que tous nos experts sont utiles à l’agrégation et sont complémentaires (hormis cart_with).

Figure 3.10: Pertes moyennes selon la perte L2

Ce graphique représente l’erreur moyenne de chaque expert, du prédicteur uniform formé de la moyenne des experts et du prédicteur agrégé MLpol. L’erreur est mesurée selon la perte L2, le carré du RMSE. On peut ainsi comparer les erreurs de chaque prédicteur. Tout d’abord, le prédicteur agrégé MLpol est meilleur que le meilleur des experts, ce qui confirme la pertinence dans notre cas de l’agrégation d’experts. On remarque également que la simple moyenne de nos experts uniform bat chaque expert. Concernant cart_with, l’expert dont le poids est le plus faible (nul sur 90% des prédictions), on voit qu’il a une très mauvaise performance par rapport aux autres experts. Cela concorde avec le fait que l’algorithme lui a donné un poids négligeable. Cependant, il faut faire attention à ne pas généraliser cette corrélation: rf_without a un poids assez faible mais une meilleure performance que rf_with, qui est l’expert avec le plus de poids.

3.5 Un modèle pour tout le graphe ?

Jusqu’ici, nous nous sommes beaucoup focalisés sur la tâche de prédiction pour une seule arête. Même si les voisins (plus précisément la circulation sur les arêtes adjacentes) font intervenir la dimension spatiale du problème de prédiction, ce sont des modèles très locaux plutôt que globaux.

On se demande alors : est-ce qu’on peut trouver un seul modèle qui décrit tout le graphe ? C’est la question que nous nous sommes posés, inspirés par la régression multivarié et l’espoir que le graphe contient plus d’information que les arêtes individuelles. Pour convertir cela en une tâche de prédiction, il faut d’abord réunir les 69 arêtes dans un grand dataframe que nous appelons “oneDF.” Chaque ligne de cet énorme jeu de données contient

• Le nom d’une arête (edgename) et la rateCar correspondante à l’heure actuelle.
• Les taux d’occupation et le nombre de voitures de toutes les autres arêtes à l’heure précédente (69 x 2 = 138 colonnes).
• Les variables statiques : la météo, l’heure, etc.

Ainsi, nous répliquons les 4 ans de données train 69 fois, ce qui résulte en un dataframe de 1.2 Go. Pourtant, cette approche nous permet d’entraîner plusieurs modèles en même temps avec un seul algorithme, l’espoir étant que le nom de l’arête en tant que variable catégorielle suffise pour déclencher la bonne prédiction.

Pour faciliter l’apprentissage, nous rajoutons par ailleurs deux variables reliées à la dimension spatiale: un indicateur si l’arête appartient au périphérique et une variable catégorielle qui divise Paris en sud, centre, nord-est et ouest.

Reste à choisir un algorithme qui arrive à gérer cette quantité de variables sans trop de calibration ni de perte de performance. Nous reprenons alors le boosting. Cependant, xgboost ne marche pas avec toutes nos variables catégorielles. L’alternative LGBM avait besoin de trop de RAM, ce qui nous motive à utiliser catboost.

Pas question de tuner les paramètres avec notre énorme dataset, nous reprenons alors les paramètres présentés dans le cours. Ensuite, nous nous forçons de réduire la taille du dataset : nous retirons quelques variables jusqu’ici inutiles comme state et ignorons deux tiers des arêtes. Enfin, catboost peut apprendre en un temps raisonnable d’une heure et demi (et il peut gérer les données qui ne rentraient pas dans la mémoire vive).

Pendant l’apprentissage, catboost calcule des scores d’importance pour chaque variable. Nous affichons les cinq variables avec les scores les plus élevés :

##                             edgename       rateCar_pontavalTOportemaillot 
##                            20.059169                            18.773107 
##                     periph_indicator rateCar_portechapelleTOporteasnieres 
##                            15.627989                             9.647768 
##           rateCar_chateletTOconcorde 
##                             6.656422

Notre idée pour le modèle est partiellement validée par ces scores : le nom des arêtes est le facteur le plus important pour déterminer la circulation et le fait d’appartenir au périphérique joue également un rôle fort. Mais est-ce que ce modèle pseudo-multivarié sert aussi à faire des prédiction ? L’erreur sur test set donne la réponse :

Une erreur moyenne de 5.07 est mauvaise en comparaison avec nos autres modèles, mais elle bat néanmoins deux modèles naïfs. Considérant le manque de calibration et la complexité du dataset, nous sommes finalement agréablement surpris par ce score. Par exemple, pour une journée sur l’arête ‘saint michel - concorde,’ on obtient

4 Conclusion

Un projet de Machine Learning avec un jeu de données aussi énorme que celui sur la circulation à Paris a deux visages.

D’une part, il faut le dompter : agréger les données de milliers de capteurs routiers, compléter les valeurs manquantes et extraire de l’information intéressante de ses données mises à jour quotidiennement. Une étape qui exige de la patience, de la minutie et de l’intuition statistique.

D’autre part, le jeu de données permet de poser un nombre presque infini de questions. Nous n’avons que gratté la superficie en étudiant plusieurs modèles visant à prédire la circulation de la prochaine heure, souvent arête par arête. A travers les voisins et à travers un modèle multivarié, nous avons essayé d’ajouter une dimension spatiale à nos problèmes. Nous sommes persuadés que la structure de graphe du réseau routier permettrait des modèles beaucoup plus poussés. Ensuite, il y a la dimension temporelle, la prédiction du prochain jour, de la prochaine semaine… Quelle précision peut-on maintenir en s’éloignant de l’heure dont on veut connaître la densité de circulation ? A partir de quel horizon la tâche devient très difficile pour le modèle naïf (heure précédente) et donc quand justifier l’utilisation de modèles plus complexes ? Réaliser des prédictions pour d’autres horisons est facilemment faisable à partir de nos fichiers de code, c’est d’ailleurs ce que l’on avait prévu un temps pour un horizon 24 (la même heure du jour suivant).

Avec toutes ces questions, nous pourrions passer des semaines, voire des mois sur le projet. Mais il faut arrêter à un moment et tirer un bilan. Nous nous sommes principalement intéressés au problème Google Maps, la prédiction de la prochaine heure, pour une ou plusieurs arêtes. Au final, c’est le modèle GAM qui remporte la meilleure performance sur une seule arête, après une calibration spécifique. Pourtant, sur l’ensemble des arêtes, nous obtenons les meilleurs scores avec random forest. C’est peu surprenant vu que ce dernier est très flexible et adaptatif aux jeux de données, tandis que GAM est plutôt rigide, adapté à un unique jeu de données.

Nous avons vu boosting et SARIMA échouer, partiellement dûe à la grande quantité de données et les restrictions computationnelles de nos ordinateurs. Tout au long du projet, il nous fallait faire des compromis entre temps d’exécution, mémoire et complétude de notre analyse. Au fur et à mesure de l’agrégation par exemple, nous nous sommes habitués à faire des choix qui réduisent la complexité du data set tout en gardant de l’interprétabilité et du réalisme.

Si nous tirons une leçon principale de ce projet, c’est le bon compromis qu’il faut souvent chercher : utiliser caret ou une validation croisée manuelle, optimiser un paramètre ou pas, stocker un jeu de donnés de 1.5 Go ou simplifier la tâche ? Notre objectif n’était pas de faire des choix parfaits, mais d’être conscients des possibles erreurs et de bien documenter nos réflexions.

Références

“Covid Index.” visité le 26.12.2021. https://www.bsg.ox.ac.uk/research/research-projects/covid-19-government-response-tracker.

“Données Du Comptage Routier.” visité le 26.12.2021. https://parisdata.opendatasoft.com/explore/dataset/comptages-routiers-permanents-historique/information/.

“Jours Fériés.” visité le 26.12.2021. https://www.data.gouv.fr/en/datasets/jours-feries-en-france/.

“Les Archives Du Calendrier Scolaire.” visité le 26.12.2021. https://www.education.gouv.fr/les-archives-du-calendrier-scolaire-12449.

“Observations Météorologiques Historiques En France.” visité le 26.12.2021. https://public.opendatasoft.com/explore/dataset/donnees-synop-essentielles-omm/table/?sort=date&q=%C3%AEle+de+france&refine.nom=ORLY&q.timerange.date=date:%5B2013-12-31T23:00:00Z+TO+2021-01-01T22:59:59Z%5D.

“Paris Data.” visité le 26.12.2021. https://parisdata.opendatasoft.com/page/home/.

“Progressive Cross-Validation.” visité le 22.02.2022. https://medium.com/@soumyachess1496/cross-validation-in-time-series-566ae4981ce4.

Van Buuren, Stef. 2018. Flexible Imputation of Missing Data. CRC press.