@@ -20,25 +20,31 @@ For running this analysis, you need the following software:
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
#+BEGIN_SRC python :results output :exports both
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 emacs-lisp :results output
#+RESULTS:
#+BEGIN_SRC emacs-lisp :results output :exports both
(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
#+BEGIN_SRC emacs-lisp :results output :exports both
(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:
...
...
@@ -63,27 +69,54 @@ This is the documentation of the data from [[https://ns.sentiweb.fr/incidence/cs
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.
In order to protect us in case the Réseau Sentinelles Web server disappears or is
modified, we make a local copy of this dataset that we store together
with our analysis. It is unnecessary and even risky to download the
data at each execution, because in case of a malfunction we might be
replacing our file by a corrupted version. Therefore we download the
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.
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
#+BEGIN_SRC python :results silent :exports both
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
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
#+BEGIN_SRC python :results output :exports both
for week, inc in data:
if len(week) != 6 or not week.isdigit():
print("Suspicious value in column 'week': ", (week, inc))
...
...
@@ -120,12 +171,14 @@ for week, inc in data:
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.
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.
This R function computes the annual incidence as defined above:
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 only when all data up to July 2019 is available.
#+BEGIN_SRC R :results silent
years <- 1986:2018
#+BEGIN_SRC R :results silent :exports both
years <- 1986:2023
#+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
#+BEGIN_SRC R :results value :exports both
annnual_inc = data.frame(year = years,
incidence = sapply(years, yearly_peak))
head(annnual_inc)
#+END_SRC
#+RESULTS:
| 1986 | 5100540 |
| 1987 | 2861556 |
| 1988 | 2766142 |
| 1989 | 5460155 |
| 1990 | 5233987 |
| 1991 | 1660832 |
** Inspection
A plot of the annual incidence:
#+BEGIN_SRC R :results output graphics :file annual-inc-plot.png
#+BEGIN_SRC R :results graphics file :file annual-inc-plot.png :exports both
A list sorted by decreasing annual incidence makes it easy to find the most important ones:
#+BEGIN_SRC R :results output
#+BEGIN_SRC R :results output :exports both
head(annnual_inc[order(-annnual_inc$incidence),])
#+END_SRC
#+RESULTS:
: year incidence
: 4 1989 5460155
: 5 1990 5233987
: 1 1986 5100540
: 28 2013 4182265
: 25 2010 4085126
: 14 1999 3897443
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 graphics :file annual-inc-hist.png
#+BEGIN_SRC R :results output graphics file :file annual-inc-hist.png :exports both
hist(annnual_inc$incidence, breaks=10, xlab="Annual incidence", ylab="Number of observations", main="")
