Introduction au guide de programmation R : De la théorie à la pratique

Voici un guide pratique de programmation en R, conçu pour accompagner les trois activités du site En avant Math !. Les contenus théoriques et les explications pédagogiques détaillées ont été présentés sur le site. L’objectif de ce guide est de fournir des programmes reproductibles, directement utilisables, afin de permettre aux apprenants d’appliquer concrètement les méthodes d’analyse des séries chronologiques en santé publique.

Les exemples s’appuient principalement sur une série classique de cas de rougeole à Londres, utilisée comme fil conducteur. Les blocs de code peuvent être adaptés à d’autres contextes (autres maladies, autres fréquences temporelles, autres régions).

Pour faire suite à ce qui a été abordé dans les 3 activités, voici ce que vous allez apprendre dans ce guide :

  • Activité 1 : L’exploration. Vous apprendrez à structurer une série temporelle, à visualiser sa dépendance interne et à décomposer son mécanisme afin d’isoler la tendance et la saisonnalité.

  • Activité 2 : La modélisation. Nous passerons à l’utilisation de modèles statistiques (ARMA, GARCH, HMM) pour modéliser la mémoire du passé et la volatilité des épidémies.

  • Activité 3 : La prévision et l’espace. Nous aborderons les modèles ARIMA/SARIMA pour gérer les séries non stationnaires et introduirons la dimension spatiale afin de comprendre comment une maladie se propage entre régions voisines.

Découvrir l’environnement de travail : R et RStudio

Pour réaliser les analyses présentées dans ce guide, nous utilisons le logiciel R, un langage de programmation puissant et gratuit, spécialisé en statistiques et en graphiques. Pour une expérience plus fluide, nous vous recommandons d’utiliser l’interface RStudio, qui permet de visualiser simultanément votre code, vos graphiques et vos résultats.

Comment essayer ces codes ?

• Copier-coller : Tous les blocs de code présentés dans ce document peuvent être copiés et collés directement dans votre script RStudio.

• Expérimentation : Nous vous encourageons vivement à ne pas seulement exécuter le code, mais à modifier les paramètres (comme les ordres \(p, d, q\) des modèles ou la plage de dates) pour observer l’impact direct sur les résultats.

Installation et chargement des outils

Pour analyser des données avec R, nous utilisons des « librairies ». Ce sont des extensions spécialisées qui contiennent des fonctions prêtes à l’emploi et qui étendent les capacités de base du logiciel.

L’utilisation de ces outils se déroule en deux étapes :

  • L’installation : Cette étape télécharge les outils depuis les serveurs officiels sur votre ordinateur. Elle ne doit être effectuée qu’une seule fois.

  • Le chargement : Cette étape active les outils pour votre session de travail actuelle. Elle doit être répétée à chaque ouverture de RStudio. 

# --- INSTALLATION ---
# Cette étape nécessite une connexion internet. 
# Enlevez le symbole '#' devant la ligne suivante pour lancer l'installation :

# install.packages(c("ggplot2", "forecast", "tseries", "GaussianHMM1d", "fGarch", "tidyverse"))

# --- CHARGEMENT ---

# ggplot2 : Indispensable pour créer des graphiques clairs et professionnels.
library(ggplot2)

# forecast : Contient les outils pour les modèles ARIMA/SARIMA et l'autocorrélation.
library(forecast)

# tseries : Contient des outils pour l'analyse statistique et les tests de stationnarité.
library(tseries)

# Note : D'autres librairies seront chargées plus loin

I. Activité 1 : Les séries chronologiques

Cette première partie se concentre sur les étapes exploratoires essentielles des séries temporelles : structuration des données, visualisation et identification des composantes temporelles.

1. Importation et préparation des données

L’analyse commence toujours par l’importation des données. Nous utilisons ici l’exemple des cas de rougeole à Londres, observés de façon bihebdomadaire. Le fichier .RData contient déjà l’objet structuré pour R. 

# On charge le fichier de données (assurez-vous qu'il est dans votre dossier de travail)
load("twentymeas.RData")

# Les données sont contenues dans une liste nommée 'twentymeas'
# On extrait spécifiquement les données de Londres
data_londres <- twentymeas$London

# On vérifie les premières lignes pour comprendre la structure (date, nombre de cas, nombre de naissance, population)
head(data_londres)

Création d’un objet de type “série chronologique”

Pour que R comprenne qu’il s’agit d’une suite temporelle, nous devons transformer les données avec la fonction ts(). 

# On crée la série :
# - data_londres$cases : les valeurs observées
# - start = 1944 : l'année de début
# - frequency = 26 : car les données sont collectées toutes les 2 semaines 
#   (soit 26 périodes par an)

ts_rougeole <- ts(data_londres$cases[1:496],
                  start=data_londres$time[1],
                  frequency=26)

2. Visualisation des données

Comme expliqué dans l’activité 1, la visualisation est la première étape pour détecter des tendances ou des cycles.

Représentation pour données discrètes (barres)

En épidémiologie, les données sont souvent discrètes (un nombre entier de cas par période). La représentation en barres reflète fidèlement cette collecte. 

ggplot(data_londres, aes(x = time, y = cases)) +
  geom_bar(stat = "identity", fill = "#2c7bb6", width = 0.02) +
  labs(title = "Nombre de cas de rougeole (représentation discrète)",
       y = "Nombre de cas", x = "Année") +
  theme_bw()

Représentation en courbe continue

Bien que les données soient discrètes, on peut représenter l’évolution des nombres de cas sous forme de points reliés par une ligne pour mieux visualiser la dynamique général (tendances à long terme, cycles saisonniers) 

ggplot(data_londres, aes(x = time, y = cases)) +
  geom_line(color = "#2c7bb6") +
  labs(title = "Évolution temporelle des cas de rougeole",
       y = "Nombre de cas", x = "Année") +
  theme_bw()

3. Autocorrélation

Comme vu dans les activités, l’autocorrélation mesure comment une valeur passée influence la valeur actuelle, comme une rangée de dominos.

Stabilisation de la variance

Les séries de maladies infectieuses ont souvent des pics très variables, ce qui est typique des modèles multiplicatifs. On applique le logarithme pour réduire ces écarts et stabiliser la série. Cela permet aussi de transformer un modèle multiplicatif en un modèle additif. 

# On applique une transformation logarithmique sur la série temporelle

log_rougeole <- log(data_londres$cases)

ggplot(data_londres, aes(x = time, y = log(cases))) +
  geom_line(color = "#2c7bb6") +
  labs(title = "Évolution temporelle des cas de rougeole",
       y = "Nombre de cas (logarithme)", x = "Année") +
  theme_bw()

On aboutit à une série plus régulière après la transformation.

Fonctions d’autocorrélation et d’autocorrélation partielle (ACF et PACF)

  • ACF : Mesure la corrélation globale entre aujourd’hui et les périodes passées (\(k\) lags).

  • PACF : Isole la relation directe en filtrant les “échos” des périodes intermédiaires.

Les fonctions suivantes nous permettent de retrouver les figures vues dans les activités et les applications Shiny. 

# Graphique de l'autocorrélation
ggAcf(log_rougeole, size = 1.5, color = "#fdae61", lag.max= 20) + 
  labs(title = "Fonction d'autocorrélation (ACF)") + theme_bw()

# Graphique de l'autocorrélation partielle
ggPacf(log_rougeole, size = 1.5, color = "#fdae61", lag.max= 20) + 
  labs(title = "Fonction d'autocorrélation partielle (PACF)") + theme_bw()

5. Décomposition

Décomposer une série, c’est comme ouvrir le mécanisme d’une montre pour observer ses rouages.

On distingue :

  1. La tendance (\(T_t\)) : l’évolution à long terme.

  2. La saisonnalité (\(S_t\)) : les cycles qui se répètent (ex: chaque année).

  3. Le bruit (\(\epsilon_t\)) : les variations aléatoires imprévisibles.

Pour le modèle additif : \(y_t = T_t + S_t + \epsilon_t\). 

ts_log_rougeole <- ts(log(data_londres$cases[1:496]),start=data_londres$time[1],frequency=26)

# On décompose la série transformée (log) avec la fonction decompose
decomp_rougeole <- decompose(ts_log_rougeole, type = "additive")

# Affichage des 4 graphiques : Série originale, Tendance, Saisonnalité et Bruit
plot(decomp_rougeole)

II. Activité 2 : Les modèles de séries chronologiques

Cette partie du guide se concentre sur la modélisation statistique. Comme nous l’avons exploré dans l’activité 2, l’objectif est de transformer nos observations en modèles mathématiques capables de capturer la dépendance temporelle (la mémoire) et de prévoir l’évolution des phénomènes de santé. Nous aborderons ici les modèles ARMA, les tests de stationnarité, ainsi que les approches pour gérer la volatilité et les changements de régimes.

1. Modélisation ARMA (AutoRegressive Moving Average)

Le modèle ARMA est le pilier de l’analyse des séries stationnaires. Il combine deux forces : la partie autorégressive (AR), qui utilise la mémoire des valeurs passées, et la partie moyenne mobile (MA), qui s’appuie sur les erreurs de prévisions précédentes pour corriger le modèle.

Simulation d’un modèle AR(2)

Dans un modèle AR(2), la valeur d’aujourd’hui dépend des deux périodes précédentes. C’est l’idée que les données portent l’écho du passé récent. Pour un modèle AR(p), la fonction d’autocorrélation partielle (PACF) s’annule généralement après le décalage (lag) p. 

# Simulation d'un modèle AR(2) : yt = 0.9yt-1 - 0.5yt-2 + et
set.seed(123)
x1 <- arima.sim(n = 250, list(ar = c(0.9, -0.5)), sd = sqrt(2))

# Conversion en dataframe pour ggplot
df_x1 <- data.frame(Temps = 1:250, Valeur = as.numeric(x1))

# Visualisation de la série
ggplot(df_x1, aes(x = Temps, y = Valeur)) +
  geom_line(color = "#2c7bb6") +
  labs(title = "Simulation d'un modèle AR(2)") +
  theme_bw()

# Visualisation des fonctions ACF et PACF
# Notez comment la PACF chute brusquement après le lag 2
ggAcf(x1 , size = 1.5, color = "#fdae61", lag.max= 20) + 
  labs(title = "ACF d'un processus AR(2)") + theme_bw()

ggPacf(x1, size = 1.5, color = "#fdae61", lag.max= 20) + 
  labs(title = "PACF d'un processus AR(2)") + theme_bw()

Simulation d’un modèle MA(2)

Le modèle à moyenne mobile (MA) repose sur une idée complémentaire : il apprend des erreurs passées pour affiner les estimations. Sa dépendance est dite de “courte portée”. Pour un modèle MA(q), c’est l’ACF qui s’annule après le décalage q. 

# Simulation d'un modèle MA(2) avec coefficients -0.2 et 0.2
x2 <- arima.sim(n = 250, list(ma = c(-0.2, 0.2)), sd = sqrt(2))
df_x2 <- data.frame(Temps = 1:250, Valeur = as.numeric(x2))

# Visualisation
ggplot(df_x2, aes(x = Temps, y = Valeur)) +
  geom_line(color = "#d7191c") +
  labs(title = "Simulation d'un modèle MA(2)") +
  theme_bw()

# Visualisation de l'ACF et de la PACF
# On observe que la dépendance s'estompe rapidement après le lag 2
ggAcf(x2, size = 1.5, color = "#fdae61", lag.max= 20) + 
  labs(title = "ACF d'un processus MA(2)") + theme_bw()

ggPacf(x2, size = 1.5, color = "#fdae61", lag.max= 20) + 
  labs(title = "PACF d'un processus MA(2)") + theme_bw()

Simulation d’un modèle ARMA(2,2)

Le modèle ARMA(2,2) combine les deux dynamiques pour une plus grande flexibilité. Il représente à la fois les effets persistants du passé et les ajustements rapides dus aux fluctuations. 

# Simulation d'un ARMA(2,2)
x3 <- 10 + arima.sim(n = 250, list(ar = c(0.9, -0.5), ma = c(-0.2, 0.2)), sd = sqrt(2))
df_x3 <- data.frame(Temps = 1:250, Valeur = as.numeric(x3))

# Visualisation
ggplot(df_x3, aes(x = Temps, y = Valeur)) +
  geom_line(color = "#1a9641") +
  labs(title = "Simulation d'un modèle ARMA(2,2)", subtitle = "Équilibre entre mémoire et ajustement") +
  theme_bw()

