Commit 2582c3da authored by Julien Quénon's avatar Julien Quénon

Exercice incidence varicelle

parent 51080c84
---
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 <http://rmarkdown.rstudio.com>.
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.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment