1-Introduction

Le flux des patients aux services d’urgence n’a cessé d’augmenter en France depuis quelques années, ce qui met ces services dans des situations de tensions parfois critiques et de plus en plus fréquentes. Afin d’améliorer la stratégie de planification des hôpitaux et d’allocations des ressources matérielles et humaines, une solution doit être établie. Sur ce plan, l’utilisation des séries chronologiques s’avère être un moyen très efficace pour étudier un tel phénomène, et d’autre part on cite les réseaux de neurones qui ont pu gagner une attention très répandue parmi les chercheurs de par leur puissance et aussi des nouvelles techniques qu’ils apportent pour les exercices de prédiction.

Ces dernières années, l’afflux des patients aux services d’urgence n’a fait qu’augmenter de façon remarquable, chose qui a déclenché l’alarme dans les directions hospitalières. Une saturation inhabituelle du nombre de patients admis aux services hospitaliers sur des périodes de l’année a rendu la prédiction des flux des futurs patients, très nécessaire pour bien les accueillir aux services d’urgences.

Dans un désir de comprendre ce phénomène, et afin d’y trouver une solution pour palier à l’augmentation du flux des patients au niveau des services d’urgences, qui nécessitent des ressources humaines et matérielles importantes mais limitées. Notre travail consiste, à étudier les flux des patients journaliers aux services d’urgence situé dans la région du Grand-Est et de prédire le futur flux des patients en se basant sur les données du passé.

Dans un premier lieu, on va privilégier l’utilisation des différentes méthodes statistiques se rapportant aux séries chronologiques pour l’exercice de prévision à savoir les modèles ARMA qui fournissent un outil d’aide à la décision utile et facile pour prédire la charge du travail des services d’urgences et permettront donc d’optimiser la planification de l’allocation des ressources hospitalières pour une meilleure gestion des situations de tension. Dans un second lieu, on optera la prédiction des flux de patient future via les algorithmes de réseau de neuronnes telles que les RNN, LSTM , XGBOOST …

2-Les série chronologique

Introduction au modèle ARIMA.

ARIMA signifie moyenne mobile intégrée auto-régressive et est spécifié par ces trois paramètres d’ordre: (p, d, q) . Le processus d’adaptation d’un modèle ARIMA est parfois appelé méthode de Box-Jenkins.

Un autorégressif (AR (p)) composant se réfère à l’utilisation de valeurs passées dans l’équation de régression pour la série Y . Le paramètre auto-régressif p spécifie le nombre de retards utilisés dans le modèle. Par exemple, AR (p) ou, de manière équivalente, ARIMA (2,0,0), est représenté par:

\[y_{t}= c + \phi_1 y_{t-1} +... +\phi_p y_{t-p} + \epsilon_{t} \] ou \(\phi_1\) et \(\phi_2\) sont les paramètres du modèle.

Le d représente le degré de différenciation dans la composante intégrée ( I (d) ) . Différencier une série implique simplement de soustraire ses valeurs actuelles et précédentes d fois. Souvent, la différenciation est utilisée pour stabiliser la série lorsque l’hypothèse de stationnarité n’est pas satisfaite, ce que nous verrons plus loin.

Une composante de moyenne mobile (MA (q)) représente l’erreur du modèle sous la forme d’une combinaison de termes d’erreur précédents e t . L’ordre q détermine le nombre de termes à inclure dans le modèle

\[y_{t}= c + \theta_1 \epsilon_{t-1}+...+\theta_q \epsilon_{t-q}+ \epsilon_t\]

Les composantes différentielles, autorégressives et moyennes mobiles constituent un modèle ARIMA non saisonnier qui peut être écrit sous forme d’équation linéaire:

\[y_t= \sum_{k=1}^{p} \phi_k \Delta y_{t-k} + \sum_{l=1}^{p} \theta_l \epsilon_{t-l} + \epsilon_{t} \]

Prédiction et évaluation des paramètres.

