---
title: "Analyse de l'incidence de la varicelle"
author: "Camille Launay"
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)
```

## Préparation des données

Les données de l'incidence du syndrome grippal 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"
```

Je télécharge une copie locale du fichier pour éviter de le télécharger à chaque fois qu'on lance le script et pour éviter les problèmes si la version en ligne est modifiée.
Je télécharge le fichier depuis Internet si la copie locale n'existe pas.
```{r}
data_file = "incidence-PAY-7.csv"
if (!file.exists(data_file)) {
    download.file(data_url, data_file, method="auto")
}
```

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)                                                       |

La première ligne du fichier CSV est un commentaire, que nous ignorons en précisant `skip=1`.
### Lecture
```{r}
data = read.csv(data_file, skip=1)
```

Regardons ce que nous avons obtenu:
```{r}
head(data)
tail(data)
```

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,]
```
_Pas de point manquant_

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

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}
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:

```{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))
}
```

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))
```

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)
```

### 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 en hiver. Le creux des incidences se trouve en été.
```{r}
with(tail(data, 200), plot(date, inc, type="l", xlab="Date", ylab="Incidence hebdomadaire"))
```

## L'incidence annuelle

### Calcul

É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 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 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.
```{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é à 1990. Nous l'enlevons donc de notre analyse. Les données se terminent en __mars 2020__, ce qui ne nous permet pas d'inclure l'année 2019.
```{r}
annees = 1991:2018
```

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)
```

### 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),])
```

Une liste triée par ordre croissant d'incidence annuelle permet de plus facilement repérer les valeurs les plus faibles:
```{r}
head(inc_annuelle[order(inc_annuelle$incidence),])
```

Enfin, un histogramme montre bien que les épidémies fortes, qui touchent environ 10% de la population française, sont assez rares: il y en eu trois au cours des 35 dernières années.
```{r}
hist(inc_annuelle$incidence, breaks=10, xlab="Incidence annuelle", ylab="Nb d'observations", main="")
```