# Diagnostic ACF/PACF
ggAcf(x3, size = 1.5, color = "#fdae61", lag.max= 20) + 
  labs(title = "ACF d'un processus ARMA(2,2)") + theme_bw()

ggPacf(x3, size = 1.5, color = "#fdae61", lag.max= 20) + 
  labs(title = "PACF d'un processus ARMA(2,2)") + theme_bw()

2. Les tests de stationnarité

Comme nous l’avons vu, la stationnarité est une condition essentielle pour garantir la fiabilité des modèles. Pour confirmer nos impressions visuelles, nous utilisons deux tests complémentaires :

  • Le test ADF (Augmented Dickey-Fuller) : Il cherche à détecter une “racine unitaire” (instabilité). Si la p-valeur est inférieure à 0,05, on rejette l’instabilité et on conclut que la série est stationnaire.

  • Le test KPSS : Il fonctionne à l’inverse. Il suppose que la série est stationnaire par défaut. Si la p-valeur est élevée, cela confirme la stationnarité.

On va appliquer ces tests sur plusieurs types de séries chronologiques issues de différents mécanismes.

Simulation des séries

Nous avons vu précédemment qu’on peut générer facilement des séries temporelles avec la fonction arima.sim(). Cela marche pour les modèles plus simples. Si on veut générer des données avec un mechanismes plus compliqué, on peut utiliser d’autres fonction, notamment celles de la librairie GaussianHMM1d qui permet de simuler des séries basées sur des modèles à changement de régime (détails dans la section 5). 

# Modèle AR(2)
set.seed(1)  # pour toujours obtenir les mêmes séries
mod.ar2 = arima.sim(n = 250, list(ar = c(0.9, -0.5)),sd = sqrt(2))


# Modèle HMM gaussien à 2 régimes (détails abordés dans la section 5)
library(GaussianHMM1d)
Q <- matrix(c(0.85, 0.25, 0.15, 0.75),2,2) 
mu <- c(10, 25)
sigma <- c(3, 7)

mod.hmm = Sim.HMM.Gaussian.1d(mu,sigma,Q,eta0=1,n=250)$x


# Marche aléatoire
set.seed(1) 
mod.marcheal = cumsum(0.5+2*rnorm(250))

# Série saisonnière
## Utilisation d'une fonction sinusoïdale pour simuler un cycle annuel (52 semaines) 

set.seed(22)
n <- 250
t <- 1:n
mod.season = 20 + 10 * sin(2 * pi * t / 52) + rnorm(n, sd = 4)


# Série saisonnière avec tendance
set.seed(11)
# Ajout d'une dérive progressive (0.2 * t) au cycle saisonnier 
set.seed(123)
mod.stend <- 20 + 10 * sin(2 * pi * t / 52) + 0.2 * t + rnorm(n, sd = 4)

# Série avec point de rupture
# Simule un changement structurel brutal (ex: introduction d'un vaccin) 
set.seed(23)
x1_rupt <- arima.sim(n = 120, list(ar = 0.7, ma = -0.2), sd = 1)
x2_rupt <- arima.sim(n = 130, list(ar = -0.4, ma = 0.2), sd = 1)
mod.rupt <- c(x1_rupt, x2_rupt + 5) # Le "+5" crée le décalage après la rupture

Visualisation de différents types de série chronologique 

# Organisation des données pour ggplot
# Nous regroupons les 6 séries simulées précédemment dans un seul tableau (format long)
df_comparaison <- data.frame(
  Temps = rep(1:250, 6),
  Valeur = c(mod.ar2, mod.hmm, mod.marcheal, mod.season, mod.stend, mod.rupt),
  Type = rep(c("1. AR(2) Stationnaire", 
               "2. HMM (Régimes)", 
               "3. Marche Aléatoire", 
               "4. Saisonnière", 
               "5. Saisonnière + Tendance", 
               "6. Point de rupture"), 
             each = 250)
)

# Création du graphique en grille (Faceting)
# Cette méthode permet de comparer visuellement la stabilité des moyennes et des variances
ggplot(df_comparaison, aes(x = Temps, y = Valeur, color = Type)) +
  geom_line(alpha = 0.8) +
  facet_wrap(~Type, scales = "free_y", ncol = 3) +
  labs(title = "Comparaison de différentes séries temporelles") +
  theme_bw() +
  theme(legend.position = "none", # On retire la légende, car les titres sont dans les facettes
        strip.background = element_rect(fill = "grey90"),
        strip.text = element_text(face = "bold"))

Stationnarité de différents types de série chronologique

Nous allons automatiser l’application des tests ADF et KPSS sur l’ensemble de nos simulations. L’objectif est de voir d’un seul coup d’œil comment les tests réagissent selon la nature de la série.

Un résultat est robuste lorsque les deux tests s’accordent. En cas de contradiction, des analyses plus approfondies ou une transformation (différenciation) peuvent être nécessaires. 

# Fonction pour extraire les résultats des tests (p-valeurs)
extraire_tests <- function(serie, nom) {
  # Test ADF (H0 : Non stationnaire)
  test_adf <- adf.test(serie)
  
  # Test KPSS (H0 : Stationnaire)
  # null = "Trend" permet de tester la stationnarité autour d'une tendance
  test_kpss <- kpss.test(serie, null = "Trend")
  
  # Retourne une ligne de tableau
  data.frame(
    Série = nom,
    ADF = round(test_adf$p.value, 3),
    KPSS = round(test_kpss$p.value, 3)
  )
}

# Application de la fonction sur nos 6 séries simulées
resultats <- rbind(
  extraire_tests(mod.ar2, "1. AR(2) Stationnaire"),
  extraire_tests(mod.hmm, "2. HMM (Régimes)"),
  extraire_tests(mod.marcheal, "3. Marche Aléatoire"),
  extraire_tests(mod.season, "4. Saisonnière"),
  extraire_tests(mod.stend, "5. Saisonnière + Tendance"),
  extraire_tests(mod.rupt, "6. Point de rupture")
)

# Ajout de l'interprétation automatique pour faciliter la lecture
# Rappel : ADF < 0.05 = Stationnaire | KPSS > 0.05 = Stationnaire

resultats$Interprétation <- ifelse(
  resultats$ADF < 0.05 & resultats$KPSS > 0.05,
  "✅ Stationnaire",
  ifelse(
    resultats$ADF >= 0.05 & resultats$KPSS <= 0.05,
    "❌ Non stationnaire",
    "⚠️ Contradictoire / À vérifier"
  )
)

