diff --git a/exercice_R_fr.org b/exercice_R_fr.org deleted file mode 100644 index bde9c941de60105b2787d9707c58ba934344d385..0000000000000000000000000000000000000000 --- a/exercice_R_fr.org +++ /dev/null @@ -1,269 +0,0 @@ -#+TITLE: Rapport computationnel sur le sujet 1 -#+AUTHOR: L. SICOT -#+LANGUAGE: fr - -#+HTML_HEAD: -#+HTML_HEAD: -#+HTML_HEAD: -#+HTML_HEAD: -#+HTML_HEAD: -#+HTML_HEAD: - -# 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: ^:{} - - -* Nom de la fiche -Sujet 1 : Concentration de CO_{2} dans l'atmosphère depuis 1958 - -* 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. - -* 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( - "C:/Disque D/Org/weekly_in_situ_co2_mlo_20221127.csv") - start_row <- grep("1958-03-29", data_raw) - data_co2 <- read.table( - "C:/Disque D/Org/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: -#+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 - -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. - -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: -: -: 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:model_trend.png]]