From 493188aa35fca17868bd876e7c9485612bd0d135 Mon Sep 17 00:00:00 2001 From: migvasc Date: Sat, 13 Mar 2021 19:07:57 +0100 Subject: [PATCH] solution to exercise 2 of module 3 --- module3/exo2/exercice_python_en.org | 396 +++++++++++++++++++++++----- 1 file changed, 331 insertions(+), 65 deletions(-) diff --git a/module3/exo2/exercice_python_en.org b/module3/exo2/exercice_python_en.org index 5782f49..6a24df9 100644 --- a/module3/exo2/exercice_python_en.org +++ b/module3/exo2/exercice_python_en.org @@ -1,6 +1,6 @@ -#+TITLE: Your title -#+AUTHOR: Your name -#+DATE: Today's date +#+TITLE: Analysis of chicken pox in france +#+AUTHOR: Miguel Felipe Silva Vasconcelos +#+DATE: 13/03/2021 #+LANGUAGE: en # #+PROPERTY: header-args :eval never-export @@ -11,84 +11,350 @@ #+HTML_HEAD: #+HTML_HEAD: -* Some explanations +#+PROPERTY: header-args :session :exports both -This is an org-mode document with code examples in R. Once opened in -Emacs, this document can easily be exported to HTML, PDF, and Office -formats. For more information on org-mode, see -https://orgmode.org/guide/. +* Foreword -When you type the shortcut =C-c C-e h o=, this document will be -exported as HTML. All the code in it will be re-executed, and the -results will be retrieved and included into the exported document. If -you do not want to re-execute all code each time, you can delete the # -and the space before ~#+PROPERTY:~ in the header of this document. +For running this analysis, you need the following software: -Like we showed in the video, Python code is included as follows (and -is exxecuted by typing ~C-c C-c~): +** Emacs 25 or higher +Older versions may suffice. For Emacs versions older than 26, org-mode must be updated to version 9.x. +** Python 3.6 or higher +We use the ISO 8601 date format, which has been added to Python's standard library with version 3.6. +#+BEGIN_SRC python :results output +import sys +if sys.version_info.major < 3 or sys.version_info.minor < 6: + print("Please use Python 3.6 (or higher)!") +#+END_SRC -#+begin_src python :results output :exports both -print("Hello world!") -#+end_src +#+RESULTS: + +#+BEGIN_SRC emacs-lisp :results output +(unless (featurep 'ob-python) + (print "Please activate python in org-babel (org-babel-do-languages)!")) +#+END_SRC + +#+RESULTS: + +** R 3.4 +We use only basic R functionality, so a earlier version might be OK, but we did not test this. + +#+BEGIN_SRC emacs-lisp :results output +(unless (featurep 'ob-R) + (print "Please activate R in org-babel (org-babel-do-languages)!")) +#+END_SRC + +#+RESULTS: + +* Data preprocessing + +The data on the incidence of influenza-like illness are available from the Web site of the [[http://www.sentiweb.fr/][Réseau Sentinelles]]. We download them as a file in CSV format, in which each line corresponds to a week in the observation period. Only the complete dataset, starting in 1984 and ending with a recent week, is available for download. The URL is: +#+NAME: data-url +https://www.sentiweb.fr/datasets/incidence-PAY-7.csv + +#+NAME: file-name +chickenpox.csv + +This is the documentation of the data from [[https://ns.sentiweb.fr/incidence/csv-schema-v1.json][the download site]]: + +| Column name | Description | +|--------------+---------------------------------------------------------------------------------------------------------------------------| +| ~week~ | ISO8601 Yearweek number as numeric (year*100 + week nubmer) | +| ~indicator~ | Unique identifier of the indicator, see metadata document https://www.sentiweb.fr/meta.json | +| ~inc~ | Estimated incidence value for the time step, in the geographic level | +| ~inc_low~ | Lower bound of the estimated incidence 95% Confidence Interval | +| ~inc_up~ | Upper bound of the estimated incidence 95% Confidence Interval | +| ~inc100~ | Estimated rate incidence per 100,000 inhabitants | +| ~inc100_low~ | Lower bound of the estimated incidence 95% Confidence Interval | +| ~inc100_up~ | Upper bound of the estimated rate incidence 95% Confidence Interval | +| ~geo_insee~ | Identifier of the geographic area, from INSEE https://www.insee.fr | +| ~geo_name~ | Geographic label of the area, corresponding to INSEE code. This label is not an id and is only provided for human reading | + +The [[https://en.wikipedia.org/wiki/ISO_8601][ISO-8601]] format is popular in Europe, but less so in North America. This may explain why few software packages handle this format. The Python language does it since version 3.6. We therefore use Python for the pre-processing phase, which has the advantage of not requiring any additional library. (Note: we will explain in module 4 why it is desirable for reproducibility to use as few external libraries as possible.) + +** Download +After downloading the raw data, we extract the part we are interested in. We first split the file into lines, of which we discard the first one that contains a comment. We then split the remaining lines into columns. + +#+BEGIN_SRC python :results output :var data_url=data-url :var file_name=file-name +from urllib.request import urlopen +import shutil +import os.path + + + +def downloadFile(): + result = urlopen(data_url) #makes the requisition for the file + out_file = open(file_name, 'wb') #tries save it to a file named influenzaincidence.csv + shutil.copyfileobj(result, out_file) #use shutil.copyfileobj if the file is large. See https://docs.python.org/dev/library/shutil.html#shutil.copyfileobj + result = result.read() + out_file.close() #close the file after downloading it + print("File downloaded!") + return result + + +def loadData(): + if os.path.isfile(file_name): + print("File Exists!") + else: + print("File not available locally... Trying to download it:") + downloadFile() + + file = open(file_name,"rb") #tries to open the file + data = file.read() + file.close() + return data + + +data = loadData() + +lines = data.decode('latin-1').strip().split('\n') +data_lines = lines[1:] +table = [line.split(',') for line in data_lines] +#+END_SRC + +#+RESULTS: +: File Exists! + +Let's have a look at what we have so far: +#+BEGIN_SRC python :results value +table[:2] +#+END_SRC + +#+RESULTS: +| week | indicator | inc | inc_low | inc_up | inc100 | inc100_low | inc100_up | geo_insee | geo_name | +| 202109 | 7 | 11766 | 8111 | 15421 | 18 | 12 | 24 | FR | France | + +** Checking for missing data +Unfortunately there are many ways to indicate the absence of a data value in a dataset. Here we check for a common one: empty fields. For completeness, we should also look for non-numerical data in numerical columns. We don't do this here, but checks in later processing steps would catch such anomalies. + +We make a new dataset without the lines that contain empty fields. We print those lines to preserve a trace of their contents. +#+BEGIN_SRC python :results output +valid_table = [] +for row in table: + missing = any([column == '' for column in row]) + if missing: + print(row) + else: + valid_table.append(row) +#+END_SRC + +#+RESULTS: + +** Extraction of the required columns +There are only two columns that we will need for our analysis: the first (~"week"~) and the third (~"inc"~). We check the names in the header to be sure we pick the right data. We make a new table containing just the two columns required, without the header. +#+BEGIN_SRC python :results silent +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)) +#+END_SRC + +Let's look at the first and last lines. We insert ~None~ to indicate to org-mode the separation between the three parts of the table: header, first lines, last lines. +#+BEGIN_SRC python :results value +[('week', 'inc'), None] + data[:5] + [None] + data[-5:] +#+END_SRC + +#+RESULTS: +| week | inc | +|--------+-------| +| 202109 | 11766 | +| 202108 | 11382 | +| 202107 | 13561 | +| 202106 | 13401 | +| 202105 | 12210 | +|--------+-------| +| 199101 | 15565 | +| 199052 | 19375 | +| 199051 | 19080 | +| 199050 | 11079 | +| 199049 | 1143 | + +** Verification +It is always prudent to verify if the data looks credible. A simple fact we can check for is that weeks are given as six-digit integers (four for the year, two for the week), and that the incidence values are positive integers. +#+BEGIN_SRC python :results output +for week, inc in data: + if len(week) != 6 or not week.isdigit(): + print("Suspicious value in column 'week': ", (week, inc)) + if not inc.isdigit(): + print("Suspicious value in column 'inc': ", (week, inc)) +#+END_SRC + +#+RESULTS: + +No problem - fine! + +** Date conversion +In order to facilitate the subsequent treatment, we replace the ISO week numbers by the dates of each week's Monday. This is also a good occasion to sort the lines by increasing data, and to convert the incidences from strings to integers. + +#+BEGIN_SRC python :results silent +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]) +#+END_SRC + +We'll look again at the first and last lines: +#+BEGIN_SRC python :results value +str_data = [(str(date), str(inc)) for date, inc in converted_data] +[('date', 'inc'), None] + str_data[:5] + [None] + str_data[-5:] +#+END_SRC + +#+RESULTS: +| date | inc | +|------------+-------| +| 1990-12-03 | 1143 | +| 1990-12-10 | 11079 | +| 1990-12-17 | 19080 | +| 1990-12-24 | 19375 | +| 1990-12-31 | 15565 | +|------------+-------| +| 2021-02-01 | 12210 | +| 2021-02-08 | 13401 | +| 2021-02-15 | 13561 | +| 2021-02-22 | 11382 | +| 2021-03-01 | 11766 | + +** Date verification +We do one more verification: our dates must be separated by exactly one week, except around the missing data point. +#+BEGIN_SRC python :results output +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"The difference between {date1} and {date2} is {date2-date1}") +#+END_SRC + +#+RESULTS: + +** Transfer Python -> R +We switch to R for data inspection and analysis, because the code is more concise in R and requires no additional libraries. + +Org-mode's data exchange mechanism requires some Python code for transforming the data to the right format. +#+NAME: data-for-R +#+BEGIN_SRC python :results silent +[('date', 'inc'), None] + [(str(date), inc) for date, inc in converted_data] +#+END_SRC + +In R, the dataset arrives as a data frame, which is fine. But the dates arrive as strings and must be converted. +#+BEGIN_SRC R :results output :var data=data-for-R +data$date <- as.Date(data$date) +summary(data) +#+END_SRC #+RESULTS: -: Hello world! +: +: date inc +: Min. :1990-12-03 Min. : 161 +: 1st Qu.:1998-06-25 1st Qu.: 7232 +: Median :2006-01-16 Median :12525 +: Mean :2006-01-16 Mean :12567 +: 3rd Qu.:2013-08-08 3rd Qu.:17134 +: Max. :2021-03-01 Max. :36298 -And now the same but in an Python session. With a session, Python's -state, i.e. the values of all the variables, remains persistent from -one code block to the next. The code is still executed using ~C-c -C-c~. -#+begin_src python :results output :session :exports both -import numpy -x=numpy.linspace(-15,15) -print(x) -#+end_src +** Inspection +Finally we can look at a plot of our data! +#+BEGIN_SRC R :results output graphics :file inc-plot.png +plot(data, type="l", xlab="Date", ylab="Weekly incidence") +#+END_SRC #+RESULTS: -#+begin_example -[-15. -14.3877551 -13.7755102 -13.16326531 -12.55102041 - -11.93877551 -11.32653061 -10.71428571 -10.10204082 -9.48979592 - -8.87755102 -8.26530612 -7.65306122 -7.04081633 -6.42857143 - -5.81632653 -5.20408163 -4.59183673 -3.97959184 -3.36734694 - -2.75510204 -2.14285714 -1.53061224 -0.91836735 -0.30612245 - 0.30612245 0.91836735 1.53061224 2.14285714 2.75510204 - 3.36734694 3.97959184 4.59183673 5.20408163 5.81632653 - 6.42857143 7.04081633 7.65306122 8.26530612 8.87755102 - 9.48979592 10.10204082 10.71428571 11.32653061 11.93877551 - 12.55102041 13.16326531 13.7755102 14.3877551 15. ] -#+end_example -Finally, an example for graphical output: -#+begin_src python :results output file :session :var matplot_lib_filename="./cosxsx.png" :exports results -import matplotlib.pyplot as plt -plt.figure(figsize=(10,5)) -plt.plot(x,numpy.cos(x)/x) -plt.tight_layout() -plt.savefig(matplot_lib_filename) -print(matplot_lib_filename) -#+end_src +A zoom on the last few years makes the peaks in winter stand out more clearly. +#+BEGIN_SRC R :results output graphics :file inc-plot-zoom.png +plot(tail(data, 200), type="l", xlab="Date", ylab="Weekly incidence") +#+END_SRC #+RESULTS: -[[file:./cosxsx.png]] -Note the parameter ~:exports results~, which indicates that the code -will not appear in the exported document. We recommend that in the -context of this MOOC, you always leave this parameter setting as -~:exports both~, because we want your analyses to be perfectly -transparent and reproducible. +* Study of the annual incidence -Watch out: the figure generated by the code block is /not/ stored in -the org document. It's a plain file, here named ~cosxsx.png~. You have -to commit it explicitly if you want your analysis to be legible and -understandable on GitLab. +** Computation of the annual incidence +Since the peaks of the epidemic happen in winter, near the transition between calendar years, we define the reference period for the annual incidence from August 1st of year /N/ to August 1st of year /N+1/. We label this period as year /N+1/ because the peak is always located in year /N+1/. The very low incidence in summer ensures that the arbitrariness of the choice of reference period has no impact on our conclusions. -Finally, don't forget that we provide in the resource section of this -MOOC a configuration with a few keyboard shortcuts that allow you to -quickly create code blocks in Python by typing ~ debut & data$date <= fin + sum(data$inc[semaines], na.rm=TRUE) + } +#+END_SRC -Now it's your turn! You can delete all this information and replace it -by your computational document. +We must also be careful with the first and last years of the +dataset. The data start in October 1984, meaning that we don't have +all the data for the peak attributed to the year 1985. We therefore +exclude it from the analysis. For the same reason, we define 2018 as +the final year. We can increase this value to 2019. (There is data for +2020, but the when the exercise was created, there was no such data. +#+BEGIN_SRC R :results silent +years <- 1991:2019 +#+END_SRC + +We make a new data frame for the annual incidence, applying the function ~yearly_peak~ to each year: +#+BEGIN_SRC R :results value +annnual_inc = data.frame(year = years, + incidence = sapply(years, yearly_peak)) +head(annnual_inc) +#+END_SRC + +#+RESULTS: +| 1991 | 553895 | +| 1992 | 834935 | +| 1993 | 642921 | +| 1994 | 662750 | +| 1995 | 651333 | +| 1996 | 564994 | + +** Inspection +A plot of the annual incidence: +#+BEGIN_SRC R :results output :file annual-inc-plot.png +plot(annnual_inc, type="p", xlab="Année", ylab="Annual incidence") +#+END_SRC + +#+RESULTS: + +** Identification of the strongest epidemics +A list sorted by decreasing annual incidence makes it easy to find the most important ones: +#+BEGIN_SRC R :results output +head(annnual_inc[order(-annnual_inc$incidence),]) +#+END_SRC + +#+RESULTS: +: year incidence +: 19 2009 841233 +: 2 1992 834935 +: 20 2010 834077 +: 26 2016 779816 +: 14 2004 778914 +: 13 2003 760765 + + +** Identification of the weakest epidemics +A list sorted by inccreasing annual incidence makes it easy to find the most important ones: +#+BEGIN_SRC R :results output +head(annnual_inc[order(annnual_inc$incidence),]) +#+END_SRC + +#+RESULTS: +: year incidence +: 12 2002 515343 +: 28 2018 539765 +: 27 2017 552906 +: 1 1991 553895 +: 6 1996 564994 +: 29 2019 584116 + + +Finally, a histogram clearly shows the few very strong epidemics, which affect about 10% of the French population, but are rare: there were three of them in the course of 35 years. The typical epidemic affects only half as many people. +#+BEGIN_SRC R :results output :file annual-inc-hist.png +hist(annnual_inc$incidence, breaks=10, xlab="Annual incidence", ylab="Number of observations", main="") +#+END_SRC + +#+RESULTS: -- 2.18.1