# Affichage du tableau récapitulatif
# On peut utiliser knitr::kable pour un rendu propre du tableau
knitr::kable(resultats, 
             caption = "Comparaison des tests de stationnarité (ADF et KPSS)")
Comparaison des tests de stationnarité (ADF et KPSS)
Série ADF KPSS Interprétation
1. AR(2) Stationnaire 0.010 0.10 ✅ Stationnaire
2. HMM (Régimes) 0.010 0.10 ✅ Stationnaire
3. Marche Aléatoire 0.532 0.01 ❌ Non stationnaire
4. Saisonnière 0.196 0.10 ⚠️ Contradictoire / À vérifier
5. Saisonnière + Tendance 0.283 0.10 ⚠️ Contradictoire / À vérifier
6. Point de rupture 0.213 0.01 ❌ Non stationnaire

Analyse des résultats

L’observation de ce tableau permet de tirer plusieurs conclusions importantes pour vos analyses :

• La concordance parfaite (Séries 1, 2 et 3) : Pour les modèles AR(2) et HMM, les deux tests s’accordent pour confirmer la stationnarité. À l’inverse, pour la marche aléatoire, l’ADF et le KPSS rejettent tous deux la stationnarité (p-valeur élevée pour l’un, faible pour l’autre). C’est le scénario idéal où la décision est simple et robuste.

• Le paradoxe de la saisonnalité (Séries 4 et 5) : Notez les résultats contradictoires pour les séries saisonnières. Ici, le test ADF ne rejette pas l’instabilité (\(p = 0,196\)), alors que le test KPSS ne rejette pas la stationnarité (\(p = 0,10\)). Les cycles réguliers (comme ceux de la grippe ou de la dengue) peuvent « tromper » les tests selon la longueur de la série ou la force du signal saisonnier. D’autres étapes sont nécessaires (à voir dans la section III)

• L’impact des ruptures (Série 6) : Le point de rupture est clairement identifié comme non stationnaire par les deux tests. Un changement brutal dans la structure des données (par exemple, suite à l’introduction d’un vaccin ou de l’émergence d’un nouveau variant) modifie la moyenne et/ou la variance, rendant la série non stationnaire. D’autres modèles peuvent donc être nécessaires (hors du cadre de ce module).

Stationnarité de la série sur les cas de rougeole

Nous appliquons maintenant les deux tests sur notre série sur les cas de rougeole, transformée en logarithme. 

# Application des tests sur la série des cas de rougeole (transformée en log)
# Ces tests permettent de valider mathématiquement la stabilité de la série
adf_result <- adf.test(log_rougeole)
kpss_result <- kpss.test(log_rougeole, null = "Trend")

# Affichage des résultats
print(adf_result)
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  log_rougeole
#> Dickey-Fuller = -7.7861, Lag order = 8, p-value = 0.01
#> alternative hypothesis: stationary
print(kpss_result)
#> 
#>  KPSS Test for Trend Stationarity
#> 
#> data:  log_rougeole
#> KPSS Trend = 0.031184, Truncation lag parameter = 6, p-value = 0.1

L’analyse des sorties nous donne une confirmation de la stationnarité de notre série de cas de rougeole :

Résultat du test ADF (p-value = 0.01) : La p-valeur est inférieure au seuil de 0,05. Nous rejetons donc l’hypothèse nulle d’instabilité. Selon ce test, la série est stationnaire. Résultat du test KPSS (p-value = 0.1) : La p-valeur est supérieure au seuil de 0,05. Nous ne rejetons donc pas l’hypothèse de stationnarité. Selon ce test, la série est également stationnaire.

Les deux tests s’accordent, nous avons la confirmation que la série des cas de rougeole (à l’échelle logarithmique) est stationnaire.

Cette étape est fondamentale. C’est parce que nous avons validé la stationnarité que nous pouvons maintenant passer à l’étape suivante : l’ajustement d’un modèle ARMA pour analyser la mémoire du passé et tenter de prévoir le nombre de cas au fil du temps.

3. Modélisation de la série de cas de rougeole

Après avoir validé la stationnarité, l’étape suivante consiste à identifier l’ordre \((p, q)\) du processus ARMA22. L’objectif est de trouver le modèle qui explique le mieux les données avec le moins de paramètres possible, respectant ainsi le principe de parcimonie.

Ajustement et sélection des modèles

En pratique, on ne connaît pas à l’avance le “meilleur” modèle. On en teste plusieurs et on compare leur performance à l’aide de critères d’information.

  • AIC (Akaike Information Criterion) : Favorise la qualité de l’ajustement.

  • BIC (Bayesian Information Criterion) : Plus sévère, il pénalise davantage la complexité inutile.

Plus la valeur de l’AIC ou du BIC est faible, meilleur est le modèle.

Pour évaluer la performance de nos modèles, nous divisons notre série de 548 observations en deux parties:

  • Période d’ajustement (1 à 496) : Utilisée pour estimer les paramètres du modèle.

  • Période de validation (497 à 548) : Les 52 dernières périodes (soit environ deux ans) serviront à comparer nos prédictions aux cas réels. 

ind0=c(1:496) ## Indicateur pour sélectionner les 496 premières observations

log_rougeole = log(data_londres$cases)
# Test de différents modèles candidats avec la fonction arima
fit_ar1    <- arima(log_rougeole[ind0], order = c(1, 0, 0))
fit_ma1    <- arima(log_rougeole[ind0], order = c(0, 0, 2))
fit_arma21 <- arima(log_rougeole[ind0], order = c(2, 0, 1))

# Comparaison des critères AIC
# Le modèle avec le AIC le plus bas est généralement préféré
comparaison_aic <- data.frame(
  Modele = c("AR(1)", "MA(1)", "ARMA(2,1)"),
  AIC = c(AIC(fit_ar1), AIC(fit_ma1), AIC(fit_arma21))
)
print(comparaison_aic)
#>      Modele      AIC
#> 1     AR(1) 348.7675
#> 2     MA(1) 822.2037
#> 3 ARMA(2,1) 176.8837

Le modèle ARMA(2,1) est le meilleur parmi les trois.

Automatisation avec auto.arima()

En pratique, explorer manuellement toutes les combinaisons peut être fastidieux. La fonction auto.arima() de la librairie forecast automatise ce processus en parcourant l’espace des modèles pour identifier celui qui minimise les critères d’information. 

