Les Suisses et le diabète
Alain Simon
2023-06-23
Information
L’objectif de ce projet sera d’établir les liens existants entre les variables et de prédire le nombre d’enregistrements (en ligne ou sur papier) de dossiers entre \(08\) et \(22h\) de patients diabetiques Suisses. Le jeu de données (lien ici) rennomé Diabet_Swiss, compte 943 lignes et 4 variables qui sont les suivantes:
VARIABLES QUALITATIVES:
V1= jours d’enregistrements des dossiers de patients diabétiques suisses. ( du “1991-04-21” au “1991-09-03” ( format année-mois-jour ))
V2= les horaires d’enregistrements: sur papier (\(8h30\),\(12h00\),\(18h00\),\(22h00\)) et en ligne (les autres)
V3= Codes des raisons pour lesquelles ont été faites les enrégistrements (variable catégorielle):
- 33 = Dose d'insuline ordinaire
- 34 = Dose d'insuline NPH
- 48 = Mesure de la glycémie non spécifiée
- 58 = Mesure de la glycémie avant le petit-déjeuner
- 60 = Mesure de la glycémie avant le déjeuner
- 62 = Mesure de la glycémie avant le souper
- 65 = Symptômes d'hypoglycémie
VARIABLE QUANTITATIVE:
- V4 le nombres de patients diabétiques suisses enrégistrés par jour.
Chargement des packages nécessaires.
library(ggplot2)
library(kableExtra)
library(tidyverse)
library(lubridate)
library(forcats)
library(fpp)
library(psych)
library(ggthemes)
library(dataedu)
library(dplyr)
library(extrafont)
library(tidytuesdayR)
library(ggtext)
library(forecast)
library(hrbrthemes)
library(systemfonts)
library(systemfonts)
library(ggThemeAssist)
library(psych)Un bref apperçu de notre jeu de données Diabet_Swiss.
## V1 V2 V3 V4
## 1 04-21-1991 9:09 58 100
## 2 04-21-1991 9:09 33 9
## 3 04-21-1991 9:09 34 13
## 4 04-21-1991 17:08 62 119
## ... <NA> <NA> ... ...
## 940 09-02-1991 23:00 48 155
## 941 09-03-1991 7:20 58 110
## 942 09-03-1991 7:20 33 9
## 943 09-03-1991 7:20 34 16Dans un premier temps nous allons renommer ces variables pour qu’elles soient plus parlantes, puis changer le format des dates. Voici un apperçu du tableau obtenu, cela correspond aux \(4\) premières et \(4\) dernères lignes de notre tableau de données Diabet_Swiss.
## dates horaires Codes nbre_enregist
## 1 1991-04-21 9:09 58 100
## 2 1991-04-21 9:09 33 9
## 3 1991-04-21 9:09 34 13
## 4 1991-04-21 17:08 62 119
## ... <NA> <NA> <NA> ...
## 940 1991-09-02 23:00 48 155
## 941 1991-09-03 7:20 58 110
## 942 1991-09-03 7:20 33 9
## 943 1991-09-03 7:20 34 16Regardons la struture de notre jeu de données.
## 'data.frame': 943 obs. of 4 variables:
## $ dates : Date, format: "1991-04-21" "1991-04-21" ...
## $ horaires : chr "9:09" "9:09" "9:09" "17:08" ...
## $ Codes : Factor w/ 7 levels "33","34","48",..: 4 1 2 6 1 3 4 1 2 1 ...
## $ nbre_enregist: int 100 9 13 119 7 123 216 10 13 2 ...Les variables ont bien été transformées pour nous faciliter les analyses.
Nous allons ensuite effectuer une analyse descriptive pour mieux comprendre nos données.
Analyse descriptive
## X1
## vars 1.00
## n 943.00
## mean 67.41
## sd 85.62
## median 15.00
## trimmed 52.34
## mad 16.31
## min 0.00
## max 343.00
## range 343.00
## skew 1.19
## kurtosis 0.18
## se 2.79
## IQR 115.00## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 7.00 15.00 67.41 122.00 343.00Cette sortie nous donne plusieurs statistiques, informations sur la variable nbre_enregist (nombre d’enregistrements des dossiers de diabétiques suisses), mais nous allons nous intéresser particulièrement à quelques une, la moyenne \(mean\) , \(sd\) l’écart-type, \(median\) la médiane, \(min\) le minimum, \(max\) le maximum.
Interpretation de median, mean, sd, min, max
Interpretation de la médiane de la variable nbre_enregist
La médiane est de \(15\), ce qui voudrait dire que la moitié des observations de la variable nbre_enregist, prend des valeurs supérieurs à 15. En d’autres termes, sur la moitié des horaires d’enregistements des dossiers (en ligne et sur papier ) observés de patients diabétiques suisses,on compte au moins \(15\) dossiers d’enregistrements pour chacune des ces horaires.
Ecart-type valeurs abérrantes
L’ecart-type est de \(85.62\), ce qui est quand même relativement élevé comme valeur, cela pourrait nous laisser comprendre que les observations de notre variables d’intéret nbre_enregist sont assez dispersés. Pour en etre sûr, nous allons commencer par effectuer la visualisation d’un Boxplot.
Nous avons bien confirmation de la disperssion des observations de notre variable d’intérêt.Nous remarquons aussi que les valeurs supérieures à environ \(300\) sont considérées ici comme extrèmes / abérrantes. Elles correspondent aux observations dont les valeurs sont soient supérieures à la valeur du 3ème quartile plus \(1.5\) fois l’intervalle interquartile ( \(Q_3 - Q_1\)),soit inférieures à la valeur du 1er quartile moins \(1.5\) fois l’intervalle interquartile. Elles correspondent aux valeurs n’appartenant pas à l’intervalle \([-165.5 , 294.5]\). Nous allons les appeler ici outlier_vals. Les voici:
outlier_vals <- boxplot.stats(Diabet_Swiss_bis$nbre_enregist)$out
outlier_vals## [1] 340 305 343 300 297 312 295 313 306 303 335 309récupérons les indexs de nos outliers_vals
outlier_idx <- which(Diabet_Swiss_bis$nbre_enregist %in% outlier_vals)
outlier_idx## [1] 23 84 135 155 214 235 242 278 504 508 789 834Grâce à ces indexs, nous obtenons les lignes suivantes constituant un nouveau tableau de nos valeurs abérrantes nommé Data_outlier, comme suit:
Data_outlier <-data.frame(Diabet_Swiss_bis[c(outlier_idx),])
Data_outlier## dates horaires Codes nbre_enregist
## 23 1991-04-24 22:09 48 340
## 84 1991-05-03 7:48 58 305
## 135 1991-05-10 17:30 62 343
## 155 1991-05-13 17:00 62 300
## 214 1991-05-23 7:35 58 297
## 235 1991-05-26 9:17 58 312
## 242 1991-05-27 8:50 58 295
## 278 1991-06-01 9:20 58 313
## 504 1991-06-30 12:06 60 306
## 508 1991-06-30 22:38 48 303
## 789 1991-08-09 7:38 58 335
## 834 1991-08-16 22:20 48 309Il est constitué de 12 lignes et de nos 4 variables.
Nous allons analyser les relations entre notre variable d’intérêt nbre_enregist et chacune des variables date et Codes.
Lien entre nombre d’enregistrements et dates
visualisons le tracer du nombre d’enregistrements en fonction des dates
Cette courbe des nombres d’enregistrements de la mesure du taux de glycémie en fonction la variables dates n’est ni monotone, ni constante, par conséquence, le nombre d’enregistrements n’est donc pas ou est faiblement correlé au temps.
Lien entre nombre d’enregistrements et les codes
Commençons d’abord par regarder toutes les modalités de la variable Codes de notre jeu de données Data_outlier présentes:
##
## 33 34 48 58 60 62 65
## 0 0 3 6 1 2 0Nous constatons que les seules modalités de la variables Codes présentes ici dans notre jeu de données Data_outlier sont:
48 = Mesure de la glycémie non spécifiée
58 = Mesure de la glycémie avant le petit-déjeuner
60 = Mesure de la glycémie avant le déjeuner
62 = Mesure de la glycémie avant le souper
Toutes liées à la mesure de la glycémie,
Camembert des modalités de la variable Codes
ggplot(data = Data_outlier) +
geom_bar(mapping = aes(x = 1, fill = Codes), position = "fill") +
coord_polar(theta = "y")
La modalité 58 associée à la mesure de la glycémie avant le petit déjeuner,représente 50 % de notre camembert . Ce qui voudrait dire ici que la moitié du nombre d’enregistrements de dossiers est associée aux patients diabétiques suisses qui prennent leur taux de glycémie avant le petit déjeuner.
Visualisation du boxplot du nombre d’enregistrements en fonction de la variable Codes
if(!require("ggstatsplot")){install.packages("ggstatsplot")}
library(ggstatsplot)
ggbetweenstats(data = Data_outlier, y = nbre_enregist,x=Codes)
Nous avons l’impression que les médianes 2 premiers groupes sont pratiquement égales et diffèrent des 2 autres.Effectuons un test de comparaison de variance pour en convaincre. Mais avant tout, vérifions les conditions de faisabilité
Conditions de faisabilité et tests d’hypothèses
Nous allons commencer par tester la normalité des groupes.
- la normalité des groupes :
test de Shapiro-Wlik:
Les échantillons des groupes ayant pour code 60 et 62 étant de tailles respectivement égales à \(1\) et \(2\) et donc pas compris entre \(3\) et \(5000\), nous ne pourrons effectuer de test de normalité sur eux. En effet un échantillon de taille 1 ou 2 ne fournit pas suffisamment d’observations pour tirer des conclusions sur la distribution des données. Le test ne se fera finalement uniquement que sur les échantillons ayant pour code 48et 58.
Nous avons pour hypothèses de test:
\(H_0\): Normalité des groupes
\(H_1\): Non-normalité des groupes.
##
## Shapiro-Wilk normality test
##
## data: .
## W = 0.86792, p-value = 0.2896##
## Shapiro-Wilk normality test
##
## data: .
## W = 0.90039, p-value = 0.3763Le résultat du test nous fournit des \(p-values\) supérieures au seuil de significativité de \(0.05\).
La conditions sur la normalité des groupes n’est pas vérifiée puisqu’il y’a des groupes sur lesquelles il nous est impossible d’effectuer un test de normalité à cause de la taille de chacune de leur échantillion inférieure à 3. Les données ne suivent une loi normale qu’au sein de deux groupes.
Nous allons ensuite tester l’indépendance de nos groupes.Pour se faire, nous allons effectuer un test du Khi-2 de Pearson
- l’indépendance des groupes :
test du Khi-deux de l’indépendance des groupes:
\(H_0\) : indépendance des groupes
\(H_1\) : dépendance des groupes
data.frame(
filter(Data_outlier, Codes == 48)$nbre_enregist,
filter(Data_outlier, Codes == 58)$nbre_enregist,
filter(Data_outlier, Codes == 60)$nbre_enregist,
filter(Data_outlier, Codes == 62)$nbre_enregist) %>% chisq.test()##
## Pearson's Chi-squared test
##
## data: .
## X-squared = 13.436, df = 15, p-value = 0.5687Le résultat du test du Chi-2 de Pearson nous donne une \(P-value\) supérieure à \(0.05\), on accepte donc l’hypothèse \(H_0\) selon laquelle les groupes sont indépendants.
Test de comparaison des variances
Nous avons des groupes indépendants qui ne suivent ou n’approchent pas toutes des lois normales. Cela nous amènera à effectuer un test de comparaison des variances. Le test adéquat est celui de Levene qui compare les variances de deux ou plus de deux groupes indépendants n’étant pas normalement distribués.
Les hypothèses de test sont les suivantes:
\(H_0\) : les variances sont égales
\(H_1\) : il y’a au moins une variance qui diffère des autres
test de Levene:
library(car)
leveneTest(nbre_enregist ~ Codes, data = Data_outlier)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.9589 0.4576
## 8La \(P-value\) égale à \(0.4576\) étant supérieure à \(0.05\), on accepte l’hypothèse \(H_0\) d’homogenéité des variances.
Tout ceci nous conduit à pouvoir éffectuer un test de comparaison des moyennes.
Test de comparaison des moyennes / médianes
- Les groupes ne suivent pas toutes des lois normales,
- Les groupes sont indépendants
- Les variances des groupes similaires
- Les effectifs sont trop faibles (n < 6)
- 4 groupes
Le test de Kruskal-Wallis est approprié lorsque les données ne suivent pas de distribution normale ou si les effectifs sont trop faibles (\(n < 6\)) et les variances similaires. Il est basé sur le rang des données et non sur leurs valeurs.
\(H_0\) : les médianes des groupes sont égales.
\(H_1\) : Les médianes des groupes ne sont pas statistiquement égales, il ya une différence entre les groupes.
test de kruskal-Wallis:
kruskal.test(nbre_enregist ~ Codes,data = Data_outlier)##
## Kruskal-Wallis rank sum test
##
## data: nbre_enregist by Codes
## Kruskal-Wallis chi-squared = 0.53846, df = 3, p-value = 0.9104Ce test nous permet de confirmer l’égalité des médianes.
Intéressons nous maintenant aux minimum du nombre de dossiers enregistrés par les patients diabétiques suisses.
Nombre jours sans enregistrements (min)
psych::headTail(filter(Diabet_Swiss_bis, nbre_enregist == min(nbre_enregist)))## dates horaires Codes nbre_enregist
## 1 1991-04-29 12:00 65 0
## 2 1991-04-30 22:00 65 0
## 3 1991-05-01 12:00 65 0
## 4 1991-05-02 13:30 65 0
## ... <NA> <NA> <NA> ...
## 48 1991-08-15 1:30 65 0
## 49 1991-08-16 11:30 65 0
## 50 1991-08-22 20:30 65 0
## 51 1991-09-01 15:15 65 0Nous avons un aperçu du tableau du nombre minimal d’enrégistrements de dossiers, \(0\) . il compte \(51\) lignes, la variable Codes à une seule modalité, la 65 représentant les symptômes d’hypoglycémie.
Le tableau nous dit que sur les \(943\) horaires d’enregistrements de dossiers observés de patients diabetiques suisses,il n’y’a eu que \(51\) sans présence de symptômes d’hypoglycémie. ces enregistrements se sont faites en ligne durant la période du \(29\) mars \(1991\) au \(1^{er}\) septembre de la même année. L’absence d’enregistrement de syptomes d’hypoglycémie pourrait signifier que ces patients n’en ont pas présentés. L’absence de symptômes d’hypoglycémie signifie qu’une personne ne présente pas les symptômes habituels associés à une baisse anormale du taux de glucose sanguin. L’hypoglycémie se caractérise par une diminution du taux de glucose dans le sang et peut être causée par divers facteurs tels que la prise de médicaments pour contrôler le diabète. On peut supposer que les patients suisses diabétiques concernés suivent plus ou moins bien leurs traitements, mais \(51\) sur \(943\) reste relativement faible.
Effectuons un test de comparaison de proportions confirmer ou infirmer cela. Un test binomial exact pour s’avoir si cette proportion est équilibrée ou pas.
test binomial exact:
binom.test( x=c(51), n=(943), p=0.5,alternative="greater")##
## Exact binomial test
##
## data: c(51) and (943)
## number of successes = 51, number of trials = 943, p-value = 1
## alternative hypothesis: true probability of success is greater than 0.5
## 95 percent confidence interval:
## 0.042489 1.000000
## sample estimates:
## probability of success
## 0.05408271La \(P-value\) étant supérieur à \(0.05\), le test nous permet donc de confirmer que cette proportion est toute sauf équilibrée. elle est meme significativement faible. Intéressons nous au maximum
Maximum d’enregistrements de dossiers (max)
# jour particulier du nombre d'enregistrements record
filter(Diabet_Swiss_bis, nbre_enregist == 343)## dates horaires Codes nbre_enregist
## 1 1991-05-10 17:30 62 343Nous n’obtenons qu’une seule ligne sur le maximun. ce qui nous permet aisément de dire que le \(10\) mai \(1991\) à \(17h30\), un nombre d’enregistrements record de \(343\) de la mesure du taux de glycémie avant le souper, a été observé. Nous pourrions à cet effet nous demander si ce jour du \(10\) mai \(1991\) représenterait peut-etre quelque chose de particulier pour ces patients diabetiques suisses.
Nous effectue ensuite le meme test précédant.
test binomial exact:
binom.test(x=c(343),n=943,p=0.5,alternative = "greater")##
## Exact binomial test
##
## data: c(343) and 943
## number of successes = 343, number of trials = 943, p-value = 1
## alternative hypothesis: true probability of success is greater than 0.5
## 95 percent confidence interval:
## 0.3377762 1.0000000
## sample estimates:
## probability of success
## 0.3637328Mêmes résultats concernant cette proportion, qui reste toujours significativement faible comme le confirme de test.
Une fois notre analyse descriptive effectuée, une analyse prédictive sera faite.
Analyse prédictive
Le but de cette analyse, sera de prédire les nombres d’enregistrements de dossiers de patients diabetiques suisses les jours à venir. Pour ce faire, nous allons dans un premier temps regrouper les données et n’avoir qu’une seule ligne par jour. Nous obtiendrons alors un nouveau data.frame constitué de \(2\) variables: dates_daily et nbre_enregist_daily, ce qui nous permettra de d’analyser une série temporelle. Voici ce que ce notre regroupement nous donne comme data.frame appelé Diabet_Suiss_daily
Diabet_Swiss_daily<- Diabet_Swiss_bis %>%
dplyr::select(dates,nbre_enregist) %>%
group_by(dates) %>%
summarize(nbre_enregist_daily = sum(nbre_enregist)) %>%
rename("dates_daily" = "dates")
Diabet_Swiss_daily## # A tibble: 136 × 2
## dates_daily nbre_enregist_daily
## <date> <int>
## 1 1991-04-21 371
## 2 1991-04-22 459
## 3 1991-04-23 417
## 4 1991-04-24 741
## 5 1991-04-25 597
## 6 1991-04-26 340
## 7 1991-04-27 547
## 8 1991-04-28 436
## 9 1991-04-29 700
## 10 1991-04-30 591
## # ℹ 126 more rows# Transformer la variable "dates" en index
rownames(Diabet_Swiss_daily) <- Diabet_Swiss_daily$dates_dailyNotre tableau comprend \(136\) dates.
Ce qui nous donne le nuage de points fonction du temps suivant:
ggplot(Diabet_Swiss_daily, aes(x =dates_daily, y = nbre_enregist_daily)) +
geom_jitter(col="red") +
#geom_freqpoly(binwidth=Diabet_Swiss_daily$dates_daily)+
labs(x = "Dates_daily", y = "Enregistremens_daily",title = "Nuage de points du nombre \n d'enregistrements au quotidien",
subtitle = "du 21 mars au 03 septembre 1991", caption = "Source: uci.edu's data")+
theme_hc(bgcolor ="darkunica")+
theme(plot.title = element_text(size = 20, color = "red"))
On constate que la plupart des points de notre nuage, ont des valeurs comprises environ entre \(250\) et \(750\).
A première vue, nous pourrions imaginer qu’en transformant ces points en une serie temporelle \(X_t\), cette dernière serait stationnaire, c’est à dire, qu’aurait des propirétés de stationnaritté suivantes: Une moyenne constante au cours du temps, une variace constante et une autocovariance ne dépendant pas du temps.
\(E(X_t) = \mu\)
\(Var(X_t)= \sigma^2\)
\(Cov(X_t,X_{t+h}) = \gamma(h))\)
Le but étant de dégager des composantes de notre série temporelle \(X_t\) et de faire des prédictions.
Avec \(X_t = \phi_1(X_{t-1}) + \phi_2(X_{t-2}) +.......+\phi_p(X_{t-p})\)
les \(\phi_i\) des coéfficients , \(i € [ 1;p]\)
Visualisons cette transformation de nos données en série temporelle.
# Transformer la variable "dates" en index
#rownames(df) <- df$dates
Diabet_Swiss_daily_ts <- ts(Diabet_Swiss_daily$nbre_enregist_daily,start = c(annee, mois, jour),
frequency = 142)
autoplot(Diabet_Swiss_daily_ts,col="red") +
labs(title="Serie temporelle du nombre \nd'entregistrements au quotidien ",
subtitle = "du 21 mars au 03 septembre 1991",
caption = "Source: uci.edu's data")+
theme(plot.title = element_text(size = 20, color = "red"))
A première vue, nous pouvons voir que notre serie temporelle est la la limite de la stationnarité . En effet nous avons une moyenne, variance presque constante avec presqu’une absence de la tendance, mais pour nous assurer de la stationnarité de notre série temporelle \(X_t\), nous allons effectuer le test de Dickey-Fuller de la non-stationnalité.
Les hypothèses de test sont les suivantes:
\(H_0\) :serie non-stationnaire.
\(H_1\): serie stationnaire.
test de Dickey-Fuller:
##
## Augmented Dickey-Fuller Test
##
## data: Diabet_Swiss_daily_ts
## Dickey-Fuller = -3.3592, Lag order = 5, p-value = 0.06435
## alternative hypothesis: stationaryLe résultat nous donne une \(P-value\) très limite de \(0.064\), à peine supérieure au seuil de significativité de \(5\)% , Néanmoins l’hypothèse \(H_0\) de non-stationnarité est acceptée. Notre série est donc non-stationnaire.
Nous allons essayer de la rendre stationnaire. La méthode la plus courante pour rendre une série stationnaire est la différenciation. Nous allons lui appliquer une différenciation d’ordre 1 en utilisant la fonction diff() sur la série. La nouvelle série temporelle sera nommée Diabet_Suiss_daily_ts_diff.
Diabet_Swiss_daily_ts_diff <- diff(Diabet_Swiss_daily_ts,differences = 1)Effectuons le même test précédant sur notre nouvelle série temporelle
adf.test(Diabet_Swiss_daily_ts_diff)##
## Augmented Dickey-Fuller Test
##
## data: Diabet_Swiss_daily_ts_diff
## Dickey-Fuller = -7.3272, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationaryLa \(p-value\) étant inférieure à \(0.05\), nous rejettons l’hypothèse \(H_0\) de non-stationnalité. Le test nous dit donc que notre nouvelle série est stationnaire, meilleure que la première. Visualisons-la
autoplot(Diabet_Swiss_daily_ts_diff,col="red") +
labs(title="Serie temporelle différenciée \n du nombre d'enregistrements au quotidien ",
subtitle = "du 21 mars au 03 septembre 1991",caption = "Source: uci.edu's data")+
#theme(plot.background = element_rect(fill = "darkblue"))
#theme(plot.background = element_rect(fill = "skyblue"))+
theme_economist()+
theme(plot.title = element_text(size = 18, color = "black"))
En effet nous voyons que nous avons maintenant une moyenne constante (environ \(0\)), variance constante (série oscillant environ entre \(400\) et \(-400\)), une absence de tendance. Ce qui est en accord avec notre test confirmant la stationarité de notre nouvelle série différenciée.
Distribution de nos deux séries
autoplot(Diabet_Swiss_daily_ts_diff,col="blue") +
labs(title="Serie temporelle différenciée \n du nombre d'enregistrements au quotidien ",
subtitle = "du 21 mars au 03 septembre 1991",caption = "Source: uci.edu's data")+
#theme(plot.background = element_rect(fill = "darkblue"))
#theme(plot.background = element_rect(fill = "skyblue"))+
theme_economist()+
theme(plot.title = element_text(size = 18, color = "black")) 
autoplot(Diabet_Swiss_daily_ts,col="red") +
labs(title="Serie temporelle du nombre d'entregistrements \n au quotidien ",subtitle = "du 21 mars au 03 septembre 1991",
caption = "Source: uci.edu's data")+
theme(plot.title = element_text(size = 20, color = "red")) 
Estimation du modèle
Précédement, pour rendre stationnaire notre série temporelle, nous l’avons différencié, cette différenciation d’ordre 1, va nous permettre de rechercher un modèle ARIMA avec \(d=1\) dans \(ARIMA(p,d,q)\) (\(p\) étant égale aux valeurs des traits noirs sortant de lintervalle délimitée par les traits bleu du \(PCAF\) (l’autocorrélation partielle) et \(q\) celles du \(ACF\) (l’autocorrélation)
acf(Diabet_Swiss_daily_ts_diff)
pacf(Diabet_Swiss_daily_ts_diff)
Nos obtenons comme valeurs possibles de p: \(0,1\) , et q: \(0,1\)
Nous allons donc regarder les modèles . \(ARIMA (1,0,1)\) . \(ARIMA (1,1,1)\) . \(ARIMA (1,0,0)\), \(ARIMA (0,0,1)\) . \(ARIMA (1,1,0)\)
Le modèle \(ARIMA(1,0,1)\) semble etre un bon modèle comme le montre la grpahe ci-dessous:
#model_serie <- tseries::arma(Diabet_Swiss_daily_ts_diff,order = c(1,0,1)) # c(p,d,q)
model_serie<- arma(x = Diabet_Swiss_daily_ts_diff, order = c(1,0,1))
#plot(Diabet_Swiss_daily_ts_diff,ylab="Diabet_Swiss_daily_ts_diff",
# main="ARIMA (1,0,1) VS Série différenciée")
#lines(model_serie$fitted.values,col="red")
autoplot(cbind(Diabet_Swiss_daily_ts_diff,model_serie$fitted.values),ylab="Diabet_Swiss_daily_ts_diff",
main="ARIMA (1,0,1) VS Série différenciée") +
theme(legend.position = "bottom") 
Nous voyons que notre modèle n’est pas si mal que cela. En effet la courbe en verte de notre modèle \(ARIMA(1,0,1)\) a plutot tendance à bien suivre celle de notre série stationaire différenciée, nous avons quelques fluctuations de part et d’autre qui dépaasent un petit peu, mais notre modèle pourrait etre accepter dans un premier temps (nous nous intéresserons après aux résidus) d’autant plus que le summary de notre modèle ci-dessous, nous donne des coéfficiants tous significatifs.
summary(model_serie)##
## Call:
## arma(x = Diabet_Swiss_daily_ts_diff, order = c(1, 0, 1))
##
## Model:
## ARMA(1,0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -507.89 -136.26 16.82 139.50 616.05
##
## Coefficient(s):
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.52668 0.07364 -7.152 8.55e-13 ***
## intercept -2.24250 17.69109 -0.127 0.899
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Fit:
## sigma^2 estimated as 42569, Conditional Sum-of-Squares = 5661690, AIC = 1826.06Avec le Summary nous constatons qu’il s’agit plutot du modèle \(ARMA (1,0)\).
choix automatique du meilleur modèle
Ici il s’agira de choisir un modèle qui minimisera le mieux le critère de sélection AIC , BIC ou HQ
library(simts)
# Appliquer la fonction select_arima
result <- select_arima(Diabet_Swiss_daily_ts)
# Afficher les résultats de sélection du modèle
headTail(result,6)## p d q sigma logLik
## 1 0 0 0 180.29 -899.44
## 2 1 0 0 179.32 -898.7
## 3 2 0 0 177.48 -897.32
## 4 3 0 0 176.14 -896.32
## 5 0 0 1 179.53 -898.86
## 6 1 0 1 175.43 -895.79
## ... ... ... ... ... ...
## 61 0 0 3 177 -896.96
## 62 1 0 3 175.13 -895.56
## 63 2 0 3 174.72 -895.25
## 64 3 0 3 169.62 -892.97
## models
## 1 467.4118, 32505.46, 239.0109, TRUE, -899.4388, 1802.878, 0, 0, 0, 0, 142, 0, 0, -96.41176, -8.411765, -50.41176, 273.5882, 129.5882, -127.4118, 79.58824, -31.41176, 232.5882, 123.5882, -123.4118, -78.41176, 190.5882, -235.4118, 35.58824, -111.4118, -310.4118, 83.58824, -346.4118, 263.5882, -104.4118, -21.41176, 53.58824, -194.4118, -65.41176, -175.4118, -155.4118, -83.41176, -215.4118, -254.4118, -151.4118, -85.41176, 178.5882, 231.5882, 54.58824, 173.5882, 204.5882, -199.4118, 206.5882, 16.58824, -142.4118, 332.5882, -53.41176, 91.58824, -208.4118, 149.5882, 108.5882, 58.58824, 272.5882, 115.5882, 16.58824, 48.58824, 111.5882, 261.5882, 307.5882, -69.41176, -180.4118, 89.58824, -191.4118, 84.58824, -101.4118, 171.5882, -56.41176, 186.5882, -102.4118, 125.5882, -92.41176, 38.58824, 104.5882, -39.41176, 524.5882, 228.5882, 2.588235, 134.5882, 292.5882, 238.5882, -44.41176, 133.5882, 38.58824, 106.5882, 244.5882, -182.4118, 367.5882, 42.58824, -18.41176, 57.58824, 198.5882, 29.58824, -30.41176, -108.4118, 232.5882, -197.4118, -292.4118, 7.588235, -3.411765, -76.41176, -363.4118, -52.41176, -66.41176, 6.588235, 161.5882, 98.58824, -108.4118, 93.58824, 55.58824, -187.4118, 78.58824, -177.4118, -22.41176, -189.4118, 101.5882, -231.4118, -234.4118, -304.4118, -93.41176, -184.4118, -15.41176, 127.5882, -279.4118, -199.4118, 56.58824, -256.4118, -270.4118, -375.4118, -296.4118, 202.5882, 553.5882, -141.4118, 4.588235, 37.58824, -117.4118, 9.588235, 115.5882, 170.5882, -47.41176, -332.4118, arima(x = xt, order = ., include.mean = include.mean), xt, 0, 0, 136, 1, -332.4118, 0, 0, 1, 0, 1
## 2 0.1046515, 467.0438, 32153.98, 0.007419486, -0.02908593, -0.02908593, 294.5332, TRUE, TRUE, -898.705, 1803.41, 1, 0, 0, 0, 142, 0, 0, -95.51647, 2.007283, -49.20205, 279.1933, 101.2862, -140.644, 93.25148, -39.41138, 236.2049, 99.57695, -136.016, -65.16713, 199.1236, -255.0277, 60.55384, -114.8067, -298.4229, 116.4027, -354.83, 300.1702, -131.6672, -10.15551, 56.15842, -199.6904, -44.73687, -168.2369, -136.7253, -66.81828, -206.3532, -231.5392, -124.4578, -69.23689, 187.8561, 213.2281, 30.6816, 168.2049, 186.7514, -220.4928, 227.7864, -4.702114, -143.8183, 347.8213, -87.8882, 97.50727, -217.6672, 171.7282, 93.26302, 47.55373, 266.7863, 87.39089, 4.82117, 47.18167, 106.8328, 250.2398, 280.5421, -101.2719, -172.8183, 108.798, -200.4579, 104.9492, -109.9346, 182.5305, -74.03931, 192.8212, -121.6091, 136.6352, -105.2253, 48.58868, 100.8793, -50.02766, 529.0421, 174.0187, -21.00445, 134.6468, 278.8328, 208.2979, -69.05096, 138.5654, 24.93744, 102.8793, 233.763, -207.6789, 387.0073, 4.448998, -22.53927, 59.84447, 192.891, 9.135097, -33.1788, -104.8997, 244.2631, -221.4231, -271.4229, 38.51897, -3.876471, -75.7253, -355.0857, -14.05077, -60.59738, 13.86774, 161.2282, 82.0072, -118.3998, 105.2631, 46.1235, -192.8997, 98.53057, -185.3067, -3.515947, -186.7369, 121.7399, -241.7137, -209.8648, -279.5508, -61.22521, -174.3067, 4.216613, 129.5305, -292.4346, -169.8415, 77.78638, -262.0044, -243.2485, -346.7834, -256.795, 233.9376, 532.7165, -199.0162, 19.7166, 37.43748, -121.016, 22.20496, 114.9142, 158.8212, -64.93466, -327.1206, arima(x = xt, order = ., include.mean = include.mean), xt, 0, 0, 136, 0.1046515, 1, -332.0438, 0, 0.1046515, 1, 0, 1
## 3 0.09007864, 0.1427166, 466.3854, 31497.8, 0.007325823, -0.0007459563, -0.03508675, -0.0007459563, 0.007289354, -0.03900352, -0.03508675, -0.03900352, 391.0841, TRUE, TRUE, TRUE, -897.3235, 1802.647, 2, 0, 0, 0, 142, 0, 0, -93.88643, 2.610137, -35.10709, 280.1172, 112.9258, -177.3431, 73.35831, -19.60979, 224.8466, 107.9074, -166.9512, -84.1457, 216.0518, -240.6016, 30.38109, -80.23293, -304.6676, 128.2374, -308.853, 283.6505, -77.92932, -48.83753, 71.20567, -195.3957, -54.75996, -140.9864, -129.4882, -43.59091, -184.9309, -222.3161, -96.96446, -34.67661, 208.6784, 228.4783, 9.026962, 136.4069, 181.9484, -241.8273, 196.1402, 27.22582, -172.6022, 343.8365, -62.25894, 49.72101, -208.2518, 156.0779, 125.6448, 28.24542, 252.6007, 83.45974, -31.93928, 31.38502, 105.6315, 245.3896, 268.8866, -133.6645, -217.2698, 116.5331, -172.9466, 89.83201, -80.92632, 169.4385, -56.60766, 167.9686, -110.3811, 108.9715, -88.32132, 29.77642, 115.0884, -53.55272, 513.9993, 187.7461, -92.08276, 102.5191, 280.8827, 193.8117, -106.8733, 104.3257, 33.68048, 84.83439, 230.2671, -218.8685, 349.9002, 36.29698, -73.92162, 53.95609, 196.8158, 4.268271, -60.63148, -109.1077, 247.4815, -202.1034, -307.036, 62.88963, 38.42412, -76.40001, -355.2544, -7.983502, -9.038282, 20.83794, 171.2602, 83.87973, -139.5664, 90.07103, 63.41749, -204.9883, 88.32406, -156.9567, -16.85921, -161.2859, 122.6361, -212.7431, -227.2774, -249.4826, -31.74891, -131.7653, 15.3186, 156.0825, -287.9178, -191.6643, 115.2151, -232.2624, -254.6032, -313.6718, -223.2155, 283.6535, 578.4296, -219.4036, -60.89245, 58.14414, -120.6651, 15.58746, 132.2685, 159.5952, -78.48709, -351.6994, arima(x = xt, order = ., include.mean = include.mean), xt, 0, 0, 136, 0.09007864, 0.1427166, 0, 1, 0, -331.3854, -6.619976, 0, 0, 0, 0, 0.09007864, 0.1427166, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0
## 4 0.07063202, 0.131957, 0.1221506, 465.9212, 31025.8, 0.007400836, -0.0006537756, -0.001172072, -0.03811937, -0.0006537756, 0.007235869, -0.0006511611, -0.04460104, -0.001172072, -0.0006511611, 0.007363384, -0.03144884, -0.03811937, -0.04460104, -0.03144884, 493.5615, TRUE, TRUE, TRUE, TRUE, -896.3191, 1802.638, 3, 0, 0, 0, 142, 0, 0, -92.75279, 2.779607, -34.50478, 291.0422, 118.9503, -165.5024, 39.07504, -35.04315, 240.8746, 102.5898, -157.9892, -113.4076, 198.3218, -222.4451, 37.65093, -105.1352, -277.4765, 116.8742, -296.7393, 315.9494, -86.52193, -5.498376, 37.6875, -181.6109, -55.12946, -150.6769, -109.6365, -40.29133, -166.5794, -208.1999, -93.82177, -13.82666, 236.6839, 249.7464, 26.10431, 118.3647, 157.8419, -242.4299, 173.4788, 4.326195, -145.4794, 316.2298, -59.13068, 69.87564, -247.4521, 159.7538, 115.3428, 57.64332, 236.8553, 76.34603, -33.69602, -0.1264178, 92.85474, 246.2752, 269.4583, -138.2798, -247.0442, 74.92487, -164.4478, 109.3301, -92.06504, 191.9767, -64.4754, 181.3244, -128.1, 116.0974, -109.5537, 42.05937, 99.72287, -39.59639, 509.8637, 184.9671, -76.95985, 41.16934, 255.8248, 200.8526, -115.3063, 70.50841, 6.875877, 92.66618, 216.6563, -217.4597, 336.1839, 11.82515, -46.63751, 9.37426, 192.7545, 11.21788, -64.73469, -133.4193, 241.6509, -194.8129, -294.9108, 26.88744, 59.7585, -40.44735, -357.4849, -15.23695, -4.41485, 63.59258, 177.295, 95.4243, -136.4962, 69.50457, 52.24749, -189.4386, 74.06489, -164.016, 3.647867, -173.0112, 140.6016, -209.8488, -207.3287, -268.7209, -11.7047, -108.0046, 48.13051, 165.428, -262.8574, -193.6235, 92.96495, -198.9581, -234.4033, -328.3825, -201.8856, 307.1001, 625.2559, -170.0324, -82.21319, -10.69019, -102.3921, 13.36727, 126.8194, 176.5072, -74.87815, -364.686, arima(x = xt, order = ., include.mean = include.mean), xt, 0, 0, 136, 0.07063202, 0.131957, 0.1221506, 0, 0, 1, 0, 0, -330.9212, 14.95989, -5.609305, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.07063202, 0.131957, 0.1221506, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0
## 5 0.08279662, 467.1511, 32230.01, 0.00600827, -0.01845828, -0.01845828, 277.5979, TRUE, TRUE, -898.8635, 1803.727, 0, 1, 0, 0, 142, 0, 0, -95.82321, -0.2443102, -50.13086, 277.9996, 106.8315, -135.9964, 91.10894, -38.69461, 236.0527, 104.3045, -131.7872, -67.23957, 196.4161, -251.4137, 56.6651, -115.8428, -300.5597, 108.7342, -355.1539, 293.2544, -128.4316, -10.5174, 54.71971, -198.6817, -48.70093, -171.1188, -140.983, -71.47818, -209.2329, -236.8273, -131.5426, -74.25982, 184.9974, 216.5317, 36.92081, 170.792, 190.7079, -214.9411, 224.6453, -1.750969, -142.0061, 344.6065, -81.68335, 98.61201, -216.3158, 167.7591, 94.95901, 50.98662, 268.6274, 93.60746, 9.09852, 48.09557, 107.8668, 252.9179, 286.9082, -92.90612, -172.4588, 104.1279, -199.7725, 101.3894, -109.5458, 180.9189, -71.13057, 192.7383, -118.1092, 135.6279, -103.3806, 47.40847, 100.9236, -47.50723, 528.7823, 185.0675, -12.47406, 135.8817, 281.5984, 215.5335, -61.99654, 138.982, 27.34166, 104.5851, 236.1896, -201.7068, 384.5495, 11.0095, -19.06265, 59.42722, 193.9285, 13.79228, -31.29305, -105.5601, 241.5889, -217.1538, -274.1715, 30.54937, -5.680483, -75.68077, -356.885, -22.60223, -64.27971, 12.17104, 160.8412, 85.5318, -115.2328, 103.3898, 47.28858, -191.0664, 94.66856, -184.9893, -6.834607, -188.5852, 117.4631, -240.8766, -214.2073, -286.4155, -69.43687, -178.402, -0.3800194, 127.8804, -289.7392, -175.1617, 71.3517, -262.0588, -248.4535, -354.58, -266.7931, 224.9385, 535.2248, -185.4659, 20.20485, 36.17601, -120.1463, 19.79661, 114.2098, 161.3927, -60.51387, -327.1408, arima(x = xt, order = ., include.mean = include.mean), xt, 0, 0, 136, 0.08279662, 1, 0, -332.1511, -27.08615, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0.08279662, 0.08279662, 0.00685528, 0, 1, 0.08279662, 0.08279662, 0.00685528
## 6 0.9079088, -0.811684, 466.0715, 30775.31, 0.006681127, -0.008362333, -0.07503975, -0.008362333, 0.01275689, 0.06629564, -0.07503975, 0.06629564, 876.0316, TRUE, TRUE, TRUE, -895.7913, 1799.583, 1, 1, 0, 0, 142, 0, 0, -92.66133, 5.844167, -37.59064, 287.3417, 112.4871, -153.5969, 70.94161, -46.01643, 223.744, 93.95422, -159.2188, -95.40938, 184.4499, -258.6069, 39.55549, -111.4899, -299.6186, 122.3463, -322.8702, 316.1541, -86.98687, 2.902536, 75.50734, -181.6533, -36.22503, -145.3036, -113.9709, -34.69695, -167.721, -194.8505, -78.46267, -11.50715, 246.9176, 269.9889, 63.59629, 175.7706, 189.7796, -230.9948, 200.2646, -8.300051, -164.086, 328.8226, -88.34808, 68.49394, -235.8467, 147.4976, 92.6206, 35.3021, 248.173, 69.66448, -31.68637, 7.931732, 74.03603, 220.4936, 249.1845, -146.2913, -236.0111, 61.94264, -222.3484, 78.01941, -114.7596, 170.6357, -73.57252, 178.2108, -127.0426, 115.5737, -112.5017, 31.29732, 95.08059, -57.06952, 514.1715, 169.7782, -67.01937, 77.96322, 233.7993, 162.8394, -128.7306, 69.54494, -26.12576, 50.4712, 188.9059, -251.0202, 329.5758, -23.51353, -76.04013, 12.7073, 156.7411, -23.3641, -76.11603, -142.4594, 215.5076, -233.5331, -302.6116, 27.56989, 12.20027, -63.288, -345.2833, -2.604525, -20.81729, 50.11042, 196.404, 111.4222, -107.3578, 104.999, 55.96811, -192.329, 92.7541, -173.3523, -1.921909, -170.5005, 135.2877, -213.7103, -197.6528, -251.8954, -21.36973, -116.8244, 57.31619, 188.2268, -242.3462, -142.3165, 122.2434, -208.4423, -206.6791, -297.5373, -196.9549, 311.9613, 622.994, -138.2217, 20.90832, 50.51691, -110.4113, 26.69177, 128.6717, 170.2088, -64.01111, -341.1996, arima(x = xt, order = ., include.mean = include.mean), xt, 0, 0, 136, 0.9079088, -0.811684, 1, 0, -331.0715, 276.9462, 0, 0, 0, 3.330669e-16, 0.9079088, 0, 1, 0, 1, -0.811684, -0.811684, 0.6588309, 0, 1, -0.811684, -0.811684, 0.6588309
## ... NA
## 61 0.0763169, 0.1127197, 0.1062847, 466.2381, 31329.94, 0.007672174, 0.0006662939, 0.0006279602, -0.03320899, 0.0006662939, 0.007693495, -0.001081994, -0.03832471, 0.0006279602, -0.001081994, 0.007290043, -0.03083915, -0.03320899, -0.03832471, -0.03083915, 384.2852, TRUE, TRUE, TRUE, TRUE, -896.9641, 1803.928, 0, 3, 0, 0, 142, 0, 0, -93.84875, 1.705632, -38.04701, 287.0752, 112.7028, -163.212, 50.02816, -27.63723, 247.5784, 103.6655, -155.119, -103.3988, 206.1199, -221.8266, 41.447, -110.3043, -281.9151, 114.3051, -310.4604, 305.5342, -103.7094, -13.76579, 35.02898, -183.337, -52.73171, -153.2711, -117.1111, -50.41925, -180.8991, -221.302, -107.5993, -31.8545, 217.8426, 231.1637, 16.95075, 124.2583, 169.7991, -227.0046, 192.7398, 10.59348, -139.645, 322.7398, -62.25374, 75.97595, -240.3214, 167.1552, 116.0191, 57.60855, 238.5218, 79.73403, -21.33212, 17.05113, 105.3907, 255.0641, 275.6044, -129.2235, -227.5516, 93.40148, -157.9821, 111.4757, -100.865, 184.6852, -69.81136, 182.9925, -126.9636, 123.2444, -105.7817, 47.43709, 100.9664, -40.04761, 512.3955, 184.4405, -63.8146, 65.38218, 276.3621, 218.0835, -97.9821, 88.28425, 20.88993, 106.6303, 225.8863, -212.7166, 348.2008, 17.15743, -35.18802, 22.50497, 200.1872, 16.68747, -55.46857, -126.1627, 247.8691, -195.0382, -290.8839, 26.60125, 49.24975, -51.07861, -366.7186, -22.72816, -16.73822, 50.57786, 163.2044, 83.38459, -137.3737, 78.50066, 57.39323, -184.866, 79.05757, -167.5335, 2.284666, -177.9307, 133.8898, -220.6426, -212.58, -276.3743, -23.73313, -127.68, 27.55565, 143.5735, -278.7308, -196.0785, 88.88483, -210.2947, -242.3679, -341.484, -219.5063, 284.766, 594.0667, -194.344, -76.63563, 3.376706, -87.20161, 25.18148, 124.3106, 168.7047, -75.80174, -357.6818, arima(x = xt, order = ., include.mean = include.mean), xt, 0, 0, 136, 0.0763169, 0.1127197, 0.1062847, 1, 0, 0, 0, -331.2381, -17.91078, -48.37435, -38.01611, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0.0763169, 0.1127197, 0.1062847, 0.0763169, 0.00582427, 0.008602419, 0.008111322, 0.1127197, 0.008602419, 0.01270573, 0.01198038, 0.1062847, 0.008111322, 0.01198038, 0.01129645, 0, 1, 0.0763169, 0.1127197, 0.1062847, 0.0763169, 0.00582427, 0.008602419, 0.008111322, 0.1127197, 0.008602419, 0.01270573, 0.01198038, 0.1062847, 0.008111322, 0.01198038, 0.01129645
## 62 0.8828741, -0.8310331, 0.06502893, -0.002400208, 466.223, 30669.95, 0.01586431, -0.01593909, 0.000122464, -0.006301754, -0.04309387, -0.01593909, 0.02356566, -0.006472375, 0.006786952, -0.004366212, 0.000122464, -0.006472375, 0.01364211, -0.007361336, -0.02676764, -0.006301754, 0.006786952, -0.007361336, 0.01118148, 0.05363356, -0.04309387, -0.004366212, -0.02676764, 0.05363356, 828.2584, TRUE, TRUE, TRUE, TRUE, TRUE, -895.564, 1803.128, 1, 3, 0, 0, 142, 0, 0, -92.65278, 2.158642, -34.79692, 287.2742, 127.3703, -154.4766, 56.4218, -44.33471, 219.5197, 103.7506, -160.5376, -108.9269, 180.1189, -247.1548, 26.20169, -104.4133, -300.9764, 114.5138, -305.584, 307.446, -61.34343, -0.7948306, 76.69785, -177.9413, -46.496, -144.4065, -117.8161, -34.69359, -163.1471, -197.6984, -80.4269, -5.967896, 253.9313, 285.2769, 70.81159, 166.4379, 185.8663, -236.0905, 174.8968, -4.520103, -172.6142, 315.7247, -73.3143, 57.01117, -236.2299, 133.5306, 103.1269, 39.30923, 247.283, 78.25801, -36.27358, -0.5579035, 70.91321, 222.0896, 256.729, -141.7559, -252.9563, 48.62772, -213.8472, 72.23651, -101.8994, 171.3688, -58.55041, 176.4861, -116.1218, 108.0258, -105.4032, 25.41835, 98.89595, -51.33088, 510.4953, 193.396, -71.68884, 61.51555, 230.1505, 167.4992, -130.5373, 54.11678, -25.35055, 47.7593, 192.0915, -241.7454, 315.4988, -3.43519, -79.82412, 8.626902, 160.2362, -13.19195, -77.75745, -144.7993, 213.1332, -216.2684, -311.9164, 21.25267, 27.45419, -52.57574, -341.237, -11.52015, -7.50908, 59.0505, 205.4445, 122.9385, -106.3656, 93.5466, 58.05308, -194.4444, 79.04809, -168.1806, -11.01106, -167.51, 130.0602, -202.0112, -206.7028, -255.6446, -24.00762, -105.6247, 60.70983, 198.597, -231.0776, -157.3891, 117.4908, -198.9139, -217.2153, -303.8289, -203.6753, 314.3971, 648.6575, -111.8991, -4.842509, 38.48594, -118.4288, 12.4547, 125.4062, 171.8002, -63.23398, -353.8344, arima(x = xt, order = ., include.mean = include.mean), xt, 0, 0, 136, 0.8828741, -0.8310331, 0.06502893, -0.002400208, 1, 0, 0, 0, -331.223, 289.5237, -22.8577, 0.8492761, 0, 0, 0, 0, 0, 3.330669e-16, -2.081668e-17, 8.673617e-19, 0, -2.081668e-17, 1.734723e-18, -8.131516e-20, 0, 8.673617e-19, -8.131516e-20, 2.541099e-21, 0.8828741, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, -0.8310331, 0.06502893, -0.002400208, -0.8310331, 0.690616, -0.05404119, 0.001994652, 0.06502893, -0.05404119, 0.004228761, -0.0001560829, -0.002400208, 0.001994652, -0.0001560829, 5.760996e-06, 0, 1, -0.8310331, 0.06502893, -0.002400208, -0.8310331, 0.690616, -0.05404119, 0.001994652, 0.06502893, -0.05404119, 0.004228761, -0.0001560829, -0.002400208, 0.001994652, -0.0001560829, 5.760996e-06
## 63 0.05928266, 0.7049654, 0.0009898226, -0.6300158, 0.08853557, 466.1766, 30525.69, 0.06130342, -0.03200443, -0.06044282, 0.01983367, -0.0002904664, 0.03026737, -0.03200443, 0.05399963, 0.03119984, -0.04999629, -0.007688039, -0.004747444, -0.06044282, 0.03119984, 0.06685424, -0.01876784, -0.004083981, -0.08162699, 0.01983367, -0.04999629, -0.01876784, 0.05288049, 0.008659458, -0.04844321, -0.0002904664, -0.007688039, -0.004083981, 0.008659458, 0.008421371, -0.003528499, 0.03026737, -0.004747444, -0.08162699, -0.04844321, -0.003528499, 805.5597, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, -895.254, 1804.508, 2, 3, 0, 0, 142, 0, 0, -92.42032, 2.746601, -37.1557, 291.1001, 122.5484, -142.1444, 47.65696, -46.27632, 221.0989, 98.66793, -151.1518, -115.208, 178.673, -250.4817, 38.45616, -120.92, -282.0814, 101.5075, -299.3455, 314.7007, -73.42149, 24.08717, 54.61331, -160.5807, -58.94025, -140.1341, -121.3863, -33.19641, -164.6512, -192.5523, -84.78377, -3.442574, 254.3189, 286.5904, 75.4978, 165.347, 178.1335, -236.3121, 172.2952, -19.61132, -159.2514, 302.1758, -71.33561, 65.12814, -247.6571, 145.2607, 84.99673, 60.34637, 233.484, 88.67981, -40.47019, 1.648424, 63.95505, 225.5694, 253.6294, -135.5669, -252.8909, 41.8934, -216.6118, 82.06797, -111.4557, 189.2519, -72.47308, 198.4311, -136.0248, 131.9782, -130.7656, 51.14331, 73.6191, -28.79858, 495.3665, 200.4126, -66.04843, 56.0501, 223.6654, 167.5927, -128.7449, 54.22694, -33.73418, 56.01245, 185.2477, -233.6694, 318.2482, -14.24963, -58.57876, -8.148101, 172.8092, -22.60931, -62.25617, -156.6587, 223.68, -227.8895, -289.3668, 1.291659, 40.43918, -54.87477, -330.7682, -14.53378, -10.33719, 67.90344, 203.0138, 128.1501, -106.1161, 93.67273, 48.46445, -188.0301, 73.22803, -172.4859, -4.052228, -177.8709, 141.802, -215.4573, -186.7194, -275.198, -8.110181, -120.8228, 81.03877, 183.3145, -214.2482, -163.9733, 114.6294, -203.3475, -207.8755, -316.3836, -195.8824, 304.3754, 655.1304, -108.301, 8.903786, 11.0553, -107.3964, -3.375481, 129.4454, 164.5213, -57.03057, -357.3217, arima(x = xt, order = ., include.mean = include.mean), xt, 0, 0, 136, 0.05928266, 0.7049654, 0.0009898226, -0.6300158, 0.08853557, 1, 0, 0, 0, -331.1766, 17.58955, 220.0691, -31.63568, 0, 0, 0, 0, 0, 1.110977e-16, -3.783866e-17, 3.374579e-18, 0, -3.783866e-17, 1.110223e-16, -1.387779e-17, 0, 3.374579e-18, -1.387779e-17, 1.734723e-18, 0.05928266, 0.7049654, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0.0009898226, -0.6300158, 0.08853557, 0.0009898226, 9.797489e-07, -0.0006236039, 8.763451e-05, -0.6300158, -0.0006236039, 0.3969199, -0.05577881, 0.08853557, 8.763451e-05, -0.05577881, 0.007838547, 0, 1, 0.0009898226, -0.6300158, 0.08853557, 0.0009898226, 9.797489e-07, -0.0006236039, 8.763451e-05, -0.6300158, -0.0006236039, 0.3969199, -0.05577881, 0.08853557, 8.763451e-05, -0.05577881, 0.007838547
## 64 -0.4436003, 0.2157841, 0.8792841, 0.533486, -0.02919572, -0.7822268, 466.1439, 28769.66, 0.01020479, 0.01134744, 0.007544682, -0.01274655, -0.01394404, -0.01134984, -0.02508999, 0.01134744, 0.01541189, 0.01091507, -0.01475812, -0.01917204, -0.01457781, -0.05125287, 0.007544682, 0.01091507, 0.008767541, -0.009733372, -0.01457419, -0.01053039, -0.05046453, -0.01274655, -0.01475812, -0.009733372, 0.01986632, 0.0220688, 0.01604608, -0.008275135, -0.01394404, -0.01917204, -0.01457419, 0.0220688, 0.03030313, 0.02124384, 0.02748648, -0.01134984, -0.01457781, -0.01053039, 0.01604608, 0.02124384, 0.01818349, 0.009969928, -0.02508999, -0.05125287, -0.05046453, -0.008275135, 0.02748648, 0.009969928, 847.5085, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, -892.9701, 1801.94, 3, 3, 0, 0, 142, 0, 0, -89.59615, 1.663622, -31.88218, 278.6398, 119.8131, -151.8623, 47.77381, -23.62039, 207.7896, 89.87418, -139.5384, -126.0547, 195.1016, -233.0694, -6.354074, -70.27731, -295.6221, 81.33125, -249.2155, 266.2464, -72.44111, 39.6978, 14.88144, -132.5044, -43.24261, -178.3516, -57.92161, -66.96828, -166.1274, -153.5157, -121.3168, 20.13408, 259.1073, 231.8279, 99.2737, 147.0999, 170.1858, -196.3492, 142.9507, 19.25297, -156.5561, 273.4074, -23.79564, 25.1018, -248.1043, 194.5385, 46.1526, 50.2708, 266.3515, 26.84492, -5.936834, 1.615977, 47.4184, 254.0008, 224.678, -155.7859, -217.4014, 35.32391, -196.8023, 75.10092, -119.6531, 188.6588, -77.21501, 167.8476, -102.0891, 89.5837, -98.39824, 36.82775, 76.93345, -34.4044, 495.5968, 174.2733, -72.13676, 53.5467, 252.7179, 149.3742, -145.4603, 82.83371, -35.0451, 42.91792, 204.8726, -261.3099, 316.4147, 12.45963, -114.0069, 23.04254, 182.7001, -61.32045, -52.80619, -134.1375, 185.0564, -185.6098, -330.9533, 27.18408, 65.81772, -111.5104, -317.9232, 20.68141, -52.46354, 89.443, 190.827, 88.94277, -74.76388, 72.84912, 61.64926, -173.7647, 51.71428, -135.322, -14.3422, -185.6089, 169.9603, -231.263, -206.3968, -213.2906, -48.10286, -95.57526, 73.59985, 162.3193, -212.7454, -160.422, 119.0476, -174.9026, -246.932, -270.884, -179.1478, 281.5609, 667.961, -159.6162, -28.13762, 111.444, -159.901, 12.33685, 185.836, 100.7605, -41.51353, -319.5147, arima(x = xt, order = ., include.mean = include.mean), xt, 0, 0, 136, -0.4436003, 0.2157841, 0.8792841, 0.533486, -0.02919572, -0.7822268, 1, 0, 0, 0, -331.1439, -104.6272, 2.966469, 248.4109, 0, 0, 0, 0, 0, 0.01225364, -0.001536538, -0.006295875, 0, -0.001536538, 0.007048276, -0.004573189, 0, -0.006295875, -0.004573189, 0.007429626, -0.4436003, 0.2157841, 0.8792841, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0.533486, -0.02919572, -0.7822268, 0.533486, 0.2846073, -0.01557551, -0.4173071, -0.02919572, -0.01557551, 0.0008523901, 0.02283768, -0.7822268, -0.4173071, 0.02283768, 0.6118788, 0, 1.012292, 0.5318958, -0.03547282, -0.7822268, 0.5318958, 0.2917316, -0.02017528, -0.4173071, -0.03547282, -0.02017528, 0.008291318, 0.02283768, -0.7822268, -0.4173071, 0.02283768, 0.6118788
## ic value minval
## 1 AIC 1802.88 FALSE
## 2 AIC 1803.41 FALSE
## 3 AIC 1802.65 FALSE
## 4 AIC 1802.64 FALSE
## 5 AIC 1803.73 FALSE
## 6 AIC 1799.58 TRUE
## ... <NA> ... <NA>
## 61 HQ 1809.85 FALSE
## 62 HQ 1810.23 FALSE
## 63 HQ 1812.79 FALSE
## 64 HQ 1811.41 FALSE# Sélectionner le meilleur modèle
best_model <- best_model(result)
best_model##
## Call:
## arima(x = xt, order = c(1, 0, 1), include.mean = include.mean)
##
## Coefficients:
## ar1 ma1 intercept
## 0.9079 -0.8117 466.0715
## s.e. 0.0817 0.1129 29.5978
##
## sigma^2 estimated as 30775: log likelihood = -895.79, aic = 1799.58Nous voyons que ces 3 critères de sélection des paramères \(p,d,q\) nous permettent d’obtenir le meilleur modèle \(ARIMA(1,0,1)\) avec un \(aic\) minimal de \(1799.58\) .
Validation du modèle
Maintenant que nous avons obtenus notre modèle \(ARIMA(1,0,1)\), nous allons nous intéresser aux résidus pour savoir si nous pourrons valider ou pas notre modèle.
Vérification des conditions sur les résidus
Nos résidus doivent etre un bruit blanc c’est à dire, vérifée les conditions suivantes:
Non-autocorrélation des rösidus
Normalité des résidus
Homocédaticité des résidus
1-Autocorrélation des résidus:
- Best modèle
acf(residuals(best_model),na.action = na.pass)
pacf(best_model$residuals,na.action = na.pass)
2- Normalité des résidus
Quant à elle se vérifie avec un test de Shapiro-wilk
shapiro.test(best_model$residuals)##
## Shapiro-Wilk normality test
##
## data: best_model$residuals
## W = 0.98272, p-value = 0.08275Une \(P-value\) de \(0.08275\) , supérieure au seuil de significativité de \(0.05\) , les résidus sont donc normalement distribués.
3 - Homocédaticité des résidus:
Nous allons vérifier cela visuellement
checkresiduals(best_model$residuals)
##
## Ljung-Box test
##
## data: Residuals
## Q* = 14.01, df = 27, p-value = 0.9812
##
## Model df: 0. Total lags used: 27Nous remarquons bien que nos résidus ne suivent aucune tendance, ce qui nous permet de conclure à une homocédaticité de nos résidus. A noter que le test de Ljung-Box, avec une \(P-value\) de \(0.9812\) largement supérieur au seuil de significativité de \(5\)%, nous a aussi permis de vérifier l’indépendance des résidus.
Suite à la vérification des conditions sur les résidus (bruit blanc), nous pouvons donc valider notre modèle.
Ainsi La formule mathématique de notre modèle \(ARIMA(1,0,1)\) est la suivante :
\(X_t \sim r * X_{t-1} + e_t + a * e_{t-1} + Cste\)
où :
\(X_t\) est la valeur de la série temporelle à l’instant \(t\).
\(r\) est le paramètre autorégressif (autoregressive parameter).
\(X_{t-1}\) est la valeur de la série temporelle à l’instant \(t-1\).
\(e_t\) est le terme d’erreur pure à l’instant \(t\).
\(a\) est le paramètre de moyenne mobile (moving average parameter).
\(e_{t-1}\) est le terme d’erreur pure à l’instant \(t-1\) \(Cste\) une constante
soit
\(X_t \sim 0.9079 * X_{t-1} + e_t -0.8117 * e_{t-1}+ 466.0715\)
Intéressons maintenant aux prédictions.
Prédictions et évaluation
C’était l’un de nos objectifs de départ, nous y sommes.
prédiction avec le modèle arima
#Mydata_predict <- predict(Best_model_serie,n.ahead = 27)
#Mydata_predict$pred
m_arima <- auto.arima(Diabet_Swiss_daily_ts)
summary(m_arima)## Series: Diabet_Swiss_daily_ts
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.9079 -0.8117 466.0715
## s.e. 0.0817 0.1129 29.5978
##
## sigma^2 = 31469: log likelihood = -895.79
## AIC=1799.58 AICc=1799.89 BIC=1811.23
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.2344868 175.4289 141.749 -21.31896 42.26898 NaN -0.04165921prediction_arima <- forecast(m_arima)
autoplot(prediction_arima)
Nous obtenons les prédictions pour les \(2\) années à venir d’un \(ARIMA(1,0,1)\) avec une moyenne de \(466.0715\) et une racine de l’erreur quadratique \(RMSE\) de \(175.4289\).
Les \(5\) premières et les \(4\) dernières lignes ci-dessous de nos prédictions. \(Lo.80\) et \(Lo.95\) les bords du plus petit intervalle de confiance,bleu. \(Hi.80\) et \(Hi.95\) ceux du plus grand, bleu ciel.
Tab_pred_arima <- headTail(data.frame(prediction_arima),5)
Tab_pred_arima## Point.Forecast Lo.80 Hi.80 Lo.95 Hi.95
## 1991.979 442.44 215.09 669.78 94.74 790.13
## 1991.986 444.61 216.22 673 95.32 793.91
## 1991.993 446.59 217.33 675.84 95.97 797.2
## 1992.000 448.38 218.42 678.35 96.68 800.08
## 1992.007 450.01 219.47 680.56 97.42 802.6
## ... ... ... ... ... ...
## 1993.951 466.07 232.82 699.33 109.34 822.81
## 1993.958 466.07 232.82 699.33 109.34 822.81
## 1993.965 466.07 232.82 699.33 109.34 822.81
## 1993.972 466.07 232.82 699.33 109.34 822.81prédictions avec le modèle tbats
m_tbats <- tbats(Diabet_Swiss_daily_ts)
summary(m_tbats)## Length Class Mode
## lambda 1 -none- numeric
## alpha 1 -none- numeric
## beta 0 -none- NULL
## damping.parameter 0 -none- NULL
## gamma.values 0 -none- NULL
## ar.coefficients 0 -none- NULL
## ma.coefficients 0 -none- NULL
## likelihood 1 -none- numeric
## optim.return.code 1 -none- numeric
## variance 1 -none- numeric
## AIC 1 -none- numeric
## parameters 2 -none- list
## seed.states 1 -none- numeric
## fitted.values 136 ts numeric
## errors 136 ts numeric
## x 136 -none- numeric
## seasonal.periods 0 -none- NULL
## y 136 ts numeric
## call 2 -none- call
## series 1 -none- character
## method 1 -none- characterprediction_tbats <- forecast(m_tbats)
autoplot(prediction_tbats)
predict(m_tbats,5)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1991.979 404.557 201.6781 642.5879 112.4601 780.4777
## 1991.986 404.557 200.8232 643.7786 111.3542 782.4028
## 1991.993 404.557 199.9734 644.9644 110.2565 784.3202
## 1992.000 404.557 199.1285 646.1454 109.1671 786.2301
## 1992.007 404.557 198.2884 647.3215 108.0856 788.1326Ici nous avons des points forcast égaux ( \(404.557\) ) inférieurs à ceux préditent par le modèle arima.
Nous allons regarder les erreurs des résidus
- Erreurs du modèle m_tbats
accuracy(m_tbats)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 5.518996 178.0785 142.7586 -19.03762 41.29801 0.722625 -0.01032437- Erreurs du modèle arima
accuracy(m_arima)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.2344868 175.4289 141.749 -21.31896 42.26898 NaN -0.04165921le \(RMSE\) du modèle arima est légèrement meilleure.
regardons cela de près
checkresiduals(m_tbats)
##
## Ljung-Box test
##
## data: Residuals from TBATS
## Q* = 15.707, df = 27, p-value = 0.9583
##
## Model df: 0. Total lags used: 27checkresiduals(m_arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,1) with non-zero mean
## Q* = 14.01, df = 25, p-value = 0.9616
##
## Model df: 2. Total lags used: 27Certes le \(RMSE\) du modèle \(ARIMA\) est meilleur mais lorsque nous regardons la fonction d’autocorrélation et la variance des résidus, c’est plutot le modèle \(tbats\) qui est meilleur, pour la normalité, les deux sont quasi-similaires, nous aurions pu faire un \(test de Kolmogorov\) pour comparer les deux distributions pour savoir si elles sont similaires, mais nous ne le ferons pas, car cela n’apportera pas grand-chose comme information. Nous allons donc choisir de garder les points forcast du modèle tbats Nousn prévoyons donc que pour les deux années à venir, le nombre d’enrégistrements de dossiers de patients diabétiques suisses soit constant au quotidien et soit d’environ \(404\).