library(tseries)
library("ff","ffbase","biglm")
library(lubridate)
library(tidyverse)
library(timetk)
library(tidyquant)
library(scales)
library(forecast)   #  forecasting pkg
library(sweep)  
# library(timekit)
library(broom)
library(tibble)
library(stringr)
library(magrittr)
library(dplyr)
library(ggthemes)
load("C:/Users/User-LENOVO/Documents/Est-Rescue/ffdata/ffdac.RData")
df.tbl<- tbl_df(ffx) #convertir en base de donnée locale 
tbl.reims<-df.tbl %>% 
  select(entree,Insee_LibBV ) %>%
  filter(Insee_LibBV == "Reims")

tbl.reims<-cbind(tbl.reims , value=rep(1,dim(tbl.reims)[1])) 

tbl.reims.count<- tbl.reims%>% 
               group_by(date= as.Date(entree , format="%d/%m/%Y" ) ) %>% 
               summarise(sum =sum(value))
tbl.reims.count%>%head()
## # A tibble: 6 x 2
##   date         sum
##   <date>     <dbl>
## 1 2013-01-01   213
## 2 2013-01-02   231
## 3 2013-01-03   216
## 4 2013-01-04   227
## 5 2013-01-05   227
## 6 2013-01-06   224
dim(tbl.reims.count)
## [1] 1826    2

Un bon point de départ est de tracer la série et de l’examiner visuellement pour y détecter d’éventuelles valeurs aberrantes, de la volatilité ou des irrégularités. Dans notre cas, l’évolution du nombre de patient par jours présentent de nombreuses fluctuations d’un jour à l’autre. Cependant, même avec cette volatilité présente, nous voyons déjà des tendances se dégager.

tbl.reims.count%>%ggplot(aes(date,sum))+geom_line()+


   labs(title = "Evolution de nombre de patient ", x = "Date", y = "Nombre de \n patient par jour",
       subtitle = "de 2013 & 2018")

Dans certains cas, l’évolution du nombre de patient admet des fluctuations suspéctes (le nombre de patient tombe au-dessous de 100 et dépasse les 300 patients le jour suivant). Ce sont des valeurs suspectes qui pourraient biaiser le modèle en faussant les résumés statistiques.

Même après la suppression des valeurs aberrantes, les données quotidiennes restent assez volatiles et irrégulière.

Pour contourner ce problème on transforme la série e par l’un des concepts les plus simples - mais aussi très utiles - de l’analyse des séries chronologiques, appelé moyenne mobile. C’est un concept intuitif qui fait la moyenne des points sur plusieurs périodes, lissant ainsi les données observées en une série plus stable.

Formellement, une moyenne mobile (MA) d’ordre m peut être calculée en prenant une moyenne des séries Y , k autour de chaque point:

\[ MA=\frac{1}{m} \sum_{j=-k}^{k} y_{t-j}\]

Notons que la moyenne mobile dans ce contexte est distincte de la composante MA (q) de la définition ARIMA ci-dessus qui fait référence aux décalages d’erreur .

tbl.reims.count%>%is.na()%>%sum() #il n'ya pas de valeurs manquantes . 
## [1] 0
y=tbl.reims.count["sum"]

tbl.reims.count<-tbl.reims.count%>%mutate(clean=tsclean(as.matrix.data.frame(y)))
tbl.reims.ma=tbl.reims.count%>%mutate(clean=tbl.reims.count$clean,ma30=ma(tbl.reims.count$clean, order=30),ma90 = as.numeric(ma(tbl.reims.count$clean, order=60)))



p1<-ggplot() +
  geom_line(data =tbl.reims.ma, aes(x = date, y = clean, colour = "Actual"))+
  geom_line(data = tbl.reims.ma, aes(x = date, y = ma30,   colour = "MA~30 jours "))  +
  geom_line(data = tbl.reims.ma, aes(x = date, y = ma90, colour = "MA~90 jours "))+
  labs(title="Moyenne Mobile",y="Nbre de patients",colour="MA") 

p2<-ggplot() + geom_line(data =tbl.reims.ma, aes(x = date, y = ma30, colour = "MA~30 jours "))

p3<-ggplot() + geom_line(data = tbl.reims.ma, aes(x = date, y = ma90, colour = "90 jours "))
p4<-ggplot() + geom_line(data = tbl.reims.ma, aes(x = date, y = clean , colour = "90 jours "))

easyGgplot2::ggplot2.multiplot(p1, p4 ,p2 , p3 ,cols=2)

Décomposition de la Série.

Les éléments constitutifs d’une analyse de série chronologique sont la saisonnalité, la tendance et le cycle. Ces composants intuitifs capturent les modèles historiques de la série. Toutes les séries ne comporteront pas les trois (ou une quelconque) de ces composantes, mais si elles sont présentes, leur décomposition peut vous aider à comprendre son comportement et à jeter les bases de la construction d’un modèle de prévision.

La composante saisonnière fait référence aux fluctuations des données liées aux cycles du calendrier.

La composante tendance est la tendance générale de la série

La composante cycle consiste en des tendances décroissantes ou croissantes non saisonnières.

Enfin, une partie de la série qui ne peut pas être attribuée à des composantes saisonnières, cycliques ou de tendance est appelée résiduelle ou erreur.

count_ma = ts(na.omit(tbl.reims.ma$ma30), frequency=30)
decomp = stl(count_ma, s.window="periodic")
deseasonal_cnt <- seasadj(decomp)
plot(decomp)

Stationarité

Le modèle ARIMA exige que la série soit immobile . Une série est dite stationnaire lorsque sa moyenne, sa variance et son autocovariance sont invariantes dans le temps. Cette hypothèse est logique: puisque ARIMA utilise des décalages de séries antérieurs pour modéliser son comportement, la modélisation de séries stables à propriétés cohérentes implique moins d’incertitude.

Le test Dickey-Fuller augmenté (ADF) est un test statistique formel pour la stationnarité. L’hypothèse nulle suppose que la série est non stationnaire. La procédure ADF teste si le changement de Y peut être expliqué par une valeur décalée et une tendance linéaire. Si la contribution de la valeur décalée à la variation de Y est non significative et qu’il existe une composante de tendance, la série est non stationnaire et l’hypothèse nulle ne sera pas rejetée.

adf.test(count_ma, alternative = "stationary") #p_value=0.01 ==> serie stationnaire. #d=0
## Warning in adf.test(count_ma, alternative = "stationary"): p-value smaller
## than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  count_ma
## Dickey-Fuller = -6.871, Lag order = 12, p-value = 0.01
## alternative hypothesis: stationary
deseasonal_cnt  %>% ggtsdisplay()

Il existe des auto-corrélations significatives aux décalages 0, 1. Les diagrammes de corrélation partielle montrent un pic important aux décalages 3 ,5 et 7. Cela suggère que nous pourrions vouloir tester des modèles avec des composantes AR ou MA d’ordre 1, 2,3,5 ou 7.

Choix du modèle & Estimation des paramètres

Il existe un certain nombre de critères de ce type permettant de comparer la qualité de l’ajustement entre plusieurs modèles. Les deux critères les plus utilisés sont les critères d’information d’Akaike (AIC) et les critères d’information de Baysian (BIC). Ces critères sont étroitement liés et peuvent être interprétés comme une estimation de la quantité d’informations perdues si un modèle donné est choisi. En comparant les modèles, on veut minimiser les AIC et BIC. la fonction auto.arima séléctionne le meilleur modèle satisfaisant les critères de validation de modèle.

auto.arima(deseasonal_cnt, seasonal=FALSE) #ARIMA(5,1,3)
## Series: deseasonal_cnt 
## ARIMA(5,1,3) 
## 
## Coefficients:
##           ar1      ar2     ar3     ar4     ar5     ma1     ma2     ma3
##       -1.2580  -0.3135  0.6311  0.4248  0.2315  2.5969  2.5559  0.9535
## s.e.   0.0265   0.0462  0.0363  0.0369  0.0332  0.0124  0.0252  0.0330
## 
## sigma^2 estimated as 0.2461:  log likelihood=-1287.44
## AIC=2592.87   AICc=2592.97   BIC=2642.31
fit<-auto.arima(deseasonal_cnt, seasonal=FALSE)
# tsdisplay(residuals(fit), lag.max=45, main='(5,1,3) Model Residuals')

