mooc-rr/module3/exo2/exercice_R_fr.org

Incidence de la varicelle

• Préparation des données
  • Récupération sur le site Sentinelles
  • Téléchargement ou lecture locale
  • Recherche de données manquantes
  • Extraction des colonnes utilisées
  • Vérification
  • Conversions
  • Vérification des dates
  • Passage Python -> R
  • Inspection
• Étude de l'incidence annuelle
  • Calcul de l'incidence annuelle
  • Inspection
  • Identification des épidémies les plus fortes

Préparation des données

Récupération sur le site Sentinelles

Les données d'incidence de la varicelle sont disponibles sur le site du Réseau Sentinelles.

https://www.sentiweb.fr/datasets/incidence-PAY-7.csv

Voici l'explication des colonnes donnée 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)

L'indication d'une semaine calendaire en format ISO-8601 est populaire en Europe, mais peu utilisée aux Etats-Unis. Ceci explique peut-être que peu de logiciels savent gérer ce format. Le langage Python le fait depuis la version 3.6. Nous utilisons donc ce langage pour la préparation de nos données, ce qui a l'avantage de ne nécessiter aucune bibliothèque supplémentaire. (Note: nous expliquerons dans le module 4 pourquoi il est avantageux pour la réproductibilité de se limiter à un minimum de bibliothèques.)

Téléchargement ou lecture locale

Si le fichier des données n'existe pas localement nous le téléchargeons sinon on utilise le fichier local existant. Nous commençons par l'extraction des données qui nous intéressent. D'abord nous découpons le contenu du fichier en lignes, dont nous jetons la première qui ne contient qu'un commentaire. Les autres lignes sont découpées en colonnes.

import os
import urllib.request

data_file = "incidence-varicelle.csv"

if not os.path.exists(data_file):
    print("Downloading from source")
    urllib.request.urlretrieve(data_url, data_file)

data = open(data_file, 'rb').read()
lines = data.decode('latin-1').strip().split('\n')
data_lines = lines[1:]
table = [line.split(',') for line in data_lines]

Regardons ce que nous avons obtenu:

table[:5]

week	indicator	inc	inc_low	inc_up	inc100	inc100_low	inc100_up	geo_insee	geo_name
202245	7	3942	1794	6090	6	3	9	FR	France
202244	7	4271	2231	6311	6	3	9	FR	France
202243	7	5863	3302	8424	9	5	13	FR	France
202242	7	3770	1950	5590	6	3	9	FR	France

Recherche de données manquantes

On vérifie s'il y a ou non des données manquantes dans le jeu de données.

valid_table = []
for row in table:
    missing = any([column == '' for column in row])
    if missing:
        print(row)
    else:
        valid_table.append(row)

valid_table[:5]

week	indicator	inc	inc_low	inc_up	inc100	inc100_low	inc100_up	geo_insee	geo_name
202245	7	3942	1794	6090	6	3	9	FR	France
202244	7	4271	2231	6311	6	3	9	FR	France
202243	7	5863	3302	8424	9	5	13	FR	France
202242	7	3770	1950	5590	6	3	9	FR	France

Extraction des colonnes utilisées

Les deux colonnes qui nous intéressent sont la 1ère ("week") et la 3ème ("inc") l'incidence.

week = [row[0] for row in valid_table]
assert week[0] == 'week'
del week[0]
inc = [row[2] for row in valid_table]
assert inc[0] == 'inc
del inc[0]
data = list(zip(week, inc))

Regardons les premières et les dernières lignes. Nous insérons None pour indiquer à org-mode la séparation entre les trois sections du tableau: en-tête, début des données, fin des données.

[('week', 'inc'), None] + data[:5] + [None] + data[-5:]

week	inc
202245	3942
202244	4271
202243	5863
202242	3770
202241	4177
199101	15565
199052	19375
199051	19080
199050	11079
199049	1143

Vérification

