---
title: "Analyse de l'incidence de la varicelle"
author: "Xavier Nardou"
date: "25 Mars 2020"
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)
```

On commence par récupérer les données sur la varicelle de 1991 à aujourd'hui.

```{r}
data_url = "https://www.sentiweb.fr/datasets/incidence-PAY-7.csv"
```

On va paramétrer pour faire en sorte que la copie locale des données soit utilisée de préférence, et qu'elle ne soit pas téléchargée à chaque fois.
```{r}

```

```{r}
data_file="incidence-PAY-7.csv"
if(!file.exists(data_file)){
  download.file(data_url, data_file, method ="auto")
  
}
```

On peut maintenant charger le fichier à partir de la copie locale qu'on a faite.

```{r}
data = read.csv(data_file, skip=1, sep=",")
```

Regardons la structure de nos données

```{r}
head(data)
tail(data)
str(data)
```
Ce sont les bonnes classes correspondantes, incidence et week sont des entiers.
Maintenant, y'a t-il des données manquantes ? 


```{r}
na_records = apply(data, 1, function (x) any(is.na(x)))
data[na_records,]
```
Il n'y a pas de données manquantes dans ce jeu de données. 

### Conversion des numéros de semaine

Utilisation de la library parsedate
```{r}
library(parsedate)
```

On remplace les semaines par les dates qui correspondent aux lundis

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

Et on applique à tous les points

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

La distribution annuelle a l'air plus étalée que la grippe. On va faire unu zoom sur les dernières années.
```{r}
with(tail(data, 200), plot(date, inc, type="l", xlab="Date", ylab="Incidence hebdomadaire"))
```
Effectivement le pic est moins abrupt, et on voit qu'il y a un creux vers le mois de Septembre, donc on va séparer les années comme ça. 

## L'incidence annuelle

### Calcul

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. 

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

On va retirer la première année (1991) et 2020 car les données vont être incomplètes pour ces années là.

```{r}
annees = 1992:2019
```

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

Pour voir les épidémies les plus faibles

```{r}
head(inc_annuelle[order(inc_annuelle$incidence),])
```

Et un petit histogramme
```{r}
hist(inc_annuelle$incidence, breaks=10, xlab="Incidence annuelle", ylab="Nb d'observations", main="")
```