Copie du rapport computationnel

module3/exo3/exercice_R_fr.org

-#+TITLE:  Votre titre
-#+AUTHOR: Votre nom
-#+DATE:   La date du jour
+#+TITLE: Rapport computationnel sur le syesujet 1
+#+AUTHOR: L. SICOT
 #+LANGUAGE: fr
-# #+PROPERTY: header-args :eval never-export
 
 #+HTML_HEAD: <link rel="stylesheet" type="text/css" href="http://www.pirilampo.org/styles/readtheorg/css/htmlize.css"/>
 #+HTML_HEAD: <link rel="stylesheet" type="text/css" href="http://www.pirilampo.org/styles/readtheorg/css/readtheorg.css"/>
@@ -11,74 +9,261 @@
 #+HTML_HEAD: <script type="text/javascript" src="http://www.pirilampo.org/styles/lib/js/jquery.stickytableheaders.js"></script>
 #+HTML_HEAD: <script type="text/javascript" src="http://www.pirilampo.org/styles/readtheorg/js/readtheorg.js"></script>
 
-* Quelques explications
+# L'option ci-dessous évite de mettre en indice les caractères suivant
+# les underscores des noms de fichiers.
+# On utilisera _{xxx} pour mettre xxx en indice.  
+#+OPTIONS: ^:{}
 
-Ceci est un document org-mode avec quelques exemples de code
-R. Une fois ouvert dans emacs, ce document peut aisément être
-exporté au format HTML, PDF, et Office. Pour plus de détails sur
-org-mode vous pouvez consulter https://orgmode.org/guide/.
 
-Lorsque vous utiliserez le raccourci =C-c C-e h o=, ce document sera
-compilé en html. Tout le code contenu sera ré-exécuté, les résultats
-récupérés et inclus dans un document final. Si vous ne souhaitez pas
-ré-exécuter tout le code à chaque fois, il vous suffit de supprimer
-le # et l'espace qui sont devant le ~#+PROPERTY:~ au début de ce
-document.
+* Nom de la fiche
+Sujet 1 : Concentration de CO_{2} dans l'atmosphère depuis 1958
 