# Identification automatique du meilleur modèle
# Cette fonction gère l'estimation et la comparaison de multiples structures
mod.rougeole <- auto.arima(log_rougeole[ind0], seasonal = FALSE)

# Affichage du modèle sélectionné
# Dans notre exemple, le ARMA(2,1) s'est révélé le plus performant
summary(mod.rougeole)
#> Series: log_rougeole[ind0] 
#> ARIMA(2,0,1) with non-zero mean 
#> 
#> Coefficients:
#>          ar1      ar2      ma1    mean
#>       1.8726  -0.9026  -0.6565  6.0149
#> s.e.  0.0298   0.0285   0.0563  0.1465
#> 
#> sigma^2 = 0.08202:  log likelihood = -83.44
#> AIC=176.88   AICc=177.01   BIC=197.92
#> 
#> Training set error measures:
#>                        ME      RMSE      MAE        MPE     MAPE     MASE
#> Training set 0.0009166998 0.2852271 0.226512 -0.2807904 4.243169 0.807111
#>                    ACF1
#> Training set 0.01585857

On retrouve les mêmes résultats tels qu’expliqués dans l’activité 2 du site.

Prédictions : Valeurs prévues du modèle et nouvelles prévisions

Une fois le modèle sélectionné, il est essentiel de vérifier sa capacité à reproduire la réalité. Cette étape de validation se fait en deux temps : d’abord en observant comment le modèle s’ajuste aux données passées (vérification de l’adéquation ou “model fitting”), puis en testant sa capacité à prédire l’avenir (“forecasting”).

Comparaison entre valeurs réelles et valeurs ajustées (1944-1945)

Pour évaluer l’adéquation du modèle, nous calculons les valeurs « ajustées ». Mathématiquement, une observation est la somme de la prédiction du modèle et d’une erreur (le résidu). En soustrayant le résidu de la valeur réelle (\(Observed - Residual\)), on obtient ce que le modèle a “compris” de la série. 

# On calcule les valeurs ajustées (Fitted values)
# On soustrait l'erreur (residuals) à la valeur observée (log_rougeole)
# pour obtenir la partie expliquée par le modèle ARMA(2,1).

ARMA_fit.52 <- log_rougeole[1:52] - residuals(mod.rougeole)[1:52]

# On prépare un tableau de données pour les deux premières années (52 points)
# Cela permet de zoomer sur la dynamique locale.

data1 = data.frame(cases= data_londres$cases[1:52],time=data_londres$time[1:52],fit=ARMA_fit.52)

# Création de la figure (cas réels vs cas prévus) pour deux ans (1944 - 1945)

ggplot(data1, aes( x = time,
                        y = log(cases),
                        color = "Cas réels")) + 
  geom_line() +
  geom_point() + 
  ggtitle("") +
  labs(y = "Nombre de cas (log) ", x = "Année") + geom_line(aes(y = fit, color = "Cas prévus")) +
  scale_color_manual(name = "",
                     values = c("Cas réels" = "#67a9cf","Cas prévus" = "#ef8a62"))+
  theme_bw()

On remarque que les valeurs prévues semblent toujours en retard par rapport aux observations. Ceci est normal, car les valeurs prévues sont basées sur les valeurs passées. Par contre, les valeurs prévues sont très proches des valeurs observées.

Prédiction de valeurs futures avec la fonction predict

Pour la suite, on peut aussi prévoir le nombre de cas de rougeole (log) pour 1963-1964, en se servant du modèle estimé pour les années antérieures. Les erreurs de prévision sont plus grandes, comme on peut s’y attendre. 

# Génération de la prédiction avec la fonction predict()
# n.ahead = 52 indique que nous voulons prévoir les 52 prochaines semaines (2 ans)
pred.rougeole <- predict(mod.rougeole, n.ahead = 52)

# Préparation du tableau de données pour la période de validation
# On inclut les bornes de confiance (Intervalle à 95%)

data2 = data.frame(cases=data_londres$cases[497:548], 
                   time=data_londres$time[497:548],
                   pred = pred.rougeole$pred, se = pred.rougeole$se, 
                   lb = pred.rougeole$pred - 2 * pred.rougeole$se, 
                   ub = pred.rougeole$pred + 2*pred.rougeole$se)

# Création de la figure 
ggplot(data2, aes(x=time, y=log(cases), color = "Cas réels")) +
  geom_line()+
  geom_point()+
  ggtitle("") +
  labs(y = "Nombre de cas (log) ", x = "Année")+
  geom_line(aes(y = pred, color = "Cas prévus"))+
  scale_color_manual(name = "", values = c("Cas réels" = "#67a9cf", "Cas prévus" = "#ef8a62"))+
  theme_bw()

4. Modèles de volatilité (ARCH et GARCH)

Certaines séries présentent des périodes calmes suivies de brusques “orages” de variabilité (pics de volatilité), comme c’est souvent le cas lors des flambées épidémiques. Les modèles ARCH et GARCH permettent de modéliser cette variance qui change dans le temps. 

# Nous utilisons la librairie fGarch pour simuler un processus GARCH(1,1)
library(fGarch)
set.seed(10) # pour avoir des resultats reproductibles

# Définition des spécifications du modèle GARCH(1,1)
# omega : variance de base
# alpha : réaction aux chocs récents (effet ARCH)
# beta  : persistance de la volatilité dans le temps (effet GARCH)
x4_garch <- garchSpec(model = list(omega = 2, alpha = 0.05, beta = 0.94))

# simulation de la série
x4_sim <- garchSim(x4_garch, n = 250)
df_x4 <- data.frame(Temps = 1:250, Valeur = as.numeric(x4_sim))

# Visualisation de la série
ggplot(df_x4, aes(x = Temps, y = Valeur)) +
  geom_line() +
  labs(title = "Série avec volatilité GARCH", subtitle = "Notez l'alternance entre calme et instabilité") +
  theme_bw()

# Diagnostic ACF/PACF
library(patchwork)

ggAcf(x4_sim, size = 1.5, color = "#fdae61", lag.max= 20) + 
  labs(title = "ACF d'un processus GARCH(1,1)") + theme_bw() | 
ggPacf(x4_sim, size = 1.5, color = "#fdae61", lag.max= 20) + 
  labs(title = "PACF d'un processus GARCH(1,1)") + theme_bw()

