Prédiction des retards à la SNCF - PARIS-REIMS
Le but de cet article est de prédire la régularité mensuelle du TGV Paris Reims pour 2018. L’article a deux buts :
– Comparer la régularité du Paris-Reims avec les autres trains pour le situer par rapport aux autres – Prédire la régularité du Paris-Reims sur les mois à venir
Pour faire court, la qualité de service offerte par SNCF est plutôt exceptionnelle avec très peu de retard rapportée à la quantité énorme de voyageurs.
1 Chargement des données
Dans sa politique d’OPEN DATA, l’Etat français et mécaniquement la SNCF ont mis en place un libre accès aux informations publiques dans leur souci de transparence, d’efficience et de communication.
Sur TGV, la régularité est calculée lors de l’arrivée du train à la dernière gare de son parcours (terminus). Ce mode de calcul, qui n’intègre pas les gares intermédiaires, propose le retard cumulé sur l’ensemble d’un trajet.
L’indicateur utilisé est celui de la « régularité composite » : un train est considéré en retard s’il arrive cinq minutes après son heure prévue pour un trajet de moins de 1h30, dix minutes pour un trajet de 1h30 à 3h, et quinze minutes pour un trajet de plus de 3h.
Un train supprimé avant 16h la veille de sa circulation ne sera pas comptabilisé. En revanche, si l’annonce de sa suppression n’a pu être faite la veille avant 16h, le train sera bien comptabilisé comme un train supprimé, que cette suppression soit totale ou partielle (le train a effectué qu’une partie de son parcours).
https://ressources.data.sncf.com/explore/?sort=modified
Cette première phase porte sur le data cleaning afin de préparer les données pour les différentes données.
setwd("~/Desktop/projet")
data<-read.csv("regularite_tgv.csv", sep = ";")
colnames(data)<-c("date","axe","origins","destinations","theorique","effectif","annule","nb_retard","regularite","commentaires")
data$origins<-as.character(data$origins)
data$destinations<-as.character(data$destinations)library(anytime)## Warning: package 'anytime' was built under R version 3.3.2data$date <- anydate(data$date)
data<-data[order(as.Date(data$date)),]1.1 Portrait robot des TGV aux départs de Paris
library(dplyr)## Warning: package 'dplyr' was built under R version 3.3.2##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unioncompa<-subset(data, origins == "PARIS EST" | origins == "PARIS MONTPARNASSE" | origins == "PARIS LYON" | origins == "PARIS NORD")
compa<-compa[,c(1,4:9)]
compa <- compa %>%
group_by(destinations) %>%
summarize(
theorique=sum(theorique),
effectif=sum(effectif),
annule=sum(annule))
compa<-compa[order(as.numeric(compa$annule)),]
compa<-compa[-c(36),]Voici le nombre de trains annulés depuis 2011 aux départs de paris (toutes les gares)
library(ggplot2)## Warning: package 'ggplot2' was built under R version 3.3.2p<- ggplot(data=compa, aes(x=destinations, y=annule)) +
geom_bar(stat="identity", fill="steelblue") +
theme_minimal()
p<- p + coord_flip()
p
compa$failure<-round(compa$annule/compa$theorique*100,2)
p<- ggplot(data=compa, aes(x=destinations, y=failure)) +
geom_bar(stat="identity", fill="red")+
theme_minimal()
p<- p + coord_flip()
p## Warning: Removed 1 rows containing missing values (position_stack).
2 Train au départ de Paris
Commençons par chercher le train qui ressemble le plus à mon PARIS-REIMS
2.1 Quels sont les TGVs au départ de Paris Gare de l’EST ?
library(dplyr)
choice<-subset(data, origins == "PARIS EST")
df<-choice[,c(3,4)]
df2 <- df %>%
group_by(origins, destinations) %>%
summarize(counts = n()) %>%
ungroup() %>%
arrange(desc(counts))
library(networkD3)## Warning: package 'networkD3' was built under R version 3.3.2name_vec <- c(unique(df2$origins), unique(df2$destinations))
n<-NROW(name_vec)
n<-n-1
nodes <- data.frame(name = name_vec, id = 0:n)
links <- df2 %>%
left_join(nodes, by = c('origins' = 'name')) %>%
rename(origin_id = id) %>%
left_join(nodes, by = c('destinations' = 'name')) %>%
rename(dest_id = id)## Warning: Column `origins`/`name` joining character vector and factor,
## coercing into character vector## Warning: Column `destinations`/`name` joining character vector and factor,
## coercing into character vectorsankeyNetwork(Links = links, Nodes = nodes, Source = 'origin_id', Target = 'dest_id',
Value = 'counts', NodeID = 'name', fontSize = 16)## Links is a tbl_df. Converting to a plain data frame.2.2 Quels sont les trains partant des autres gares parisiennes?
library(dplyr)
choice<-subset(data, origins == "PARIS NORD" | origins == "PARIS MONTPARNASSE" | origins == "PARIS LYON")
df<-choice[,c(3,4)]
df2 <- df %>%
group_by(origins, destinations) %>%
summarize(counts = n()) %>%
ungroup() %>%
arrange(desc(counts))
library(networkD3)
name_vec <- c(unique(df2$origins), unique(df2$destinations))
n<-NROW(name_vec)
n<-n-1
nodes <- data.frame(name = name_vec, id = 0:n)
links <- df2 %>%
left_join(nodes, by = c('origins' = 'name')) %>%
rename(origin_id = id) %>%
left_join(nodes, by = c('destinations' = 'name')) %>%
rename(dest_id = id)## Warning: Column `origins`/`name` joining character vector and factor,
## coercing into character vector## Warning: Column `destinations`/`name` joining character vector and factor,
## coercing into character vectorsankeyNetwork(Links = links, Nodes = nodes, Source = 'origin_id', Target = 'dest_id',
Value = 'counts', NodeID = 'name')## Links is a tbl_df. Converting to a plain data frame.2.3 Quels sont les TGVs qui ressemblent le plus aux caractéristiques de mon TGV Paris-Reims ?
Rappelons que le dataset contient 4 variables quantitatives :
- Un nombre de trains théorique : nombre de trains programmé / qui devrait partir
- Un nombre de trains effectif : nombre de trains qui arrive à bon quai / qui sont réellement partis
- Un nombre de train en retard
- Un % de régularité
compa2<-compa[,c(2:4)]
rownames(compa2)<-compa$destinations## Warning: Setting row names on a tibble is deprecated.# Compute dissimilarity matrix
d <- dist(compa2, method = "euclidean")
# Hierarchical clustering using Ward's method
res.hc <- hclust(d, method = "ward.D2" )
# Cut tree into 4 groups
grp <- cutree(res.hc, k = 4)
# Visualize
plot(res.hc, cex = 0.6) # plot tree
2.3.1 Clustering
## Warning: package 'factoextra' was built under R version 3.3.2## Welcome! Related Books: `Practical Guide To Cluster Analysis in R` at https://goo.gl/13EFCZ## Warning: Setting row names on a tibble is deprecated.
## Warning: argument frame is deprecated; please use ellipse instead.## Warning: argument frame.type is deprecated; please use ellipse.type
## instead.
2.4 Choix de ma destination de prédilection ? Reims
Ici, l’échantillon est réduit au train qui ont pour destination REIMS au départ de Paris EST.
df<-subset(data, origins =="REIMS")3 Train prévu vs. Train à l’heure
Ce graphe décrit la différence entre les trains programmés (qui devraient partir théoriquement) et les trains effectifs (qui arriveront effectivement à Reims)
library(ggplot2)
ggplot(df, aes(date)) +
geom_point(aes(y = theorique,colour = "Nombre de trains théoriques")) +
geom_line(aes(y = effectif, colour = "Nombre de trains effectifs"))
4 Etude d’une série temporelle : Prévision des trains qui arriveront à l’heure
Un modèle ARIMA cherche à déterminer chaque valeur de la série en fonction des valeurs qui la précède yt = f(yt-1, yt-2, …)
— Les processus autorégressifs supposent que chaque point peut être prédit par la somme pondérée d’un ensemble de points précédents, plus un terme aléatoire d’erreur.
— Le processus d’intégration suppose que chaque point présente une différence constante avec le point précédent.
— Les processus de moyenne mobile supposent que chaque point est fonction des erreurs entachant les points précédant, plus sa propre erreur.
Un modèle ARIMA est étiqueté comme modèle ARIMA (p,d,q), dans lequel: - p est le nombre de termes auto-régressifs - d est le nombre de différences - q est le nombre de moyennes mobiles.
df_effectif<-df$effectif
plot(df_effectif, main = "Nombre de Trains arrivés à l'heure",t = "l", col = "blue", xlab = "temps", ylab = "Nombre de trains")
Calcul de l’écart type :
sd(df_effectif)## [1] 26.36837L’estimation des modèles ARIMA suppose que l’on travaille sur une série stationnaire. Ceci signifie que la moyenne de la série est constante dans le temps, ainsi que la variance. La meilleure méthode pour éliminer toute tendance est de différencier, c’est-à-dire de remplacer la série originale par la série des différences adjacentes. Une série temporelle qui a besoin d’être différenciée pour atteindre la stationnarité est considérée comme une version intégrée d’une série stationnaire (d’où le terme Integrated).
data_diff1 <- diff(df_effectif, differences=1)
plot.ts(data_diff1, t = "l",col = "red")
sd(data_diff1)## [1] 20.57428acf(df_effectif)
Pourquoi la différentiation semble utile ?
- Ecart type significativement réduit on passe de 26.36837 à 20.57428
- Le graphique de la fonction d’auto-corrélation présente une régression exponentielle qui tend vers 0
Box.test(df_effectif, lag = 1, type = c("Box-Pierce", "Ljung-Box"), fitdf = 0)##
## Box-Pierce test
##
## data: df_effectif
## X-squared = 33.806, df = 1, p-value = 6.088e-09Compute the Box–Pierce or Ljung–Box test statistic for examining the null hypothesis of independence in a given time series. These are sometimes known as ‘portmanteau’ tests.
If p-value < 0.05: You can reject the null hypothesis assuming a 5% chance of making a mistake. So you can assume that your values are showing dependence on each other.
If p-value > 0.05: You don’t have enough statistical evidence to reject the null hypothesis. So you can not assume that your values are dependent. This could mean that your values are dependent anyway or it can mean that your values are independent. But you are not proving any specific possibility, what your test actually said is that you can not assert the dependence of the values, neither can you assert the independence of the values.
In general, what is important here is to keep in mind that p-value < 0.05 lets you reject of the null-hypothesis, but a p-value > 0.05 does not let you confirm the null-hypothesis.
acf(data_diff1)
pacf(data_diff1) 
The partial correlogram shows that the partial autocorrelations at lags 1, 3 and 13 exceed the significance bounds, are negative. The partial autocorrelations tail off to zero after lag 3.
4.1 TROUVER AR ? Auto-corrélation partielle
Une extinction brutale de l’auto-corrélation partielle associée à un déclin plus progressif de l’auto-corrélation constitue la signature caractéristique d’un processus auto-régressif. — Plus particulièrement, l’auto-corrélation partielle de décalage k est égale au coefficient AR(k) estimé dans un modèle contenant k termes AR. On pourrait d’ailleurs déterminer les coefficients AR par régression multiple, en prédisant (yt-yt-1) à partir de k échantillons représentant les k décalages. Le décalage auquel l’auto-corrélation partielle disparaît indique le nombre de termes autorégressifs à inclure. — Généralement, ce pattern est associé à une auto-corrélation de décalage 1 positive, signe que la série demeure sous-différenciée. Une légère sous-différenciation peut donc être compensée par l’ajout d’un terme auto-régressif.
4.2 TROUVER MA ? Auto-corrélation
La fonction d’auto-corrélation joue pour les processus de moyenne mobile le même rôle que la fonction d’auto-corrélation partielle pour les processus auto-régressifs.
Si l’auto-corrélation est significative au décalage k mais plus au décalage k+1, ceci indique que k termes de moyenne mobile doivent être ajoutés au modèle.
Une signature MA est généralement associée à une auto-corrélation négative au décalage 1, signe que la série est sur-différenciée. Une légère sur-différenciation peut donc être compensée par l’ajout d’un terme de moyenne mobile.
4.3 TROUVER LE BON MODELE ? # (IMPORTANT CAR conclusions)
- MA s’occupe de la sur-différentiation
- AR s’occupe de la sous-différentiation
Du coup, si AR alors on différencie au niveau 1 => sous-différencié et on corrige avec des termes AR Au contraire, si MA alors on différencie un peu plus (d’ordre 2) et on corrige en ajoutant des termes MA
Conclusions : deux modèles alternatifs soient possibles: un premier avec 0 ou 1 ordre de différenciation combiné avec des termes AR, et un autre avec le niveau de différenciation supérieur, combiné à des termes MA.
4.4 Règle from DUKE
If the partial autocorrelation function (PACF) of the differenced series displays a sharp cutoff and/or the lag-1 autocorrelation is positive–i.e., if the series appears slightly “underdifferenced”–then consider adding one or more AR terms to the model. The lag beyond which the PACF cuts off is the indicated number of AR terms.
If the autocorrelation function (ACF) of the differenced series displays a sharp cutoff and/or the lag-1 autocorrelation is negative–i.e., if the series appears slightly “overdifferenced”–then consider adding an MA term to the model. The lag beyond which the ACF cuts off is the indicated number of MA terms.
library(forecast)## Warning: package 'forecast' was built under R version 3.3.2arima(df_effectif, c(4, 1, 0))##
## Call:
## arima(x = df_effectif, order = c(4, 1, 0))
##
## Coefficients:
## ar1 ar2 ar3 ar4
## -0.2916 -0.2479 -0.0749 -0.5676
## s.e. 0.0997 0.1061 0.1049 0.1077
##
## sigma^2 estimated as 257: log likelihood = -294.45, aic = 598.89arima(df_effectif, c(0, 1, 4))##
## Call:
## arima(x = df_effectif, order = c(0, 1, 4))
##
## Coefficients:
## ma1 ma2 ma3 ma4
## -0.2854 -0.0846 0.0861 -0.4765
## s.e. 0.1004 0.1189 0.0980 0.1109
##
## sigma^2 estimated as 286.3: log likelihood = -297.96, aic = 605.92arima(df_effectif, c(1, 1, 2))##
## Call:
## arima(x = df_effectif, order = c(1, 1, 2))
##
## Coefficients:
## ar1 ma1 ma2
## -0.5150 0.1412 -0.5707
## s.e. 0.2108 0.2236 0.1352
##
## sigma^2 estimated as 329.9: log likelihood = -302.64, aic = 613.27arima(df_effectif, c(4, 1, 4))##
## Call:
## arima(x = df_effectif, order = c(4, 1, 4))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 ma2 ma3 ma4
## 0.0377 -0.1585 0.0848 -0.8716 -0.4001 0.0679 -0.2214 0.5526
## s.e. 0.1242 0.1182 0.1283 0.1088 0.1781 0.1685 0.2218 0.2400
##
## sigma^2 estimated as 239.2: log likelihood = -292.89, aic = 603.77| Criteres | ARIMA(4, 1, 0) | ARIMA(0, 1, 4) | ARIMA(1, 1, 2) | ARIMA(1, 1, 2) | ARIMA(4,1,4) | |
|---|
| AIC | 598.89 | 605.92 | 613.27 | 613.27 | 603.77 | |
Nous allons maintenant tester les performances des 3 modèles. Pour ce faire, nous enlevons les 10 derniers points de notre série Yt, pour ensuite comparer nos 10 valeurs prédites avec nos 10 valeurs réèlles.
fit <- Arima(df_effectif, order=c(4,1,0))
summary(fit)## Series: df_effectif
## ARIMA(4,1,0)
##
## Coefficients:
## ar1 ar2 ar3 ar4
## -0.2916 -0.2479 -0.0749 -0.5676
## s.e. 0.0997 0.1061 0.1049 0.1077
##
## sigma^2 estimated as 272.6: log likelihood=-294.45
## AIC=598.89 AICc=599.83 BIC=610.14
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -1.576101 15.91729 11.8052 -1.465414 6.383638 0.8607959
## ACF1
## Training set -0.077155Acf(residuals(fit))
Box.test(residuals(fit), lag=24, fitdf=4, type="Ljung")##
## Box-Ljung test
##
## data: residuals(fit)
## X-squared = 15.017, df = 20, p-value = 0.7754plot(forecast(fit))
Démarche entreprise :
— 1. Plot the data. Identify any unusual observations. — 2. If necessary, transform the data (using a Box-Cox transformation) to stabilize the variance. — 3. If the data are non-stationary: take first differences of the data until the data are stationary. — 4. Examine the ACF/PACF: Is an AR(pp) or MA(qq) model appropriate? — 5. Try your chosen model(s), and use the AICc to search for a better model. — 6. Check the residuals from your chosen model by plotting the ACF of the residuals, and doing a portmanteau test of the residuals. If they do not look like white noise, try a modified model. — 7. Once the residuals look like white noise, calculate forecasts.
For step 6 : The ACF plot of the residuals from the ARIMA(3,1,1) model shows all correlations within the threshold limits indicating that the residuals are behaving like white noise. A portmanteau test returns a large p-value, also suggesting the residuals are white noise.
Acf(residuals(fit))
Box.test(residuals(fit), lag=24, fitdf=4, type="Ljung")##
## Box-Ljung test
##
## data: residuals(fit)
## X-squared = 15.017, df = 20, p-value = 0.77544.4.1 Dynamic regression models
library(forecast)
fit1 <- Arima(df_effectif, order=c(4,1,0), include.drift=TRUE)
fit2 <- Arima(df_effectif, order=c(4,1,1), include.drift=TRUE)
par(mfrow=c(2,1))
plot(forecast(fit2), ylim=c(120,250))
plot(forecast(fit1), ylim=c(120,250))
Voici le nombre de train qui devrait réellement arrivé dans les mois à venir :
library(forecast)
forecast(fit2)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 72 146.4609 125.5146 167.4072 114.42633 178.4955
## 73 144.1621 120.0799 168.2444 107.33154 180.9927
## 74 159.2747 133.6806 184.8689 120.13184 198.4176
## 75 145.7261 118.4964 172.9557 104.08189 187.3702
## 76 181.5910 154.1873 208.9946 139.68070 223.5013
## 77 182.0352 154.4694 209.6009 139.87691 224.1934
## 78 165.4263 137.0142 193.8384 121.97371 208.8789
## 79 170.9900 141.9612 200.0188 126.59430 215.3858
## 80 149.0520 117.2773 180.8268 100.45672 197.6473
## 81 147.3380 114.3419 180.3341 96.87487 197.80125 Mise en place d’un modèle à Réseau de Neurones
fit <- nnetar(df_effectif)
plot(forecast(fit,h=20))
6 Prediction du nombre de trains annulés pour 2018 sur Paris-Reims
library(plotly)## Warning: package 'plotly' was built under R version 3.3.2##
## Attaching package: 'plotly'## The following object is masked from 'package:ggplot2':
##
## last_plot## The following object is masked from 'package:stats':
##
## filter## The following object is masked from 'package:graphics':
##
## layoutp <- plot_ly(df, x = ~date, y = ~annule, name = 'Trains Annulés', type = 'scatter', mode = 'lines') %>%
add_trace(y = ~nb_retard, name = 'Trains en retard', mode = 'lines+markers')
psd(df$annule)## [1] 2.11643df_annule_diff <- diff(df$annule, differences=1)
plot.ts(df_annule_diff, t = "l",col = "green")
sd(df_annule_diff)## [1] 3.102126acf(df$annule)
acf(df$annule)
pacf(df$annule) 
fit <- Arima(df$annule, order=c(0,0,0))
summary(fit)## Series: df$annule
## ARIMA(0,0,0) with non-zero mean
##
## Coefficients:
## mean
## 0.6761
## s.e. 0.2494
##
## sigma^2 estimated as 4.479: log likelihood=-153.47
## AIC=310.94 AICc=311.12 BIC=315.47
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -9.992007e-14 2.101473 0.952192 -Inf Inf 0.8545313
## ACF1
## Training set -0.06030234Acf(residuals(fit))
Box.test(residuals(fit), lag=24, fitdf=4, type="Ljung")##
## Box-Ljung test
##
## data: residuals(fit)
## X-squared = 18.388, df = 20, p-value = 0.5618plot(forecast(fit))
Voici le nombre de train qui devrait être annulé dans les mois à venir sur le trajet PARIS-REIMS :
Soit moins d’un train par mois devrait être en retard dans les mois à venir
forecast(fit)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 72 0.6760563 -2.036258 3.38837 -3.47207 4.824183
## 73 0.6760563 -2.036258 3.38837 -3.47207 4.824183
## 74 0.6760563 -2.036258 3.38837 -3.47207 4.824183
## 75 0.6760563 -2.036258 3.38837 -3.47207 4.824183
## 76 0.6760563 -2.036258 3.38837 -3.47207 4.824183
## 77 0.6760563 -2.036258 3.38837 -3.47207 4.824183
## 78 0.6760563 -2.036258 3.38837 -3.47207 4.824183
## 79 0.6760563 -2.036258 3.38837 -3.47207 4.824183
## 80 0.6760563 -2.036258 3.38837 -3.47207 4.824183
## 81 0.6760563 -2.036258 3.38837 -3.47207 4.8241836.1 Mise en place d’un modèle à Réseau de Neurones
Le but est de prédire le nombre de trains annulés à venir
fit <- nnetar(df$annule)
plot(forecast(fit,h=20))
Conclusion : Mon train Paris-Reims a très peu de chance d’être en retard. La SNCF offre une qualité de service plutôt exceptionnelle.
Voici le nombre de trains qui devraient prendre le départ de PARIS GARDE DE L’EST en direction de REIMS
knitr::include_graphics("image1.png")