From 2582c3dae79b66faf0920cd5da303b2b3f764833 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Julien=20Qu=C3=A9non?= Date: Tue, 18 May 2021 23:31:03 +0200 Subject: [PATCH] Exercice incidence varicelle --- module3/exo2/exercice_fr.Rmd | 183 ++++++++++++++++++++++++++++++++--- 1 file changed, 167 insertions(+), 16 deletions(-) diff --git a/module3/exo2/exercice_fr.Rmd b/module3/exo2/exercice_fr.Rmd index 7eece5e..0f8b488 100644 --- a/module3/exo2/exercice_fr.Rmd +++ b/module3/exo2/exercice_fr.Rmd @@ -1,33 +1,184 @@ --- -title: "Votre titre" -author: "Votre nom" -date: "La date du jour" -output: html_document +title: "Analyse de l'incidence de la varicelle" +author: "Julien Quénon" +date: "18/05/2021" +output: + pdf_document: + toc: true + html_document: + toc: true + theme: journal +documentclass: article +classoption: a4paper +header-includes: +- \usepackage[french]{babel} +- \usepackage[upright]{fourier} +- \hypersetup{colorlinks=true,pagebackref=true} --- - ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` -## Quelques explications +### Préparation des données + +Les données de l'incidence de la varicelle sont disponibles du site Web du [Réseau Sentinelles](http://www.sentiweb.fr/). 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 1984 et se termine avec une semaine récente. L'URL est: +```{r} +data_url = "https://www.sentiweb.fr/datasets/incidence-PAY-7.csv" +``` + +On effectue une copie locale (appelée `data_local`des données, accessibles depuis l'adresse `data_url`: +```{r} +data_local = "data_varicelle.csv" +``` +Une fois cette copie locale réalisée, on ne veut plus avoir à télécharger le fichier à chaque exécution du script, car on pourrait rencontrer des problèmes (écrasement des données et/ou url invalide). +Pour éviter cela, on vérifie si la copie locale existe (est non vide). Si oui,on n'a rien à faire et si non, on a besoin de télécharger les données. On utilise pour cela la fonction `!file.exists` (qui vérifie que le [fichier n'existe pas](https://statisticsglobe.com/check-existence-of-local-file-in-r)) +```{r} +if (!file.exists(data_local)){ # si le fichier data_local n'existe PAS + download.file(data_url, data_local, method = "auto") # alors on le télécharge à l'adresse data_url et on le stocke dans data_local +} +``` + +Voici l'explication des colonnes donnée sur le [sur le site d'origine](https://ns.sentiweb.fr/incidence/csv-schema-v1.json): + +| 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) | + +### Lecture des données +On lit les données contenues dans la copie **locale** des données (objet `data_local`): +```{r} +data = read.csv(data_local, skip = 1) +``` +La première ligne du fichier CSV est un commentaire, que nous ignorons en précisant `skip=1`. + +Regardons ce que nous avons obtenu: +```{r} +head(data) +tail(data) +``` +On a des données qui vont de la 49ème semaine de l'année 1990 à la 18ème semaine de l'année 2021 + +Y a-t-il des points manquants dans nos données ? +```{r} +na_records = apply(data, 1, function (x) any(is.na(x))) +data[na_records,] +``` +On n'a pas de données manquantes -Ceci est un document R markdown que vous pouvez aisément exporter au format HTML, PDF, et MS Word. Pour plus de détails sur R Markdown consultez . +Les deux colonnes qui nous intéressent sont `week` et `inc`. Vérifions leurs classes: +```{r} +class(data$week) +class(data$inc) +``` +Ce sont des entiers, tout va bien ! + +### Conversion des numéros de semaine -Lorsque vous cliquerez sur le bouton **Knit** ce document sera compilé afin de ré-exécuter le code R et d'inclure les résultats dans un document final. Comme nous vous l'avons montré dans la vidéo, on inclue du code R de la façon suivante: +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](https://en.wikipedia.org/wiki/ISO_8601). En `R`, il est géré par la bibliothèque [parsedate](https://cran.r-project.org/package=parsedate): -```{r cars} -summary(cars) +```{r} +library(parsedate) ``` -Et on peut aussi aisément inclure des figures. Par exemple: +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: -```{r pressure, echo=FALSE} -plot(pressure) +```{r} +convert_week = function(w) { + ws = paste(w) + iso = paste0(substring(ws, 1, 4), "-W", substring(ws, 5, 6)) + as.character(parse_iso_8601(iso)) +} ``` -Vous remarquerez le paramètre `echo = FALSE` qui indique que le code ne doit pas apparaître dans la version finale du document. Nous vous recommandons dans le cadre de ce MOOC de ne pas utiliser ce paramètre car l'objectif est que vos analyses de données soient parfaitement transparentes pour être reproductibles. +Nous appliquons cette fonction à tous les points, créant une nouvelle colonne `date` dans notre jeu de données: +```{r} +data$date = as.Date(convert_week(data$week)) +``` -Comme les résultats ne sont pas stockés dans les fichiers Rmd, pour faciliter la relecture de vos analyses par d'autres personnes, vous aurez donc intérêt à générer un HTML ou un PDF et à le commiter. +Vérifions qu'elle est de classe `Date`: +```{r} +class(data$date) +``` + +Les points sont dans l'ordre chronologique inverse, il est donc utile de les trier: +```{r} +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: +```{r} +all(diff(data$date) == 7) +``` -Maintenant, à vous de jouer! Vous pouvez effacer toutes ces informations et les remplacer par votre document computationnel. +### Inspection + +Regardons enfin à quoi ressemblent nos données ! +```{r} +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. Le creux des incidences se trouve en fin d'été/début d'automne. +```{r} +with(tail(data, 200), plot(date, inc, type="l", xlab="Date", ylab="Incidence hebdomadaire")) +``` + +## L'incidence annuelle + +### Calcul + +On utilise la recommendation faite par les instructeurs du MOOC RR et on choisit le 1er septembre comme début de chaque période annuelle. Du reste, le pic de l'épidémie se situe à cheval entre deux années civiles. Nous définissons donc la période de référence entre deux minima de l'incidence, du 1er septembre de l'année $N$ au 1er septembre 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 la varicelle est très faible en automne, 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, mais il s'agit ici d'une précision inutile puisqu'on a vu avant qu'il n'yavait aucune données manquantes (on le laisse ici par habitude) +```{r} +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. De même, pour une exécution en mai 2021, les données se terminent le 29 mars 2021, c'est-à-dire avant le 1er septembre 2021, ce qui ne nous permet pas d'inclure non plus l'année 2021. +```{r} +annees = 1991:2020 +``` + +Nous créons un nouveau jeu de données pour l'incidence annuelle, en applicant la fonction `pic_annuel` à chaque année: +```{r} +inc_annuelle = data.frame(annee = annees, + incidence = sapply(annees, pic_annuel)) +head(inc_annuelle) +tail(inc_annuelle) +``` +On observe une importante diminution de l'incidence en 2020 par rapport aux années précédentes. On pourrait avancer quelques pistes d'explications, mais ce n'est pas le sujet ici :) + +### Inspection + +Voici les incidences annuelles en graphique: +```{r} +plot(inc_annuelle, type="p", xlab="Année", ylab="Incidence annuelle") +``` + +### Identification des épidémies les plus fortes + +Une liste triée par ordre décroissant d'incidence annuelle permet de plus facilement repérer les valeurs les plus élevées : +```{r} +head(inc_annuelle[order(-inc_annuelle$incidence),]) +``` +C'est donc l'année 2009 qui a vu survenir l'épidémie de varicelle la plus forte en France. + + +Une liste triée cette fois-ci par ordre croissant d'incidence annuelle permet de plus facilement repérer les valeurs les plus faibles : 2020 puis 2002 +```{r} +head(inc_annuelle[order(inc_annuelle$incidence),]) +``` +C'est donc l'année 2002 qui a vu survenir l'épidémie de varicelle la plus faible en France. -- 2.18.1