title: "Concentration de CO2 dans l'atmosphère depuis 1958"
author: "Thomas Pérot"
date: "12 juin 2023"
output: pdf_document
editor_options:
markdown:
wrap: sentence
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Quelques explications
## Introduction
Dans ce document nous allons présenter et analyser des données d'évolution de la concentration en CO2 atmosphérique provenant
du site de l'observatoire de Mauna Loa situé sur l'île d'Hawaï.
*C. D. Keeling, S. C. Piper, R. B. Bacastow, M. Wahlen, T. P. Whorf, M. Heimann, and
H. A. Meijer, Exchanges of atmospheric CO2 and 13CO2 with the terrestrial biosphere and
oceans from 1978 to 2000. I. Global aspects, SIO Reference Series, No. 01-06, Scripps
Institution of Oceanography, San Diego, 88 pages, 2001.*
## Chargement des données
Les données sont disponibles sur le site suivant : [Scripps CO2 Program](https://scrippsco2.ucsd.edu/data/atmospheric_co2/primary_mlo_co2_record.html).
Les premières 57 lignes du fichier donnent des informatins sur les données.
Les lignes 2 à 19 indiquent l'encroit des mesures, la source des données et d'autres informations sur les données.
Les lignes 23 à 37 indiquent comment citer les données.
Les lignes 41 à 56 décrivent le tableau de données.
Le tableau de données commence à la ligne 58 avec l'entête des colonnes.
Ensuite il y a deux lignes qui sont des informations sur les colonnes, soit par rapport au format, soit par rapport à l'unité, soit sur comment ont été obtenues les valeurs.
Les données commencent réellement à partir de la ligne 61 avec janvier 1958
data = read.csv("monthly_in_situ_co2_mlo.csv",skip=57)
head(data)
```
Dans un premier temps nous allons donc récuperer l'entête des colonnes, puis charger les données sans l'entête et renommer les colonnes conformément à l'entête.
Enfin, la description du fichier nous indique que les valeurs -99.99 correspondent à des données manquantes.
Cette conversion renvoie le 15 du mois ce qui est conforme à la description dans le fichier "The monthly values have been adjusted to 24:00 hours on the 15th of each month"
## Représentations graphique
### Réprésentation de l'évolution de la concentration de CO2 atmosphérique en fonction du temps
D'après la description des données contenu dans le fichier csv la colonne 5 CO2 correspond à la concentration
en CO2 : "Column 5 below gives monthly Mauna Loa CO2 concentrations in micro-mol CO2 per mole (ppm)".
Ce ne sont pas des données brutes mais des données recalculées à une échelle mensuelle.
mtext("Mauna Loa Observatory, Hawaii", side=3,line=2,cex=1.25,adj=0)
mtext("Monthly average carbone dioxide concentration", side=3,line=1,cex=1.25,adj=0)
mtext("Data from Scripps CO2 Program", side=3,lin=0,cex=.75,adj=0)
```
Nous voyons clairement deux structures dans les données, une forte tendance à l'augmentation au cours du temps et une oscillation des valeurs autour de cette tendance à une échelle intra-annuelle.
Ceci est un document R markdown que vous pouvez aisément exporter au format HTML, PDF, et MS Word. Pour plus de détails sur R Markdown consultez <http://rmarkdown.rstudio.com>.
## Caractérisation de l'oscillation périodique et de la tendance
Lorsque vous cliquerez sur le bouton **Knit** ce document sera compilé afin de ré-exécuter le code R et d'inclure les résultats dans un document final. Comme nous vous l'avons montré dans la vidéo, on inclue du code R de la façon suivante:
### Transformation des données en série temporelle
```{r cars}
summary(cars)
Pour caractériser ces deux composantes, nous allons considérer les données comme une série temporelle.
Nous allons transformer les données en un objet ts. Pour cela nous devons utiliser une série de données sans valeurs manquantes.
Plutôt que la colonne 5 nous allons donc utiliser la colonne 9 "CO2.1" qui correspond aux mêmes données que la colonne 5
mais où les données manquantes ont été remplacées en utilisant une "smoothed version of the data generated from a stiff cubic spline function plus 4-harmonic functions with linear gain". Cette procédure est décrite en détail dans le rapport cité en introduction de ce document.
Nous pouvons le vérifier en superposant les données remplies et les données avec données manquantes :
mtext("Mauna Loa Observatory, Hawaii", side=3,line=2,cex=1.25,adj=0)
mtext("Monthly average carbone dioxide concentration", side=3,line=1,cex=1.25,adj=0)
mtext("Data from Scripps CO2 Program", side=3,lin=0,cex=.75,adj=0)
points(data$Date,data$CO2,type="l",col=1)
legend("topleft",c("CO2 data with NA","CO2 data filled"),lty=1,col=c(1,2),bty="n")
```
Et on peut aussi aisément inclure des figures. Par exemple:
Nous constatons qu'il y a très peu de données manquantes et qu'elles se situent toutes entre 1958 et 1964 (nous les avions déjà identifiées plus haut).
mtext("Mauna Loa Observatory, Hawaii", side=3,line=2,cex=1.25,adj=0)
mtext("Monthly average carbone dioxide concentration", side=3,line=1,cex=1.25,adj=0)
mtext("Data from Scripps CO2 Program", side=3,lin=0,cex=.75,adj=0)
points(data$Date,data$CO2,type="l",col=1)
legend("topleft",c("CO2 data with NA","CO2 data filled"),lty=1,col=c(1,2),bty="n")
```
Vous remarquerez le paramètre `echo = FALSE` qui indique que le code ne doit pas apparaître dans la version finale du document. Nous vous recommandons dans le cadre de ce MOOC de ne pas utiliser ce paramètre car l'objectif est que vos analyses de données soient parfaitement transparentes pour être reproductibles.
Il nous faut également supprimer les premières et dernières données (manquantes)
allant de janvier 1958 à février 1958 puis de juin 2023 à décembre 2023.
```{r}
data.trunc<-data[!is.na(data$CO2.1),]
head(data.trunc[,c("Date","CO2","CO2.1")],10)
tail(data.trunc[,c("Date","CO2","CO2.1")],10)
```
Nous pouvons maintenant transformer la colonne 9 en un objet ts.
Nous commençons par trier les données pour être sûr que les lignes sont bien ordonnées selon la date.
La série temporelle correspond à des données mensuelles sur plusieurs années.
La première donnée est le mois de mars 1958. La fin est le mois de mai 2023.
Nous créons notre objet ts en indiquant le début et la fréquence de la série.
Comme les résultats ne sont pas stockés dans les fichiers Rmd, pour faciliter la relecture de vos analyses par d'autres personnes, vous aurez donc intérêt à générer un HTML ou un PDF et à le commiter.
Maintenant, à vous de jouer! Vous pouvez effacer toutes ces informations et les remplacer par votre document computationnel.
### Décomposition de la série temporelle
La fonction stl du package stats permet de décomposer la série temporelle en utilisant un modèle loess pour ajuster la tendance.
Nous avons fixé le paramètre s.window à 13 ce qui signifie que la fonction utilise une fenêtre de 13 annés consécutives pour estimer chaque valeur de la composante saisonnière. Nous avons préféré cela plutôt que la valeur par défaut ("periodic") qui considère que la composante saisonnière ne varie pas au cours du temps.
La fonction renvoie la composante saisonnière, la tendance et les résidus.
```{r}
require(stats)
#decom.stl<-stl(data.ts,s.window="periodic")
decom.stl<-stl(data.ts,s.window=13)
summary(decom.stl)
plot(decom.stl)
```
### Caractérisation de la saisonalité
Nous pouvons récupérer les différentes composantes de la série temporelle et représenter la seasonalité seule
mtext("Variation saisonnière de la concentration en CO2 (ppm)",line=1)
mtext("hors tendance",line=0)
```
Nous pouvons également donner quelques statistiques de cette composante saisonnières.
Dans une année le CO2 atmosphérique peut varier de `r round(min(data.trunc$seasonSTL),2)` ppm à `r round(max(data.trunc$seasonSTL),2)` ppm par rapport à la moyenne annuelle. L'écart-type de la composante saisonnière est de `r round(sd(data.trunc$seasonSTL),2)` ppm.
Nous voyons sur le graphe précédent que cette variation intra-annuelle augmente au cours du temps.
La variance de la composante saisonnière n'est pas constante au cours du temps.
### Caractérisation de la tendance
Nous pouvons superposer la tendance aux données brutes ce qui correspond au même type de données que la colonne 6 seasonally : "Column 6 gives the same data after a seasonal adjustment to remove the quasi-regular seasonal cycle"
ce modèle est simpliste car par exemple nous n'avons pas tenu compte des observations répétées par années, c'est à dire de la présence de pseudoréplications, ni de l'autocorrélation temporelle.
Néanmoins, même s'il n'est pas parfait, ce modèle très simple semble bien convenir aux données.
### Prédiction jusqu'en 2025 à l'aide du modèle simple
legend("topleft",c("tendance observée","valeurs prédites","intervalle de confiance des prédictions"),col=c("red","black","black"),lty=c(1,1,2),bty="n",cex=0.75)
```
```{r}
modeleSimple0<-lm(data=data.trunc, CO2.1 ~ time + I((time)^2))