Le modèle séléctioné est le modèle ARIMA(5,1,3) : \[\Delta y_t= -1.26\Delta y_{t-1} -0.31 \Delta y_{t-2} + 0.63 \Delta y_{t-3} + 0.42 \Delta y_{t-4} + 0.23 \Delta y_{t-5} + 2.6\epsilon_{t-1} + 2.6\epsilon_{t-2} + 0.95\epsilon_{t-3} + \mu_t \]

Verification des Résidus

Pouvons-nous faire confiance à ce modèle? Nous pouvons commencer par examiner les graphiques ACF et PACF pour déterminer les résidus du modèle. Si les paramètres et la structure de l’ordre du modèle sont correctement spécifiés, nous ne nous attendons à aucune autocorrélation significative.

checkresiduals(fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(5,1,3)
## Q* = 1016.9, df = 52, p-value < 2.2e-16
## 
## Model df: 8.   Total lags used: 60

Les résidus de ce modèle sont illustrés à la figure ci dessous. Il existe quelques pics significatifs dans ACF et le modèle échoue au test de Ljung-Box. Le modèle peut toujours être utilisé pour la prévision, mais les intervalles de prévision peuvent ne pas être précis en raison des résidus corrélés.

Prévision

autoplot(forecast(fit) , h=90)

Correction du modèle

Les modèles ARIMA peuvent être ajustés à la fois aux données saisonnières et non saisonnières. ARIMA saisonnier nécessite une spécification plus complexe de la structure du modèle, bien que le processus de détermination (P, D, Q) soit similaire à celui consistant à choisir des paramètres d’ordre non saisonniers.

#Forcast ARIMA(p,d,q)(P,D,Q)= (1,1,2)(2,0,1)[30] 
fit_w_seasonality = auto.arima(deseasonal_cnt, seasonal=TRUE)
fit_w_seasonality
## Series: deseasonal_cnt 
## ARIMA(1,1,2)(2,0,1)[30] 
## 
## Coefficients:
##          ar1     ma1      ma2     sar1     sar2     sma1
##       0.9526  0.1991  -0.7311  -0.1584  -0.1036  -0.6839
## s.e.  0.0102  0.0219   0.0212   0.0323   0.0288   0.0240
## 
## sigma^2 estimated as 0.1495:  log likelihood=-854.22
## AIC=1722.44   AICc=1722.5   BIC=1760.88
checkresiduals(fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(5,1,3)
## Q* = 1016.9, df = 52, p-value < 2.2e-16
## 
## Model df: 8.   Total lags used: 60
seas_fcast <- forecast(fit_w_seasonality, h=30)
autoplot(forecast(seas_fcast) , h=90)

Validation du modèle: hold_up cross validation

hold0 <- window(ts(deseasonal_cnt), start=1700)

fit_no_holdout0 = arima(ts(deseasonal_cnt[-c(1701:1796)]), order=c(5,1,3))

fcast_no_holdout0 <- forecast(fit_no_holdout0,h=96)

plot(fcast_no_holdout0, main="Modèle ARIMA(5,1,3) ")
lines(ts(deseasonal_cnt))

hold1 <- window(ts(deseasonal_cnt), start=1700)
fit_no_holdout1 = arima(ts(deseasonal_cnt[-c(1701:1796)]), order=c(1,1,2) , seasonal = c(2,0,1))
fcast_no_holdout1 <- forecast(fit_no_holdout1,h=96)
plot(fcast_no_holdout1, main="Modèle ARIMA(1,1,2)(2,0,1) ")
lines(ts(deseasonal_cnt))

taux.erreur0<- mean( (fcast_no_holdout0$mean - hold0)^2   )
taux.erreur1<- mean( (fcast_no_holdout1$mean - hold1)^2  )

modele<-c("ARIMA(5,1,3)" , "ARIMA(1,1,2)(2,0,1)")
erreur<-c(paste(round(taux.erreur0,3),"%") , paste(round(taux.erreur1,3) ,"%") ) 

tab<-cbind(modele,erreur)
knitr::kable(tab)
modele erreur
ARIMA(5,1,3) 34.158 %
ARIMA(1,1,2)(2,0,1) 35.973 %