From 4d9aeefee133eeee20654105e9b570f9be7499ef Mon Sep 17 00:00:00 2001 From: feee8ecec645c53ceeda57e948ea51be Date: Wed, 8 Sep 2021 12:09:05 +0000 Subject: [PATCH] Upload a draft submission on subject 2. --- module3/exo3/subject2_exercice_en.Rmd | 191 ++++++++++++++++++++++++++ 1 file changed, 191 insertions(+) create mode 100644 module3/exo3/subject2_exercice_en.Rmd diff --git a/module3/exo3/subject2_exercice_en.Rmd b/module3/exo3/subject2_exercice_en.Rmd new file mode 100644 index 0000000..c6078da --- /dev/null +++ b/module3/exo3/subject2_exercice_en.Rmd @@ -0,0 +1,191 @@ +--- +title: "Purchasing power of English workers from the 16th to the 19th century" +author: "Eleni Gkiouzepi" +output: + html_document: + toc: true + theme: journal + pdf_document: + toc: true +documentclass: article +classoption: a4paper +header-includes: +- \usepackage[english]{babel} +- \usepackage[upright]{fourier} +- \hypersetup{colorlinks=true,pagebackref=true} +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +## Introduction + +[William Playfair]("https://en.wikipedia.org/wiki/William_Playfair") was one of the pioneers of the graphical presentation of data, being credited in particular with the invention of the histogram. One of his famous graphs, taken from his book "[A Letter on Our Agricultural Distresses, Their Causes and Remedies]("https://books.google.fr/books/about/A_Letter_on_Our_Agricultural_Distresses.html?id=aQZGAQAAMAAJ")", shows [the evolution of the wheat price and average salaries from 1565 to 1821]("https://fr.wikipedia.org/wiki/William_Playfair#/media/File:Chart_Showing_at_One_View_the_Price_of_the_Quarter_of_Wheat,_and_Wages_of_Labour_by_the_Week,_from_1565_to_1821.png"). +
+ +## Data preprocessing + +Playfair did not publish his raw data because reproducibility was not yet a topic in his days. The numbers obtained from a scan of the graph are available for download, the CSV version being most convenient. + +```{r} +data_url = "https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/HistData/Wheat.csv" +``` + +This is the documentation of the data from [the Vincent Arel-Bundock github site](https://vincentarelbundock.github.io/Rdatasets/doc/HistData/Wheat.html): + +| Column name | Description | +|--------------+---------------------------------------------------------------------------------------------------------------------------| +| `` | observation number | +| `Year` | Year, in intervals of 5 from 1565 to 1821: a numeric vector | +| `Wheat` | Price of Wheat (Shillings/Quarter bushel equal 15 British pounds or about 6,8 kg): a numeric vector | +| `Wages` | Weekly wage (Shillings): a numeric vector | + +* _Important notice_: Until 1971, the pound sterling was divided into 20 shillings, and a shilling into 12 pence. + +### If the local file does not exist, download the data and put them into the local file + +```{r} +destfile = "Wheat.csv" + +if(!file.exists(destfile)){ + res <- tryCatch(download.file(data_url, + destfile, + method="auto"), + error=function(e) 1) +} +``` + +### Read the local CSV file. + +The first line of the CSV file is a comment, which we ignore with `skip=1`. + +```{r} +data = read.csv(destfile, skip=1) +``` + +Let's have a look at what we got: +```{r} +head(data) +tail(data) +``` + +Are there missing data points? +```{r} +na_records = apply(data, 1, function (x) any(is.na(x))) +data[na_records,] +``` + +Well, there three values missing. The weekly wages from 1815 to 1821. + +The three relevant columns for us are `Year`, `Wheat` and `Wages`. Let's verify their classes: +```{r} +class(data$Year) +class(data$Wheat) +class(data$Wages) +``` + + +### Conversion of the week numbers + +Date handling is always a delicate subject. There are many conventions that are easily confused. Our dataset uses the [ISO-8601](https://en.wikipedia.org/wiki/ISO_8601) week number format, which is popular in Europe but less so in North America. In `R`, it is handled by the library [parsedate](https://cran.r-project.org/package=parsedate): + +```{r} +if(!require(parsedate)) install.packages("parsedate") + require(parsedate) +``` + +In order to facilitate the subsequent treatment, we replace the ISO week numbers by the dates of each week's Monday. This function does it for one value: + +```{r} +convert_week = function(w) { + ws = paste(w) + iso = paste0(substring(ws, 1, 4), "-W", substring(ws, 5, 6)) + as.character(parse_iso_8601(iso)) +} +``` + +We apply it to all points, creating a new column `date` in our data frame: +```{r} +data$date = as.Date(convert_week(data$week)) +``` + +Let's check that is has class `Date`: +```{r} +class(data$date) +``` + +The points are in inverse chronological order, so it's preferable to sort them: +```{r} +data = data[order(data$date),] +``` + +That's a good occasion for another check: our dates should be separated by exactly seven days: +```{r} +all(diff(data$date) == 7) +``` + +### Inspection + +Finally we can look at a plot of our data! +```{r} +plot(data$date, data$inc, type="l", xlab="Date", ylab="Weekly incidence") +``` + +A zoom on the last few years makes the peaks in winter stand out more clearly. +```{r} +with(tail(data, 200), plot(date, inc, type="l", xlab="Date", ylab="Weekly incidence")) +``` + +## Annual incidence + +### Computation + +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 September 1st of year $N$ to September 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. + +The argument `na.rm=True` in the sum indicates that missing data points are removed. This is a reasonable choice since there is only one missing point, whose impact cannot be very strong. +```{r} +yearly_peak = function(year) { + start = paste0(year-1,"-09-01") + end = paste0(year,"-09-01") + semaines = data$date > start & data$date <= end + sum(data$inc[semaines], na.rm=TRUE) + } +``` + +We must also be careful with the first and last years of the dataset. The data start in December 1990, meaning that we don't have all the data for the peak attributed to the year 1990. We therefore exclude it from the analysis. For the same reason, we define 2020 as the final year. We can increase this value to 2021 only when all data up to August 2021 is available. +```{r} +years = 1991:2020 +``` + +We make a new data frame for the annual incidence, applying the function `yearly_peak` to each year: +```{r} +annnual_inc = data.frame(year = years, + incidence = sapply(years, yearly_peak)) +head(annnual_inc) +``` + +### Inspection + +A plot of the annual incidences: +```{r} +plot(annnual_inc, type="p", xlab="Year", ylab="Annual incidence") +``` + +### Identification of the strongest and the weakest epidemics + +A list sorted by decreasing annual incidence makes it easy to find the most important ones: +```{r} +head(annnual_inc[order(-annnual_inc$incidence),]) +``` + +A list sorted by decreasing annual incidence makes it easy to find the least important ones: +```{r} +head(annnual_inc[order(annnual_inc$incidence),]) +``` + +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. +```{r} +hist(annnual_inc$incidence, breaks=10, xlab="Annual incidence", + ylab="Number of observations", main="") +``` -- 2.18.1