# Test de stationnarité
adf.test(x4_sim)
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  x4_sim
#> Dickey-Fuller = -6.0332, Lag order = 6, p-value = 0.01
#> alternative hypothesis: stationary
kpss.test(x4_sim, null = "Trend")
#> 
#>  KPSS Test for Trend Stationarity
#> 
#> data:  x4_sim
#> KPSS Trend = 0.085092, Truncation lag parameter = 5, p-value = 0.1

5. Modèles à changement de régimes (HMM)

Le modèle de Markov caché (Hidden Markov Model, HMM) est particulièrement puissant en épidémiologie. Il suppose que le système bascule entre plusieurs états “cachés” (états latents), par exemple entre une phase endémique (calme) et une phase épidémique (active). 

library(GaussianHMM1d)

# Paramètres : Régime 1 (Calme) vs Régime 2 (Épidémique)
Q <- matrix(c(0.8, 0.3, 0.2, 0.7), 2, 2) # Matrice de transition
mu <- c(-0.3, 0.7) ; sigma <- c(0.2, 0.6)

# Simulation
x5_hmm <- Sim.HMM.Gaussian.1d(mu, sigma, Q, eta0=1, n = 250)$x
df_x5 <- data.frame(Temps = 1:250, Valeur = x5_hmm)


fit_hmm <- EstHMM1d(x5_hmm, reg = 2)
regimes_estimes <- EstRegime(1:250, x5_hmm, fit_hmm$lambda, fit_hmm$eta)$B

df_x5$Regime <- as.factor(regimes_estimes)
# Visualisation de la série
hmm_vis <- ggplot(df_x5, aes(x = Temps, y = Valeur)) +
  geom_line() +
  labs(title = "Série avec changement de régime", subtitle = "Notez l'alternance des phases") +
  theme_bw()

# Visualisation avec identification des régimes (Modélisation)
hmm_regime <- ggplot(df_x5, aes(x = Temps, y = Valeur, color = Regime)) +
  geom_point() +
  scale_color_manual(values = c("black", "red"), labels = c("Calme", "Épidémique")) +
  labs(title = "Série HMM : Détection des phases épidémiques", subtitle = "Les points rouges indiquent le passage au régime actif") +
  theme_bw()

library(patchwork)
hmm_vis/hmm_regime

Dans le cas HMM, les régimes ne sont pas observables, mais ils peuvent être prévus. Dans cet exemple, le régime 2 (en rouge) est le régime ayant la plus petite moyenne.

III. Activité 3 : Les covariables et les effets spatiaux

Cette dernière partie du guide traite des séries non stationnaires, soumises à des tendances ou à des cycles saisonniers. Nous allons apprendre à stabiliser ces séries par la différenciation et à intégrer des composantes saisonnières ainsi que des facteurs externes (covariables et espace).

1. Les modèles ARIMA : Gérer la non stationnarité

Comme nous l’avons vu, la plupart des séries épidémiologiques peuvent présenter des tendances à long terme, ce qui rend la série non stationnaire. Le modèle ARIMA (AutoRegressive Integrated Moving Average) résout ce problème par la différenciation.

Simulation d’un modèle ARIMA(2,1,2)

# Simulation d'un processus ARIMA(2,1,2) :
# n = 250 : Génère 250 points d'observation.
# order = c(2, 1, 2) : Définit la structure du modèle (p, d, q).
#   - p = 2 (AR) : Dépendance aux 2 valeurs précédentes.
#   - d = 1 (I)  : le modèle intègre une différenciation de premier ordre pour simuler une tendance.
#   - q = 2 (MA) : Dépendance aux 2 erreurs de prévision précédentes.
# ar = c(0.7, -0.2) : Définit la force de la mémoire du passé (coefficients AR).
# ma = c(0.5, 0.3)  : Définit la correction basée sur les chocs passés (coefficients MA).


set.seed(123)
x_arima <- arima.sim(n = 250, list(order = c(2, 1, 2), 
                                   ar = c(0.7, -0.2), 
                                   ma = c(0.5, 0.3)))

x_arima_ts <- ts(x_arima, frequency = 26)

# Visualisation et décomposition de la série simulée 
plot(decompose(x_arima_ts))

Un examen de la fonction d’autocorrélation (ACF) et l’application de tests statistiques permettent de confirmer la nécessité d’une transformation. 

# Observation de l'ACF de la série brute
# Un ACF qui décroît très lentement confirme le besoin de différenciation 
ggAcf(x_arima, size = 1.5, color = "#fdae61", lag.max= 20) + 
  labs(title = "ACF d'une série ARIMA(2,1,2)") + theme_bw() |
ggPacf(x_arima, size = 1.5, color = "#fdae61", lag.max= 20) + 
  labs(title = "PACF d'une série ARIMA(2,1,2)") + theme_bw()

# test de stationnarité
adf.arima <- adf.test(x_arima)
kpss.arima <- kpss.test(ts(x_arima), null = "Trend")

# Affichage des résultats
print(adf.arima)
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  x_arima
#> Dickey-Fuller = -2.5124, Lag order = 6, p-value = 0.36
#> alternative hypothesis: stationary
print(kpss.arima)
#> 
#>  KPSS Test for Trend Stationarity
#> 
#> data:  ts(x_arima)
#> KPSS Trend = 0.26103, Truncation lag parameter = 5, p-value = 0.01
## Les résultats des deux tests confirment la non stationnarité

Le principe de la différenciation La différenciation consiste à calculer la variation entre deux moments successifs (\(\nabla y_t = y_t - y_{t-1}\)). Cela permet de supprimer une tendance linéaire pour rendre la série stationnaire. On va appliquer une différenciation à cette série simulée et voir le comportement de la série différenciée. 

# Application de la différenciation avec la fonctin diff()
diff_arima <- diff(x_arima)
diff_arima_ts <- ts(diff_arima, frequency = 26)

# Visualisation et décomposition de la série après différenciation
plot(decompose(diff_arima_ts))

# Observation de l'ACF de la série après différenciation
ggAcf(diff_arima, size = 1.5, color = "#fdae61", lag.max= 20) + 
  labs(title = "ACF après différenciation") + theme_bw() |
ggPacf(diff_arima, size = 1.5, color = "#fdae61", lag.max= 20) + 
  labs(title = "PACF après différenciation") + theme_bw()

# test de stationnarité
adf.diff_arima <- adf.test(diff_arima)
kpss.diff_arima <- kpss.test(ts(diff_arima), null = "Trend")

# Affichage des résultats
print(adf.diff_arima)
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  diff_arima
#> Dickey-Fuller = -5.1024, Lag order = 6, p-value = 0.01
#> alternative hypothesis: stationary
print(kpss.diff_arima)
#> 
#>  KPSS Test for Trend Stationarity
#> 
#> data:  ts(diff_arima)
#> KPSS Trend = 0.038785, Truncation lag parameter = 5, p-value = 0.1

Nous pouvons voir qu’après une première différenciation, la série est devenue stationnaire.

Application aux données de la rougeole

Nous appliquons maintenant ces concepts à la série chronologique des cas de rougeole à Londres afin de déterminer si une différenciation est requise pour optimiser la modélisation.

Effet de la différenciation sur la série de rougeole 

# Visualisation de la série de rougeole différenciée (logarithme)
# La fonction diff() applique yt - yt-1
log_rougeole_diff <- diff(log_rougeole[ind0])

df_diff <- data.frame(
  Temps = data_londres$time[ind0][-1], # On enlève le premier point à cause de la différence
  Valeur = as.numeric(log_rougeole_diff)
)

# Comparaison visuelle : Série brute vs Série différenciée

p_brute <- ggplot(data_londres, aes(x = time, y = log(cases))) +
  geom_line(color = "#2c7bb6") +
  labs(title = "Évolution temporelle des cas de rougeole",
       y = "Nombre de cas (logarithme)", x = "Année") +
  theme_bw()

p_diff <- ggplot(df_diff, aes(x = Temps, y = Valeur)) +
  geom_line(color = "#2c7bb6") +
  labs(title = "Série de la rougeole après une différenciation (d=1)",
       subtitle = "La différenciation stabilise la série avec des fluctuations plus régulières",
       y = "Variation du log(Cas)", x = "Année") +
  theme_bw()

library(patchwork)
p_brute / p_diff

Comparaison de modèles pour les données de rougeole : ARMA vs ARIMA

Pour sélectionner le modèle le plus performant, nous comparons le modèle ARMA(2,1) et le modèle ARIMA(2,1,1) (incluant une différenciation). Le critère de sélection est l’AIC (Akaike Information Criterion); la valeur la plus faible indique le meilleur compromis entre ajustement et parcimonie. 

# 1. Modèle ARMA(2,1) sur les données en log (Stationnaire par hypothèse)
arma_rougeole <- arima(log_rougeole[ind0], order = c(2, 0, 1))

# 2. Modèle ARIMA(2,1,1) sur les données en log (inclus une différenciation d=1)
arima_rougeole <- arima(log_rougeole[ind0], order = c(2, 1, 1))

# Comparaison des critères AIC (le plus faible est le meilleur)
cat("AIC du modèle ARMA(2,1) :", AIC(arma_rougeole), "\n")
#> AIC du modèle ARMA(2,1) : 176.8837
cat("AIC du modèle ARIMA(2,1,1) :", AIC(arima_rougeole), "\n")
#> AIC du modèle ARIMA(2,1,1) : 220.5567

Bien que le modèle ARIMA traite explicitement la tendance par différenciation, nous constatons ici que l’AIC du modèle ARMA(2,1) est inférieur à celui du modèle ARIMA(2,1,1). Cela indique que pour cette série spécifique, la transformation logarithmique initiale a suffi à stabiliser la série, rendant la différenciation supplémentaire superflue. Ce constat corrobore les résultats obtenus précédemment via la fonction auto.arima vu dans la partie II, section 3.1.

2. Le modèle SARIMA : Intégrer la saisonnalité

Alors que les modèles ARIMA corrigent les tendances à long terme, les modèles SARIMA (Seasonal AutoRegressive Integrated Moving Average) permettent de capturer les cycles périodiques qui se répètent à intervalles réguliers. En santé publique, ces modèles sont essentiels pour représenter les maladies saisonnières comme la grippe, la bronchiolite ou la dengue.

Le modèle est noté \((p, d, q)(P, D, Q)_s\), où les lettres majuscules représentent la composante saisonnière et \(s\) la longueur du cycle (par exemple, \(s=26\) pour des données bihebdomadaires).

Simulation d’un modèle saisonnier

Pour comprendre le comportement d’une série saisonnière, nous simulons un processus intégrant un cycle annuel de 26 périodes (bihebdomadaire).

Nous allons utiliser la fonction sarima.sim() de la librairie astsa. Cette fonction permet de spécifier les paramètres \((p, d, q)\) et \((P, D, Q)_s\) de manière très proche de la notation mathématique. 

# --- Installation et chargement ---
# install.packages("astsa") # À faire une seule fois
library(astsa)


# Simulation d'un processus SARIMA(1,1,1)(1,1,1)[26]
# ar/ma : parties non saisonnières | sar/sma : parties saisonnières
# d/D : ordres de différenciation | S : période saisonnière

set.seed(123)
n <- 250

x_sarima <- sarima.sim(ar = 0.7, d = 1, ma = -0.4, 
                       sar = 0.8, D = 1, sma = -0.5, 
                       S = 26, n = n)

# Conversion en objet Time Series
x_sarima_ts <- ts(as.numeric(x_sarima), frequency = 26)

# --- Visualisation et Décomposition ---
plot(decompose(x_sarima_ts))

# --- Diagnostic de la structure ---
# L'ACF doit montrer des pics aux lags saisonniers (26, 52...)
ggAcf(x_sarima_ts, color = "#fdae61", lag.max = 60) + 
  labs(title = "ACF de la série SARIMA simulée",
       subtitle = "Pics réguliers significatifs confirmant le cycle de 26 périodes") + 
  theme_bw() |
ggPacf(x_sarima_ts, color = "#fdae61", lag.max = 60) + 
  labs(title = "PACF de la série SARIMA simulée") + 
  theme_bw()

La présence d’un cycle se manifeste par des pics récurrents dans la fonction d’autocorrélation (ACF) aux multiples de la fréquence \(s\).

Les résultats des deux tests confirment par la suite la non stationnarité 

# test de stationnarité
adf.sarima <- adf.test(x_sarima)
kpss.sarima <- kpss.test(ts(x_sarima), null = "Trend")

