#+TITLE: Analyse de l'incidence du syndrôme grippal # Définit automatiquement les en-têtes des sessions (une même session # partagée et du code et résultats exportés #+PROPERTY: header-args :session :exports both # On définit une variable data-url en org-mode #+NAME: data-url https://www.sentiweb.fr/datasets/incidence-PAY-3.csv Nous allons charger un fichier CSV de données via la librairie Pandas (en python). On récupère l'URL du CSV via une variable stockée en org-mode: data-url(voir dans l'en-tête l'opérande :var) #+begin_src python :results silent :var data_url=data-url import pandas as pd table = pd.read_csv(data_url, encoding='latin', header=1, index_col=0) #+end_src #+RESULTS: #+begin_example week indicator inc ... inc100_up geo_insee geo_name 0 202015 3 0 ... 0.0 FR France 1 202014 3 0 ... 0.0 FR France 2 202013 3 0 ... 0.0 FR France 3 202012 3 8321 ... 17.0 FR France 4 202011 3 101704 ... 166.0 FR France ... ... ... ... ... ... ... ... 1845 198448 3 78620 ... 176.0 FR France 1846 198447 3 72029 ... 163.0 FR France 1847 198446 3 87330 ... 195.0 FR France 1848 198445 3 135223 ... 308.0 FR France 1849 198444 3 68422 ... 213.0 FR France [1850 rows x 10 columns] #+end_example La lecture semble s'être bien passée. Regardons les premières lignes du jeu: #+begin_src python :results value table[:5] #+end_src #+RESULTS: : indicator inc inc_low inc_up inc100 inc100_low inc100_up geo_insee geo_name : week : 202015 3 0 0.0 0.0 0 0.0 0.0 FR France : 202014 3 0 0.0 0.0 0 0.0 0.0 FR France : 202013 3 0 0.0 0.0 0 0.0 0.0 FR France : 202012 3 8321 5873.0 10769.0 13 9.0 17.0 FR France : 202011 3 101704 93652.0 109756.0 154 142.0 166.0 FR France Concentrons-nous sur les colonnes décrivant les semaines(week) et l'estimation du nombre de cas grippaux (inc) #+begin_src python :results output :exports both table[['inc']].head() #+end_src #+RESULTS: : inc : week : 202015 0 : 202014 0 : 202013 0 : 202012 8321 : 202011 101704 Maintenant, vérifions que les dates sont écrite sous la forme d'un nombre entier comportant 6 chiffres. Vérifier exactement le format ISO-8601 requiert de complexifier davantage l'expression régulière(ce qui n'est pas l'objectif ici) #+begin_src python :results value :exports both week_correctly_written = table.index.astype(str).str.match(r'^\d{6}$') table[week_correctly_written == False] #+end_src #+RESULTS: : Empty DataFrame : Columns: [indicator, inc, inc_low, inc_up, inc100, inc100_low, inc100_up, geo_insee, geo_name] : Index: [] L'Index retourné est vide. On conclut que les semaines sont toutes bien formatées dans ce fichier. On convertit les index initialement au format ISO vers le format datetime #+begin_src python :results output :exports both table.index = pd.to_datetime(table.index.astype(str) + ':1', format='%G%V:%u') #+end_src #+RESULTS: Observons la table réindexée au format datetime et triée selon les dates croissantes #+begin_src python :results value :exports both table.sort_index(inplace=True) table[['inc']] #+end_src #+RESULTS: #+begin_example inc week 1984-10-29 68422 1984-11-05 135223 1984-11-12 87330 1984-11-19 72029 1984-11-26 78620 ... ... 2020-03-09 101704 2020-03-16 8321 2020-03-23 0 2020-03-30 0 2020-04-06 0 [1850 rows x 1 columns] #+end_example On regarde si les enregistrements successifs sont bien espacés de 7 jours: #+begin_src python :results value :exports both is_weekly_spaced = table.index.to_series().diff() == pd.Timedelta('7 days') # On affiche les enregistrements mal espacés s'il y a (le premier est évité car est l'origine) table.loc[is_weekly_spaced == False][1:] #+end_src #+RESULTS: : Empty DataFrame : Columns: [indicator, inc, inc_low, inc_up, inc100, inc100_low, inc100_up, geo_insee, geo_name] : Index: [] On voit que les semaines sont bien espacées :) Formattons un peu les données afin de pouvoir les transmettre à du code R qui nous facilitera l'analyse descriptive. #+NAME: data-for-R #+begin_src python :results silent :exports both [('date', 'inc'), None] + table['inc'].reset_index().astype({'week': str}).values.tolist() #+end_src #+begin_src R :results output :var data=data-for-R data$date <- as.Date(data$date) summary(data) #+end_src #+RESULTS: : : date inc : Min. :1984-10-29 Min. : 0 : 1st Qu.:1993-09-07 1st Qu.: 4990 : Median :2002-07-18 Median : 15987 : Mean :2002-07-18 Mean : 61908 : 3rd Qu.:2011-05-28 3rd Qu.: 50572 : Max. :2020-04-06 Max. :1001824 Observons maintenant nos données #+BEGIN_SRC R :results output graphics :file inc-plot.png :exports both plot(data, type="l", xlab="Date", ylab="Incidence hebdomadaire") #+END_SRC #+RESULTS: [[file:inc-plot.png]] Regardons de plus près la saisonnalité qu'on devine: #+BEGIN_SRC R :results output graphics :file inc-plot-zoom.png :exports both plot(tail(data, 200), type="l", xlab="Date", ylab="Incidence hebdomadaire") #+END_SRC #+RESULTS: [[file:inc-plot-zoom.png]] On observe qu'il y a bien des pics de contamination qui reviennent chaque début d'année, en hiver. Ecrivons une fonction qui permet de calculer le cumul des contaminations estimées sur une année complète démarrant pendant un creux stable dans le temps: le premier aout. Nous souhaitons comparer les contaminations d'une année à l'autre, il est donc important de prendre comme référence un point d'origine stable au cours des années. #+BEGIN_SRC R :results silent :exports both pic_annuel = function(annee) { debut = paste0(annee-1,"-08-01") fin = paste0(annee,"-08-01") semaines = data$date > debut & data$date <= fin sum(data$inc[semaines], na.rm=TRUE) } #+END_SRC Afin de ne pas sous-estimer le nombre de contaminations, on commence à calculer à partir de la première année pleine: 1985. 2020 risque d'être légèrement sous-estimée bien que la plupart des contaminations grippales aient lieu avant avril. Cependant, la pandémie de covid-19 risque de contre-balancer les chiffres. Il faudrait être davantage vigilant sur la saisie des données. Il peut avoir eu confusion entre la grippe et ce coronavirus au moment de la détection. #+BEGIN_SRC R :results silent :exports both annees <- 1985:2020 #+END_SRC Construisons maintenant la dataframe avec les contaminations annuelles #+BEGIN_SRC R :results value :exports both inc_annuelle = data.frame(annee = annees, incidence = sapply(annees, pic_annuel)) head(inc_annuelle) #+END_SRC #+RESULTS: | 1985 | 5779196 | | 1986 | 5100540 | | 1987 | 2861556 | | 1988 | 2766142 | | 1989 | 5460155 | | 1990 | 5233987 | Affichons ces données #+BEGIN_SRC R :results output graphics :file annual-inc-plot.png :exports both plot(inc_annuelle, type="p", xlab="Année", ylab="Incidence annuelle") #+END_SRC #+RESULTS: [[file:annual-inc-plot.png]] Identifions maintenant les années les plus touchées #+BEGIN_SRC R :results output :exports both head(inc_annuelle[order(-inc_annuelle$incidence),]) #+END_SRC #+RESULTS: : annee incidence : 1 1985 5779196 : 5 1989 5460155 : 6 1990 5233987 : 2 1986 5100540 : 29 2013 4182265 : 26 2010 4085126 Il est intéressant de voir que les années 1985, 1989, 1990 et 1986 ont été bien plus marquée que les suivantes. Voici un histogramme illustrant que les épidémies fortes(touchant plus de 7-8% de la population) sont assez rares. #+BEGIN_SRC R :results output graphics :file annual-inc-hist.png :exports both hist(inc_annuelle$incidence, breaks=5, xlab="Incidence annuelle", ylab="Nb d'observations", main="") #+END_SRC #+RESULTS: [[file:annual-inc-hist.png]]