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:
| ~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.
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
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.
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.
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 ~<p~, ~<P~ or ~<PP~
followed by ~Tab~.
This R function computes the annual incidence as defined above:
#+BEGIN_SRC R :results silent
yearly_peak = function(year) {
debut = paste0(year-1,"-09-01")
fin = paste0(year,"-09-01")
semaines = data$date > 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
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="")