-Comme nous vous l'avons montré dans la vidéo, on inclut du code
-R de la façon suivante (et on l'exécute en faisant ~C-c C-c~):
+* Contexte
+Plusieurs autres fiches étaient intéressantes et m'auraient permis de découvrir
+de nouveaux packages, fonctions ou méthodologies mais ne pouvant consacrer
+autant de temps que je l'aurais voulu à cet exercice, j'ai choisi le sujet 1
+qui semblait se rapprocher le plus de ce que je connaissais.
 
-#+begin_src R :results output :exports both
-print("Hello world!")
-#+end_src
+* Objectifs
+Les trois objectifs de cette fiche sont les suivants :
+1. Réalisez un graphique qui vous montrera une oscillation périodique
+   superposée à une évolution systématique plus lente.
+2. Séparez ces deux phénomènes. Caractérisez l'oscillation périodique.
+   Proposez un modèle simple de la contribution lente, estimez ses paramètres
+   et tentez une extrapolation jusqu'à 2025 (dans le but de pouvoir valider
+   le modèle par des observations futures).
+3. Déposer dans FUN votre résultat
+
+* Source du fichier
+- Les données utilisées pour cette fiche sont issues du site
+[[https://scrippsco2.ucsd.edu/data/atmospheric_co2/mlo.html]]
+- Le nom du fichier, téléchargé le 27/11/2022, est :
+weekly_in_situ_co2_mlo.csv
+
+* Analyse
+Dans un premier temps, on charge les librairies dont nous aurons besoin :
+  - lubridate pour la gestion des dates
+  - MASS pour la fonction de régression rlm
+  - stplus pour la décomposition des effets saisonniers
+
+On ouvre ensuite le fichier avec readLines pour comprendre
+sa structure. Pour supprimer les premières lignes, on cherche la
+position de la date de départ avec grep. La valeur de la ligne est stockée
+dans la variable start_row.
+
+On peut ensuite ouvrir à nouveau le fichier en se limitant aux données qui
+nous intéressent, renommer les colonnes et convertir la colonne Date en
+classe Date.
+   
+#+BEGIN_SRC R :session :results none :exports both
+  library(lubridate)
+  library(MASS)
+  library(stlplus)
+  data_raw <- readLines(
+	 "weekly_in_situ_co2_mlo_20221127.csv")
+  start_row <- grep("1958-03-29", data_raw)
+  data_co2 <- read.table(
+	"weekly_in_situ_co2_mlo_20221127.csv",
+    sep = ",", dec = ".", as.is = T, header = F, skip = (start_row - 1))
+  names(data_co2) <- c("Date", "MicromolCO2")
+  data_co2$Date <- as.Date(data_co2$Date)
+#+END_SRC
+
+Les données se présentent sous deux colonnes : la date et la mesure de CO_{2}.
+#+BEGIN_SRC R :session :results output  :exports both
+head(data_co2)
+tail(data_co2)
+#+END_SRC
 
 #+RESULTS:
-: [1] "Hello world!"
+#+begin_example
+        Date MicromolCO2
+1 1958-03-29      316.19
+2 1958-04-05      317.31
+3 1958-04-12      317.69
+4 1958-04-19      317.58
+5 1958-04-26      316.48
+6 1958-05-03      316.95
+           Date MicromolCO2
+3294 2022-09-24      414.82
+3295 2022-10-01      415.12
+3296 2022-10-08      414.85
+3297 2022-10-15      415.31
+3298 2022-10-22      415.60
+3299 2022-10-29      416.08
+#+end_example
 
-Voici la même chose, mais avec une session R (c'est le cas le
-plus courant, R étant vraiment un langage interactif), donc une
-persistance d'un bloc à l'autre (et on l'exécute toujours en faisant
-~C-c C-c~).
+Comme vous téléchargerez votre fichier plus tard que moi, pour avoir les
+mêmes résultats, pensez à filtrer les dates pour avoir un jeu de données
+qui se termine le 29/10/2022.
 
-#+begin_src R :results output :session *R* :exports both
-summary(cars)
-#+end_src
+Les mesures sont normalement hebdomadaires. Pour vérifier que c'est bien le
+cas, on peut calculer l'écart entre chaque date 
+#+BEGIN_SRC R :session :results output  :exports both
+table(difftime(tail(as.Date(data_co2$Date), -1), head(as.Date(data_co2$Date), -1)))
+#+END_SRC
 
 #+RESULTS:
-:      speed           dist       
-:  Min.   : 4.0   Min.   :  2.00  
-:  1st Qu.:12.0   1st Qu.: 26.00  
-:  Median :15.0   Median : 36.00  
-:  Mean   :15.4   Mean   : 42.98  
-:  3rd Qu.:19.0   3rd Qu.: 56.00  
-:  Max.   :25.0   Max.   :120.00
-
-Et enfin, voici un exemple de sortie graphique:
-#+begin_src R :results output graphics :file "./cars.png" :exports results :width 600 :height 400 :session *R* 
-plot(cars)
+: 
+:    7   14   21   28   35   42   63  133 
+: 3269   17    5    2    2    1    1    1
+
+On note que plusieurs mesures sont manquantes. Pour conserver un écart de 7
+jours, nous allons ajouter des dates avec une mesure NA.
+
+#+BEGIN_SRC R :session :results output  :exports both
+  range_date <-range(data_co2$Date)
+  seq_date <- seq(range_date[1], range_date[2], by = "week")
+  df_na <- data.frame(Date = seq_date,
+		      MicromocCO2 = rep(NA, length(seq_date)))
+  data_co2 <- merge(df_na, data_co2, all.x = TRUE, by = 'Date')[, c(1, 3)]
+#+END_SRC
+
+On vérifie que les dates manquantes sont bien remplacées par une date
+avec une mesure NA.
+#+BEGIN_SRC R :session :results output  :exports both
+head(data_co2, 15)
+#+END_SRC
+
+Le graphe résultant, qui répond à l'objectif 1, montre bien la superposition
+d'une oscillation périodique à une évolution systématique plus lente :
+#+BEGIN_SRC R :exports both :session :results output graphics file :file co2_1_plot.png
+    plot(data_co2$Date, data_co2$MicromolCO2,
+	 type = "l", las = 1, 
+	 xlab = "", ylab = expression(paste("Concentration de ",CO[2]," dans l'atmosphère (ppm)")))
 #+end_src
 
+En zoomant, on peut observer les données manquantes :
+
+#+BEGIN_SRC R :exports both :session :results output graphics file :file co2_2_plot.png
+	    plot(data_co2$Date, data_co2$MicromolCO2,
+		 type = "l", las = 1,
+		 xlab = "",
+		 ylab = expression(paste("Concentration de ",CO[2]," dans l'atmosphère (ppm)")),
+		 ylim = c(311, 329), xlim = c(-4683.2972, -146.7261))
+#+end_src
+
+Afin de décomposer la partie oscillatoire de la tendance globale, on utilise
+la fonction stlplus du package éponyme. La fonction a besoin d'un objet de
+classe ts ou un vecteur de données numériques (concentration en CO_{2}) et
+un vecteur avec les dates à considérer (dates hebdomadaires).
+
+J'ai aussi essayé d'utiliser la fonction ts mais les arrondis faisaient que les
+dates générées ne correspondaient pas exactement aux dates de la colonne
+data_co2$Date.
+
+Pour générer le vecteur temps, j'ai donc remplacé le code ci-dessous :
+#+BEGIN_EXAMPLE
+  c_date <- ts(data_co2$MicromolCO2, 
+   freq=365.25/7, 
+   start=decimal_date(ymd("1958-03-29"))).
+#+END_EXAMPLE
+par celui-ci qui utilise des fonctions du package lubridate :
+#+BEGIN_SRC R :session :results none :exports both
+c_date <- decimal_date(ymd(data_co2$Date))
+#+end_src
+
+On utilise ensuite ce vecteur avec la fonction stlplus. J'ai dû jouer avec les
+paramètres, en particulier s.window, afin de réduire les résidus sans faire
+apparaître de schémas répétitifs et obtenir la décomposition ci-dessous.
+Le paramètre n.p correspond à la périodicité des données (52 valeurs par an).
+
+Le premier graphe correspond aux données brutes. Le deuxième et le troisième graphes
+sont la réponse à la première partie de l'objectif 2, à savoir la décomposition
+de la partie oscillatoire périodique et de la tendance à long terme.
+Le dernier graphe correspond aux résidus entre les données brutes et l'ajustement
+de l'effet saisonnier et la tendance.
+
+#+BEGIN_SRC R :exports both :session :results output graphics file :file decompose.png
+data_decompose <- stlplus(data_co2$MicromolCO2, c_date,
+n.p = 52, s.window = 6, s.degree = 2)
+plot(data_decompose, scales = list(relation = "free"))
+#+end_src
+
+De manière à trouver la période du signal oscillant, on cherche les changements
+de signe du signal. Chaque changement de signe indique un passage à 0 du signal.
+La différence de temps entre deux passages à 0 donne la demi-période du signal.
+Il suffit de doubler la médiane des demi-périodes pour avoir une valeur médiane
+de la période du signal oscillant.
+#+BEGIN_SRC R :session :results output  :exports both
+  df_season <- data.frame(Date = c_date,
+			  Signal = data_decompose$data$seasonal)
+  c_sign <- sign(head(df_season$Signal, -1) * tail(df_season$Signal, -1))
+  c_sign <- c(c_sign[1], c_sign)
+  row_minus <- which(c_sign < 0)
+  print(paste(
+    "La période de l'oscillation est de ",
+    round(median(diff(df_season$Date[row_minus]))*2, 2),
+    " année.", sep = ""))
+#+END_SRC
+
+#+RESULTS:
+: [1] "La période de l'oscillation est de 0.96 année."
+
+Pour déterminer l'amplitude du signal, on découpe le signal par période et,
+pour chaque période, on recherche les valeurs min et max.
+On donne ensuite les valeurs médianes min et max de l'amplitude.
+#+BEGIN_SRC R :session :results output  :exports both
+  cut_date <- cut(df_season$Date,
+			breaks = df_season$Date[row_minus[seq(1, length(row_minus), by = 2)]])
+  interval_date <- split(df_season, cut_date)
+  df_min_max <- do.call("rbind",
+			lapply(interval_date,
+			       function(X){ data.frame(Date = X$Date[1],
+						       val_min = min(X$Signal),
+						       val_max = max(X$Signal))}))
+  med_min <- median(df_min_max$val_min)
+  med_max <- median(df_min_max$val_max)
+  print(paste("La valeur minimale médiane du signal oscillant est de ",
+	      round(med_min, 2), " ppm.", sep = ""))
+  print(paste("La valeur maximale médiane du signal oscillant est de ",
+	      round(med_max, 2), " ppm.", sep = ""))
+#+END_SRC
+
+#+RESULTS:
+: [1] "La valeur minimale médiane du signal oscillant est de -3.68 ppm."
+: [1] "La valeur maximale médiane du signal oscillant est de 3.14 ppm."
+Le signal oscillant a donc une période proche d'un an et est relativement symétrique.
+
+Le dernier objectif consiste à modéliser la contribution lente et à proposer
+une extrapolation jusqu'en 2025.
+Après avoir crée un data.frame avec la contribution lente, on ajuste les données
+avec un polynôme d'ordre 2 (pour être tout à fait honnête, j'ai trouvé ce modéle
+sur internet à cette adresse : [[https://mgimond.github.io/ES218/Week10a.html]]).
+
+On affiche ensuite les paramètres du modèle et la prévision pour 2025.
+#+BEGIN_SRC R :session :results output  :exports both
+  df_trend <- data.frame(Date = c_date,
+			 Trend = data_decompose$data$trend)			  
+  fit_trend <- rlm(Trend ~ Date + I(Date^2), df_trend)
+  print(fit_trend)
+#+END_SRC
+
+D'après le modèle, la prévision de concentration de CO_{2} en ppm pour 2025 sera
+de :
+#+BEGIN_SRC R :session :results output  :exports both
+  trend_2025 <- round(predict(fit_trend, data.frame(Date = 2025)) ,1)
+  print(trend_2025[[1]])
+#+END_SRC
+
+#+RESULTS:
+: [1] 424.2
+
+Ci-dessous, la représentation graphique de la courbe de tendance, du modèle
+et la prévision en 2025.
+#+BEGIN_SRC R :exports both :session :results output graphics file :file model_trend.png
+  plot(df_trend, type = "l", las = 1, cex.axis = 1.2, cex.lab = 1.2,
+       xlab = "",
+       ylab = expression(paste(CO[2], " concentration trend  (ppm)")),
+       xlim = c(1956, 2026),
+       ylim = c(300, 450)
+       )
+  points(df_trend$Date, predict(fit_trend), type = "l", col ="red")
+  points (2025, trend_2025, pch = 20, cex = 1.6, col = "red")
+  legend(x = "topleft", c("Trend", "Model", "Forecast in 2025"),
+	 col = c("black", "red", "red"),
+	 pch = c(45, 45, 20), pt.cex = 2, cex = 1.2, inset = 0.02, bg = "white", border = T)
+#+END_SRC
+
 #+RESULTS:
-[[file:./cars.png]]
-
-Vous remarquerez le paramètre ~:exports results~ 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 changer ce paramètre
-(indiquer ~both~) car l'objectif est que vos analyses de données soient
-parfaitement transparentes pour être reproductibles.
-
-Attention, la figure ainsi générée n'est pas stockée dans le document
-org. C'est un fichier ordinaire, ici nommé ~cars.png~. N'oubliez pas
-de le committer si vous voulez que votre analyse soit lisible et
-compréhensible sur GitLab.
-
-Enfin, pour les prochains exercices, nous ne vous fournirons pas
-forcément de fichier de départ, ça sera à vous de le créer, par
-exemple en repartant de ce document et de le commiter vers
-gitlab. N'oubliez pas que nous vous fournissons dans les ressources de
-ce MOOC une configuration avec un certain nombre de raccourcis
-claviers permettant de créer rapidement les blocs de code R (en
-faisant ~<r~ ou ~<R~ suivi de ~Tab~).
-
-Maintenant, à vous de jouer! Vous pouvez effacer toutes ces
-informations et les remplacer par votre document computationnel.
+[[file:model_trend.png]]