From 19dfed82424fba956d1fbb9d6ac30ac2285b00cb Mon Sep 17 00:00:00 2001 From: 42eaffb07ab321885c173115fa1a5bf9 <42eaffb07ab321885c173115fa1a5bf9@app-learninglab.inria.fr> Date: Sun, 4 Dec 2022 09:35:39 +0000 Subject: [PATCH] Copie du rapport computationnel --- module3/exo3/exercice_R_fr.org | 309 ++++++++++++++++++++++++++------- 1 file changed, 247 insertions(+), 62 deletions(-) diff --git a/module3/exo3/exercice_R_fr.org b/module3/exo3/exercice_R_fr.org index 1bb8f61..a4c8375 100644 --- a/module3/exo3/exercice_R_fr.org +++ b/module3/exo3/exercice_R_fr.org @@ -1,8 +1,6 @@ -#+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: #+HTML_HEAD: @@ -11,74 +9,261 @@ #+HTML_HEAD: #+HTML_HEAD: -* 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 ~