Les données de l’incidence de la varicelle sont disponibles du site +Web du Réseau Sentinelles. Nous +les récupérons sous forme d’un fichier en format CSV dont chaque ligne +correspond à une semaine de la période demandée. Nous téléchargeons +toujours le jeu de données complet, qui commence en 1991 et se termine +avec une semaine récente. L’URL est:
+data_url = "http://www.sentiweb.fr/datasets/incidence-PAY-7.csv"
+Voici l’explication des colonnes donnée sur le sur le site +d’origine:
+Nom de colonne | +Libellé de colonne | +
---|---|
week |
+Semaine calendaire (ISO 8601) | +
indicator |
+Code de l’indicateur de surveillance | +
inc |
+Estimation de l’incidence de consultations en nombre de cas | +
inc_low |
+Estimation de la borne inférieure de l’IC95% du nombre de cas de +consultation | +
inc_up |
+Estimation de la borne supérieure de l’IC95% du nombre de cas de +consultation | +
inc100 |
+Estimation du taux d’incidence du nombre de cas de consultation (en +cas pour 100,000 habitants) | +
inc100_low |
+Estimation de la borne inférieure de l’IC95% du taux d’incidence du +nombre de cas de consultation (en cas pour 100,000 habitants) | +
inc100_up |
+Estimation de la borne supérieure de l’IC95% du taux d’incidence du +nombre de cas de consultation (en cas pour 100,000 habitants) | +
geo_insee |
+Code de la zone géographique concernée (Code INSEE) http://www.insee.fr/fr/methodes/nomenclatures/cog/ | +
geo_name |
+Libellé de la zone géographique (ce libellé peut être modifié sans +préavis) | +
La première ligne du fichier CSV est un commentaire, que nous
+ignorons en précisant skip=1
.
Si le fichier n’existe pas déjà sur l’ordinateur, on le télécharge +depuis la source indiquée plus haut. Cette vérification nous évite des +téléchargements répétés et inutiles, et nous permet de travailler sans +connexion internet (si le fichier existe).
+library(fs)
+
+file_name = "incidence_varicelle.csv"
+if (!file_exists(file_name)) {
+ download.file(url=data_url, destfile=file_name, method="auto")
+}
+data = read.csv(file_name, skip=1)
+Regardons ce que nous avons obtenu:
+head(data)
+## week indicator inc inc_low inc_up inc100 inc100_low inc100_up geo_insee
+## 1 202245 7 3895 1762 6028 6 3 9 FR
+## 2 202244 7 4271 2231 6311 6 3 9 FR
+## 3 202243 7 5863 3302 8424 9 5 13 FR
+## 4 202242 7 3770 1950 5590 6 3 9 FR
+## 5 202241 7 4177 2219 6135 6 3 9 FR
+## 6 202240 7 4883 1472 8294 7 2 12 FR
+## geo_name
+## 1 France
+## 2 France
+## 3 France
+## 4 France
+## 5 France
+## 6 France
+tail(data)
+## week indicator inc inc_low inc_up inc100 inc100_low inc100_up
+## 1662 199102 7 16277 11046 21508 29 20 38
+## 1663 199101 7 15565 10271 20859 27 18 36
+## 1664 199052 7 19375 13295 25455 34 23 45
+## 1665 199051 7 19080 13807 24353 34 25 43
+## 1666 199050 7 11079 6660 15498 20 12 28
+## 1667 199049 7 1143 0 2610 2 0 5
+## geo_insee geo_name
+## 1662 FR France
+## 1663 FR France
+## 1664 FR France
+## 1665 FR France
+## 1666 FR France
+## 1667 FR France
+Y a-t-il des points manquants dans nos données ?
+na_records = apply(data, 1, function (x) any(is.na(x)))
+data[na_records,]
+## [1] week indicator inc inc_low inc_up inc100
+## [7] inc100_low inc100_up geo_insee geo_name
+## <0 lignes> (ou 'row.names' de longueur nulle)
+Les deux colonnes qui nous intéressent sont week
et
+inc
. Vérifions leurs classes:
class(data$week)
+## [1] "integer"
+class(data$inc)
+## [1] "integer"
+Ce sont des entiers, tout va bien !
+La gestion des dates est toujours un sujet délicat. Il y a un grand
+nombre de conventions différentes qu’il ne faut pas confondre. Notre
+jeux de données utilise un format que peu de logiciels savent traiter:
+les semaines en format ISO-8601. En
+R
, il est géré par la bibliothèque parsedate:
library(parsedate)
+Pour faciliter le traitement suivant, nous remplaçons ces semaines +par les dates qui correspondent aux lundis. Voici une petite fonction +qui fait la conversion pour une seule valeur:
+convert_week = function(w) {
+ ws = paste(w)
+ iso = paste0(substring(ws, 1, 4), "-W", substring(ws, 5, 6))
+ as.character(parse_iso_8601(iso))
+}
+Nous appliquons cette fonction à tous les points, créant une nouvelle
+colonne date
dans notre jeu de données:
data$date = as.Date(convert_week(data$week))
+Vérifions qu’elle est de classe Date
:
class(data$date)
+## [1] "Date"
+Les points sont dans l’ordre chronologique inverse, il est donc utile +de les trier:
+data = data[order(data$date),]
+C’est l’occasion pour faire une vérification: nos dates doivent être +séparées d’exactement sept jours:
+all(diff(data$date) == 7)
+## [1] TRUE
+Regardons enfin à quoi ressemblent nos données !
+plot(data$date, data$inc, type="l", xlab="Date", ylab="Incidence hebdomadaire")
+Un zoom sur les dernières années montre mieux la localisation des +pics en hiver. Le creux des incidences se trouve en été.
+with(tail(data, 200), plot(date, inc, type="l", xlab="Date", ylab="Incidence hebdomadaire"))
+Étant donné que le pic de l’épidémie se situe en hiver, à cheval
+entre deux années civiles, nous définissons la période de référence
+entre deux minima de l’incidence, du 1er août de l’année \(N\) au 1er août de l’année \(N+1\). Nous mettons l’année \(N+1\) comme étiquette sur cette année
+décalée, car le pic de l’épidémie est toujours au début de l’année \(N+1\). Comme l’incidence de syndrome
+grippal est très faible en été, cette modification ne risque pas de
+fausser nos conclusions. L’argument na.rm=True
dans la
+sommation précise qu’il faut supprimer les points manquants. Ce choix
+est raisonnable car il n’y a qu’un seul point manquant, dont l’impact ne
+peut pas être très fort.
pic_annuel = function(annee) {
+ debut = paste0(annee-1,"-09-01")
+ fin = paste0(annee,"-09-01")
+ semaines = data$date > debut & data$date <= fin
+ sum(data$inc[semaines], na.rm=TRUE)
+ }
+Nous devons aussi faire attention aux premières et dernières années +de notre jeux de données. Les données commencent en décembre 1990, ce +qui ne permet pas de quantifier complètement le pic attribué à 1991. +Nous l’enlevons donc de notre analyse. Par contre, pour une exécution en +novembre 2022, les données se terminent après le 1er septembre 2018, ce +qui nous permet d’inclure cette année.
+annees = 1991:2022
+Nous créons un nouveau jeu de données pour l’incidence annuelle, en
+applicant la fonction pic_annuel
à chaque année:
inc_annuelle = data.frame(annee = annees,
+ incidence = sapply(annees, pic_annuel))
+head(inc_annuelle)
+## annee incidence
+## 1 1991 553895
+## 2 1992 834935
+## 3 1993 642921
+## 4 1994 662750
+## 5 1995 651333
+## 6 1996 564994
+Voici les incidences annuelles en graphique:
+plot(inc_annuelle, type="p", xlab="Année", ylab="Incidence annuelle")
+Une liste triée par ordre décroissant d’incidence annuelle permet de +plus facilement repérer les valeurs les plus élevées:
+incidence_annuelle_triee = inc_annuelle[order(-inc_annuelle$incidence),]
+head(incidence_annuelle_triee)
+## annee incidence
+## 19 2009 841233
+## 2 1992 834935
+## 20 2010 834077
+## 26 2016 779816
+## 14 2004 778914
+## 13 2003 760765
+De la même manière, nous pouvons trouver l’année avec l’incidence la +plus faible:
+tail(incidence_annuelle_triee)
+## annee incidence
+## 1 1991 553895
+## 27 2017 552906
+## 28 2018 539765
+## 12 2002 515343
+## 31 2021 377933
+## 30 2020 221183
+Enfin, un histogramme semble indiquer que les épidémies sont en +majorité regroupées autour d’incidences “moyennes” entre 600 000 et 650 +000 personnes.
+hist(inc_annuelle$incidence, breaks=10, xlab="Incidence annuelle", ylab="Nb d'observations", main="")
+