# Affichage des résultats
print(adf.sarima)
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  x_sarima
#> Dickey-Fuller = -1.4802, Lag order = 6, p-value = 0.7946
#> alternative hypothesis: stationary
print(kpss.sarima)
#> 
#>  KPSS Test for Trend Stationarity
#> 
#> data:  ts(x_sarima)
#> KPSS Trend = 0.67209, Truncation lag parameter = 5, p-value = 0.01

Application d’une différenciation

L’application d’une différenciation saisonnière permet de rendre la série plus stable. Toutefois, on peut remarquer une tendance encore importante. 

# On soustrait la valeur au temps t par celle au temps t-26
diff_seasonal <- diff(x_sarima, lag = 26)

# Visualisation après stabilisation saisonnière
plot(decompose(ts(diff_seasonal, frequency = 26)))

On peut par la suite appliquer une différenciation non saisonnière supplémentaire. 

diff_sarima <- diff(diff_seasonal)

# Visualisation après stabilisation saisonnière
plot(decompose(ts(diff_sarima, frequency = 26)))

# test de stationnarité
adf.diff_sarima  <- adf.test(diff_sarima )
kpss.diff_sarima  <- kpss.test(ts(diff_sarima ), null = "Trend")

# Affichage des résultats
print(adf.diff_sarima )
#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  diff_sarima
#> Dickey-Fuller = -6.4353, Lag order = 6, p-value = 0.01
#> alternative hypothesis: stationary
print(kpss.diff_sarima )
#> 
#>  KPSS Test for Trend Stationarity
#> 
#> data:  ts(diff_sarima)
#> KPSS Trend = 0.042732, Truncation lag parameter = 4, p-value = 0.1

La série devient encore plus stable et la stationnarité est confirmée par les tests.

Application aux données de la rougeole

Nous allons maintenant déterminer si l’ajout d’une composante saisonnière améliore la modélisation de la rougeole à Londres, dont le cycle épidémique est fortement marqué. 

sarima_rougeole <- Arima(log_rougeole[ind0], order = c(1,1,1), 
                         seasonal = list(order = c(1,1,1),period = 26))

# Comparaison des critères AIC (le plus faible est le meilleur)
cat("AIC du modèle ARMA(2,1) :", AIC(arma_rougeole), "\n")
#> AIC du modèle ARMA(2,1) : 176.8837
cat("AIC du modèle ARIMA(2,1,1) :", AIC(arima_rougeole), "\n")
#> AIC du modèle ARIMA(2,1,1) : 220.5567
cat("AIC du modèle SARIMA(1,1,1)(1,1,1)26 :", AIC(sarima_rougeole), "\n")
#> AIC du modèle SARIMA(1,1,1)(1,1,1)26 : 103.7598

Le modèle SARIMA est donc le modèle le plus adapté pour notre série de cas de rougeole.

Validation et diagnostic final

Après avoir identifié le modèle optimal, l’étape de validation est cruciale pour s’assurer que toute l’information temporelle a été correctement capturée par les paramètres du modèle. Un modèle est considéré comme robuste si ses résidus ne présentent plus aucune structure temporelle résiduelle et se comportent comme un « bruit blanc ».

La fonction checkresiduals() est un outil puissant qui nous permet d’effectuer un diagnostic complet en une seule commande. 

# Vérification des résidus du modèle SARIMA
# Les pics de l'ACF des résidus doivent rester dans les limites de significativité
checkresiduals(sarima_rougeole)

#> 
#>  Ljung-Box test
#> 
#> data:  Residuals from ARIMA(1,1,1)(1,1,1)[26]
#> Q* = 9.2583, df = 6, p-value = 0.1596
#> 
#> Model df: 4.   Total lags used: 10

Interprétation des sorties de la fonction checkresiduals()

Test statistique (Ljung-Box) : Ce test sert à vérifier l’autocorrélation des résidus. Avec une \(p\text{-value}\) de 0,1596 (\(> 0,05\)), nous ne rejetons donc pas l’hypothèse nulle d’absence d’autocorrélation des résidus. Les résidus sont considérés comme un bruit blanc (non autocorrélés). Ce qui valide la qualité de notre modèle.

Diagnostics visuels : Les résidus oscillent aléatoirement autour de zéro (pas de structure apparente). La quasi-totalité des pics de l’ACF se situe dans les limites de significativité, confirmant l’absence de corrélation résiduelle. Par rapport à l’histogramme, les erreurs suivent approximativement une loi normale (courbe en cloche).

3. Covariables et effets spatiaux

Comme nous l’avons exploré dans l’Activité 3, les phénomènes de santé publique ne dépendent pas uniquement de leur propre passé. Ils sont influencés par leur contexte environnemental et leur localisation géographique.

L’intégration de facteurs externes : Les covariables

Les phénomènes de santé sont souvent influencés par des covariables:

Dépendantes du temps : température ou pluviométrie (reflétant des conditions changeantes).

Indépendantes du temps : caractéristiques fixes comme la densité de population.

Exemple de mise en œuvre (à titre indicatif) : Dans R, ces variables externes sont intégrées via l’argument xreg (external regressors) des fonctions de modélisation. 

# Note : Ce code nécessite des séries de données externes (ex: relevés météo).
# On ajuste ici un modèle ARIMA en incluant une covariable de température.

# fit_cov <- auto.arima(l_rougeole_ts, xreg = temperature_hebdo)
# summary(fit_cov)

La dimension géographique : Les effets spatiaux

Le modèle STAR (Spatio-Temporal AutoRegressive) reconnaît que l’évolution d’une série de cas dans une région dépend de son propre passé et de la situation dans les régions voisines.

On parle de pression épidémique externe pour décrire l’influence exercée par le voisinage immédiat sur une localité donnée.

La modélisation spatio-temporelle permet par exemple de générer des cartes de prévision pour visualiser les zones de forte incidence et anticiper la propagation régionale.

Note: La génération de cartes épidémiologiques (comme l’exemple ci-dessus pour l’État de New York) nécessite des fichiers géographiques complexes et des matrices de proximité. Bien que ces outils soient essentiels pour la planification sanitaire, leur mise en œuvre avancée dépasse le cadre de ce module d’introduction.

Félicitations ! Vous avez maintenant parcouru l’ensemble des étapes de modélisation, de l’exploration initiale à la prise en compte de la dimension spatiale. Ces outils R constituent la base technique pour une surveillance épidémiologique moderne et rigoureuse.