Il est toujours prudent de vérifier si les données semblent crédibles. Nous savons que les semaines sont données par six chiffres (quatre pour l'année et deux pour la semaine), et que les incidences sont des nombres entiers positifs.

for week, inc in data:
    if len(week) != 6 or not week.isdigit():
        print("Valeur suspecte dans la colonne 'week': ", (week, inc))
    if not inc.isdigit():
        print("Valeur suspecte dans la colonne 'inc': ", (week, inc))

Pas de valeurs suspectes détectées !.

Conversions

Pour faciliter les traitements suivants, nous remplaçons les numéros de semaine ISO par les dates qui correspondent aux lundis. A cette occasion, nous trions aussi les données par la date, et nous transformons les incidences en nombres entiers.

import datetime
converted_data = [(datetime.datetime.strptime(year_and_week + ":1" , '%G%V:%u').date(),
                  int(inc))
                  for year_and_week, inc in data]
converted_data.sort(key = lambda record: record[0])

Regardons de nouveau les premières et les dernières lignes:

str_data = [(str(date), str(inc)) for date, inc in converted_data]
[('date', 'inc'), None] + str_data[:5] + [None] + str_data[-5:]

date	inc
1990-12-03	1143
1990-12-10	11079
1990-12-17	19080
1990-12-24	19375
1990-12-31	15565
2022-10-10	4177
2022-10-17	3770
2022-10-24	5863
2022-10-31	4271
2022-11-07	3942

Vérification des dates

Nous faisons encore une vérification: nos dates doivent être séparées d'exactement une semaine, sauf autour du point manquant.

dates = [date for date, _ in converted_data]
for date1, date2 in zip(dates[:-1], dates[1:]):
    if date2-date1 != datetime.timedelta(weeks=1):
        print(f"Il y a {date2-date1} entre {date1} et {date2}")

Il n'y a pas d'écart inférieur ou supérieur à 7 jours entre nos données.

Passage Python -> R

Nous passons au langage R pour inspecter nos données, parce que l'analyse et la préparation de graphiques sont plus concises en R, sans nécessiter aucune bibliothèque supplémentaire.

Nous utilisons le mécanisme d'échange de données proposé par org-mode, ce qui nécessite un peu de code Python pour transformer les données dans le bon format.

[('date', 'inc'), None] + [(str(date), inc) for date, inc in converted_data]

data$date <- as.Date(data$date)
summary(data)

      date                 inc       
 Min.   :1990-12-03   Min.   :  161  
 1st Qu.:1998-11-26   1st Qu.: 7135  
 Median :2006-11-20   Median :12359  
 Mean   :2006-11-20   Mean   :12420  
 3rd Qu.:2014-11-13   3rd Qu.:16984  
 Max.   :2022-11-07   Max.   :36298

Inspection

Regardons enfin à quoi ressemblent nos données !

plot(data, type="l", xlab="Date", ylab="Incidence hebdomadaire")

Un zoom sur les dernières années montre mieux la situation des pics en hiver. Le creux des incidences se trouve en été.

plot(tail(data, 100), type="l", xlab="Date", ylab="Incidence hebdomadaire")

Étude de l'incidence annuelle

Calcul de l'incidence annuelle

É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 du syndrome grippal est très faible en été, cette modification ne risque pas de fausser nos conclusions.

Voici une fonction qui calcule l'incidence annuelle en appliquant ces conventions.

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 octobre 1984, ce qui ne permet pas de quantifier complètement le pic attribué à l'année 1985. Nous le supprimons donc de notre analyse. Pour la même raison, nous arrêtons en 2018. Nous devons attendre les données pour juillet 2019 avant d'augmenter la dernière année à 2019.

annees <- 1991:2021

inc_annuelle = data.frame(annee = annees,
                          incidence = sapply(annees, pic_annuel))
head(inc_annuelle)

1991	553895
1992	834935
1993	642921
1994	662750
1995	651333
1996	564994

tail(inc_annuelle)

2016	779816
2017	552906
2018	539765
2019	584116
2020	221183
2021	377933

Inspection

Voici les incidences annuelles en graphique.

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:

head(inc_annuelle[order(-inc_annuelle$incidence),])

   annee incidence
19  2009    841233
2   1992    834935
20  2010    834077
26  2016    779816
14  2004    778914
13  2003    760765

L'année où l'incidence a été la plus forte :

inc_annuelle[which.max(inc_annuelle$incidence),]

   annee incidence
19  2009    841233

L'année où l'incidence a été la plus faible :

inc_annuelle[which.min(inc_annuelle$incidence),]

   annee incidence
30  2020    221183

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.

hist(inc_annuelle$incidence, breaks=10, xlab="Incidence annuelle", ylab="Nb d'observations", main="")