diff --git a/module3/exo3/exercice_en.Rmd b/module3/exo3/exercice_en.Rmd index 13b258ddd0da29bc3bf08c64b6a1db742f6d5409..00128ba918441fce3154ed30e538e2376d5119b4 100644 --- a/module3/exo3/exercice_en.Rmd +++ b/module3/exo3/exercice_en.Rmd @@ -1,33 +1,175 @@ --- -title: "Your title" -author: "Your name" -date: "Today's date" -output: html_document +title: "CO2 concentration in the atmosphere since 1958" +author: "Konrad Hinsen" +output: + pdf_document: + toc: true + html_document: + toc: true + theme: journal +documentclass: article +classoption: a4paper +header-includes: +- \usepackage[french]{babel} +- \usepackage[upright]{fourier} +- \hypersetup{colorlinks=true,pagebackref=true} --- - ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` -## Some explanations +## Data preprocessing + +in 1958, Charles David Keeling started measuring the CO2 concentration in the atmosphere at Mauna Loa observatory, Hawaii, USA.[Réseau Sentinelles](http://www.sentiweb.fr/). 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: +```{r} +data_url = "https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/monthly/monthly_in_situ_co2_mlo.csv" +``` +### Download + +```{r} +data_file = "monthly_in_situ_co2_mlo.csv" +if (!file.exists(data_file)) { + download.file(data_url, data_file, method="auto") +} +``` + +This is the documentation of the data from [the download site]: + +| Column name | Description | +|--------------+-------------------------------------------------------------------------------------------------------------------------------| +| Columns 1-4 give the dates in several redundant formats. | +| Column 5 below gives monthly Mauna Loa CO2 concentrations in micro-mol CO2 per mole (ppm), reported on the 2012 SIO manometric mole | +| fraction scale. This is the standard version of the data most often sought. The monthly values have been adjusted to 24:00 hours on | +| the 15th of each month. | +| Column 6 gives the same data after a seasonal adjustment to remove the quasi-regular seasonal cycle. The adjustment involves subtracting | +| from the data a 4-harmonic fit with a linear gain factor. | +| Column 7 is a smoothed version of the data generated from a stiff cubic spline function plus 4-harmonic functions with linear gain. | +| Column 8 is the same smoothed version with the seasonal cycle removed. | +| Column 9 is identical to Column 5 except that the missing values from Column 5 have been filled with values from Column 7. | +| Column 10 is identical to Column 6 except missing values have been filled with values from Column 8. Missing values are denoted by -99.99 | +|--------------+-------------------------------------------------------------------------------------------------------------------------------| + +### Lecture + + +The first line of the CSV file is a comment, which we ignore with `skip=1`. +```{r} +data = read.csv(data_file, 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,] +``` + +The two relevant columns for us are `week` and `inc`. Let's verify their classes: +```{r} +class(data$week) +class(data$inc) +``` +Integers, fine! + +### 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} +library(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: -This is an R Markdown document that you can easily export to HTML, PDF, and MS Word formats. For more information on R Markdown, see . +```{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")) +``` -When you click on the button **Knit**, the document will be compiled in order to re-execute the R code and to include the results into the final document. As we have shown in the video, R code is inserted as follows: +## Annual incidence -```{r cars} -summary(cars) +### 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 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. + +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) { + debut = paste0(year-1,"-08-01") + fin = paste0(year,"-08-01") + semaines = data$date > debut & data$date <= fin + sum(data$inc[semaines], na.rm=TRUE) + } ``` -It is also straightforward to include figures. For example: +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. +```{r} +years = 1986:2018 +``` -```{r pressure, echo=FALSE} -plot(pressure) +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) ``` -Note the parameter `echo = FALSE` that indicates that the code will not appear in the final version of the document. We recommend not to use this parameter in the context of this MOOC, because we want your data analyses to be perfectly transparent and reproducible. +### Inspection -Since the results are not stored in Rmd files, you should generate an HTML or PDF version of your exercises and commit them. Otherwise reading and checking your analysis will be difficult for anyone else but you. +A plot of the annual incidences: +```{r} +plot(annnual_inc, type="p", xlab="Année", ylab="Annual incidence") +``` + +### Identification of the strongest epidemics -Now it's your turn! You can delete all this information and replace it by your computational document. +A list sorted by decreasing annual incidence makes it easy to find the most 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="") +```