diff --git a/module3/exo3/Module3Exo3Sujet1.Rmd b/module3/exo3/Module3Exo3Sujet1.Rmd new file mode 100644 index 0000000000000000000000000000000000000000..41299ce529d5049f18314cdacc675b5475e1ffa4 --- /dev/null +++ b/module3/exo3/Module3Exo3Sujet1.Rmd @@ -0,0 +1,417 @@ +--- +title: "Concentration de CO2 dans l'atmosphère depuis 1958" +author: "Thomas Pérot" +date: "12 juin 2023" +output: + pdf_document: + toc: true +editor_options: + markdown: + wrap: sentence +bibliography: references.bib +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +## 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ï. Ces données sont décrites dans le document suivant : + +*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 informations sur les données. +Les lignes 2 à 19 indiquent l'endroit d'où proviennent les 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 + +```{r} + +#data_url = "https://scrippsco2.ucsd.edu/assets/data/ +#atmospheric/stations/in_situ_co2/monthly/monthly_in_situ_co2_mlo.csv" +data_csv = "monthly_in_situ_co2_mlo.csv" +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. +Nous allons les remplacer par des NA et nous afficherons l'ensemble des valeurs manquantes. + +```{r} +header = read.csv("monthly_in_situ_co2_mlo.csv",skip=57,nrows=1) +data = read.csv("monthly_in_situ_co2_mlo.csv",skip=60,header=FALSE) +names(data)<-names(header) +data[data==-99.99]<-NA + +na_records = apply(data, 1, function (x) {any(is.na(x))}) +data[na_records,] + +``` + +### Conversion des colonnes en date + +**Problème rencontré** : les colonnes 3 et 4 sont des dates mais le format de la date n'est pas indiqué. +Après avoir recherché les différentes possibilités voici la description des quatre premières colonnes : + +- colonne 1 : année ; +- colonne 2 : mois ; +- colonne 3 : date en nombre de jours depuis une origine, la date d'origine étant 1899-12-30 (comme dans Excel ce qui était indiqué par le mot Excel ligne 60 du fichier) ; +- colonne 4 : date en année décimal. + +Nous allons convertir la colonne 3 en date en indiquant la date 1899-12-30 comme origine. + +```{r} +class(data$Date) +data$Date<-as.Date(data$Date, origin = "1899-12-30") +class(data$Date) +head(data[,c("Yr","Mn","Date","CO2")]) +tail(data[,c("Yr","Mn","Date","CO2")]) + +``` + +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". +Nous n'avons pas réussi à faire de même avec la colonne 4. + +## Représentations graphiques des données + +### 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 contenue 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)". +D'après le document cité en introduction, ce ne sont pas des données brutes mais des données recalculées à une échelle mensuelle. +Nous pouvons représenter ces données sur un graphe semblable à celui du site [Scripps CO2 Program](https://scrippsco2.ucsd.edu/data/atmospheric_co2/primary_mlo_co2_record.html). + +```{r} +plot(data$Date,data$CO2,type="l",main="",xlab="Date (month)", + ylab="CO2 concentration (ppm)") +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. + +## Caractérisation de l'oscillation périodique et de la tendance + +### Transformation des données en série temporelle + +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 : + +```{r} +plot(data$Date,data$CO2.1,type="l",main="",xlab="Date (month)", + ylab="CO2 concentration (ppm)",col=2) +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") +``` + +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). + +```{r} +plot(data$Date,data$CO2.1,type="l",main="",xlab="Date (month)", + ylab="CO2 concentration (ppm)",col=2, + xlim=as.Date(c("1958-01-01","1970-01-01")),ylim=c(300,340)) +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") +``` + +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. + +```{r} +data.trunc<-data.trunc[order(data.trunc$Date),] +data.ts<-ts(data=data.trunc$CO2.1,start=c(1958,3),frequency=12) +plot(data.ts) +``` + +### 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) +``` + +Nous pouvons récupérer les différentes composantes de la série temporelle et représenter les données avec la tendance. + +```{r} +data.trunc$seasonSTL<-as.numeric(decom.stl[1]$time.series[,1]) +data.trunc$trendSTL<-as.numeric(decom.stl[1]$time.series[,2]) +data.trunc$residSTL<-as.numeric(decom.stl[1]$time.series[,3]) + +plot(data.trunc$Date,data.trunc$CO2.1,type="l",main="",xlab="Date (month)", + ylab="CO2 concentration (ppm)") +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.trunc$Date,data.trunc$trendSTL,col="blue",type="l") +legend("topleft",c("CO2 concentration", + "Trend from stl function"),lty=1,col=c(1,"blue"),bty="n") +``` + +### Caractérisation de la saisonalité + +Nous pouvons également représenter la saisonnalité seule. + +```{r} +plot(data.trunc$Date,data.trunc$seasonSTL,type="l",xlab="Date (month)", + ylab="CO2 concentration (ppm)",ylim=c(-4,4)) +mtext("Variation saisonnière de la concentration en CO2 (ppm)",line=1) +mtext("hors tendance",line=0) +``` + +Nous pouvons 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. +La variance de la composante saisonnière n'est pas constante au cours du temps. Nous voyons sur le graphe précédent que cette variation intra-annuelle augmente légèrement au cours du temps. + +Pour mieux voir cet effet saisonnier nous pouvons le représenter en superposant toutes les années dans un graphe annuel avec la fonction ggseasonplot() du packages forecast. Pour mieux identifier où se situe les pics et les creux nous allons faire cette représentation à partir de la composante saisonnière obtenue avec la fonction stl(). + +```{r} +#install.packages("forecast") +library(forecast) +ggseasonplot(decom.stl[1]$time.series[,1]) +``` +Nous voyons qu'au sein d'une année le pic de concentration de CO2 a lieu vers le mois de mai et les valeurs les plus basses sont vers le mois de septembre octobre. Les données reflètent les variations de la concentration de CO2 atmosphérique dans l'hemisphère nord. Dans l'hemisphère nord, à partir du printemps, les plantes captent le CO2 atmosphérique expliquant une partie de ce cycle saisonnier (voir le document cité en introduction pour plus de détails). + +### Caractérisation de la tendance + +La tendance obtenue avec la fonction stl() nous permet de donner quelques chiffres sur l'évolution de la concentration en CO2 atmosphérique. +Entre 1958 et 2023 la concentration de CO2 atmosphérique est passé de `r round(data.trunc[data.trunc$Date==as.Date("1958-03-15"),]$trendSTL,2)` ppm +à `r round(data.trunc[data.trunc$Date==as.Date("2023-05-15"),]$trendSTL,2)` ppm, soit une augmentation de `r +round(100*(data.trunc[data.trunc$Date==as.Date("2023-05-15"),]$trendSTL-data.trunc[data.trunc$Date==as.Date("1958-03-15"),]$trendSTL)/data.trunc[data.trunc$Date==as.Date("1958-03-15"),]$trendSTL,1)`% en 65 ans. Cette augmentation est principalement lié à l'utilisation de carburant fossile par l'homme (voir le document cité en introduction pour plus de détails). + +La tendance obtenue avec la fonction stl() 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". Nous pouvons le vérifier en superposant les deux variables sur un même graphe. + +```{r} +plot(data.trunc$Date,data.trunc$CO2.1,type="l",main="",xlab="Date (month)", + ylab="CO2 concentration (ppm)") +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.trunc$Date,data.trunc$seasonally,col="red",type="l") +points(data.trunc$Date,data.trunc$trendSTL,col="blue",type="l") +legend("topleft",c("CO2 concentration","Column seasonally from file", + "Trend from stl() function"),lty=1,col=c(1,"red","blue"),bty="n") + +``` + +Cependant, en zoomant nous voyons que la tendance que nous obtenons avec la fonction stl() correspond à un lissage plus fort que la colonne 6 "seasonally". +Les méthodes utilisées ne sont pas les mêmes. + +```{r} + +plot(data.trunc$Date,data.trunc$CO2.1,type="l",main="",xlab="Date (month)", + ylab="CO2 concentration (ppm)",xlim=as.Date(c("1980-01-01","1990-01-01")), + ylim=c(330,360)) +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.trunc$Date,data.trunc$seasonally,col="red",type="l") +points(data.trunc$Date,data.trunc$trendSTL,col="blue",type="l") +legend("topleft",c("CO2 concentration","Column seasonally from file", + "Trend from stl function"),lty=1,col=c(1,"red","blue"),bty="n") + +``` + +## Prédiction de la concentration de CO2 jusqu'en 2025 + +### Modélisatin de la tendance avec un modèle simple + +Lorsque l'on zoome sur les données comme sur le graphe précédent, un modèle linéaire semble convenir pour représenter la tendance d'augmentation du CO2 dans le temps. Cependant en observant l'ensemble des données il apparaît clairement que la concentration de CO2 augmente de plus en plus vite. +Un premier modèle simple pourrait être d'ajuster un modèle linéaire ayant une composante linéaire et une composante quadratique (un polynôme de degré 2). +Pour simplifier nous allons ajuster ce modèle sur la tendance extraite précédemment avec la fonction stl(). +Pour ajuster le modèle nous allons convertir le champ date en numérique et utiliser cette variable comme variable explicative du modèle. +Nous nommons cette variable "time". + +```{r} +data.trunc$time<-as.numeric(data.trunc$Date) +head(data.trunc) +``` + +Nous ajustons le modèle et récupérons les valeurs ajustées et les résidus. + +```{r} +modeleSimple0<-lm(data=data.trunc, trendSTL ~ time + I((time)^2)) + +data.trunc$modeleSimple0Fit<-fitted(modeleSimple0) +data.trunc$modeleSimple0Res<-resid(modeleSimple0) + +par(mfrow=c(1,2)) +plot(data.trunc$Date,data.trunc$trendSTL,type="l",main="",xlab="Date (month)", + ylab="CO2 concentration (ppm)") +mtext("Modèle simple", side=3,line=2,cex=0.75,adj=0) +mtext("Valeur ajustée", side=3,line=1,cex=0.75,adj=0) +points(data.trunc$Date,data.trunc$modeleSimple0Fit,col="blue",type="l") +legend("topleft",c("tendance","tendance ~ time + time^2"),lty=1, + col=c("black","blue"),bty="n",cex=0.75) +plot(data.trunc$Date,data.trunc$modeleSimple0Res,type="l",main="", + xlab="Date (month)",ylab="CO2 concentration (ppm)") +mtext("Modèle simple", side=3,line=2,cex=0.75,adj=0) +mtext("Résidus", side=3,line=1,cex=0.75,adj=0) +abline(h=0) +``` + +Ce modèle est simpliste car, par exemple, nous n'avons pas tenu compte des observations répétées par année, 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 + +Pour prédire la concentration de CO2 en 2025, nous allons utiliser le modèle précédent et l'appliquer jusqu'en 2025. +Au delà de mai 2023, nous n'avons plus de données. Il s'agira d'extrapolations à l'aide du modèle simple. +Nous utilisons la fonction predict.lm() qui permet d'obtenir l'intervalle de confiance des prédictions. + +Nous commençons par construire un data.frame avec une colonne "time" sur l'ensemble de la période de mars 1958 à décembre 2025. +Puis nous utilisons ce data.frame pour prédire les valeurs de CO2 à l'aide du modèle simple et de la fonction predict.lm(). +Nous créons un data.frame spécifiquement pour les données au-delà de mai 2023 qui correspond aux valeurs prédites extrapolées. + +```{r} +data.trunc$time<-as.numeric(data.trunc$Date) +myNewData<-data.frame(time=as.numeric(seq.Date(from=as.Date("1958-03-15"), + to=as.Date("2025-12-15"),by="month"))) +predictions<-data.frame(predict.lm(modeleSimple0,myNewData, + interval="prediction",type="response")) +myNewData<-cbind(myNewData,predictions) +myNewDataExtrapolation<-myNewData[myNewData$time > as.numeric(as.Date("2023-05-15")),] +``` + +Le graphe suivant représente la tendance "observée" sur la période connue, les valeurs de la tendance obtenues par ajustement du modèle simple et enfin les valeurs de la tendance prédite au-delà de mai 2023. + +```{r} +plot(data.trunc$Date,data.trunc$trendSTL,type="l",main="",xlab="Date (month)", + ylab="CO2 concentration (ppm)",col="red",xlim=as.Date(c("2000-01-01","2026-01-01")), + ylim=c(360,440)) +points(myNewData$time,myNewData$fit,col=1,type="l") +points(myNewData$time,myNewData$lwr,col=1,type="l",lty=2) +points(myNewData$time,myNewData$upr,col=1,type="l",lty=2) +points(myNewDataExtrapolation$time,myNewDataExtrapolation$fit,col="blue",type="l") +points(myNewDataExtrapolation$time,myNewDataExtrapolation$lwr,col="blue",type="l",lty=2) +points(myNewDataExtrapolation$time,myNewDataExtrapolation$upr,col="blue",type="l",lty=2) +legend("topleft",c("tendance observée","valeurs ajustées (modèle simple)", + "valeurs prédites (modèle simple)","intervalle de confiance des prédictions"),col=c("red","black","blue","black"),lty=c(1,1,1,2),bty="n",cex=0.75) + +``` + +Nous pouvons estimer à l'aide du modèle simple la concentration moyenne attendue en CO2 pour l'année 2025. + +```{r} +myNewData2025<-myNewDataExtrapolation[myNewDataExtrapolation$time >= + as.numeric(as.Date("2025-01-01")) & +myNewDataExtrapolation$time <= as.numeric(as.Date("2025-12-31")),] +dim(myNewData2025) +round(mean(myNewData2025$fit),2) +round(mean(myNewData2025$fit)-mean(myNewData2025$lwr),2); +round(mean(myNewData2025$fit)-mean(myNewData2025$upr),2) +``` + +**D'après le modèle simple, en première approximation la valeur de concentration moyenne en CO2 en 2025 serait de `r round(mean(myNewData2025$fit),2)` ppm plus ou moins `r round(mean(myNewData2025$fit)-mean(myNewData2025$lwr),2)` ppm.** + +### Modélisation à l'aide d'un modèle stochastique + +Pour aller un peu plus loin nous allons tenter d'utiliser un modèle stochastique de type SARIMA pour faire les prédictions. +En effet, nous avons vu précédemment que les données présentaient à la fois une tendance et une saisonnalité. +Nous pouvons également vérifier que les résidus présentent une autocorrélation et une autocorrélation partielle à l'aide des fonctions acf et pcaf. + +```{r} +par(mfrow=c(1,2)) +acf(data.trunc$residSTL) +pacf(data.trunc$residSTL) + +``` + +Pour simplifier le processus de construction du modèle nous allons utiliser la fonction auto.arima() du package forecast qui permet de choisir le type de modèle approprié et d'estimer les différents paramètres du modèle stochastique. + +```{r} +#install.packages("forecast") +library(forecast) + +modelSarima<-auto.arima(data.ts,stepwise=FALSE) +summary(modelSarima) +``` + +Nous pouvons vérifier qu'il n'y a plus d'autocorrélation dans les résidus du modèle SARIMA. + +```{r} +par(mfrow=c(1,2)) +acf(resid(modelSarima)) +pacf(resid(modelSarima)) +``` + +Pour prédire la concentration de CO2 jusqu'en 2025 nous allons utiliser le modèle SARIMA et la fonction forecast() du package forecast. +Pour aller jusqu'à décembre 2025 nous allons faire des prédictions sur 31 valeurs supplémentaires (de juin 2023 à décembre 2025). + +```{r} +predictionsSARIMA<-data.frame(forecast(modelSarima,h=31)) +tail(predictionsSARIMA) +``` + +Nous pouvons représenter les prédictions à l'aide de la fonction plot.forecast() du package forecast qui représente les prédictions avec les intervalles de confiance à 80 et 95%. La fonction représente également les données qui ont servi à construire le modèle. Pour plus de lisibilité nous choisissons de ne représenter que les 50 dernières valeurs passées avec le paramètre include=50. + +```{r} +plot(forecast(modelSarima,h=31),include=50) +``` + +De la même façon que pour le modèle simple nous pouvons estimer la concentration de CO2 en 2025. +Nous récupérons les dates dans le data frame des prédictions du modèle SARIMA et nous créons un data.frame uniquement pour 2025. +Les dates sont contenues dans le nom des lignes avec un format utilisant une abréviation du mois en anglais suivi de l'année. +Nous devons donc fixer la langue en anglais pour récupérer la date puis la remettre en français. + +```{r} +Sys.setlocale("LC_TIME", "English") +predictionsSARIMA$Date<-as.Date(paste("15 ",row.names(predictionsSARIMA),sep=""),format="%d %b %Y") +Sys.setlocale("LC_TIME", "French") +predictionsSARIMA2025<-predictionsSARIMA[predictionsSARIMA$Date >= as.Date("2025-01-01") & +predictionsSARIMA$Date <= as.Date("2025-12-31"),] +``` + +**D'après le modèle SARIMA, la valeur de concentration moyenne en CO2 en 2025 serait de `r round(mean(predictionsSARIMA2025$Point.Forecast),2)` ppm plus ou moins `r round(mean(predictionsSARIMA2025$Point.Forecast)-mean(predictionsSARIMA2025$Lo.95),2)` ppm.** + +Le modèle SARIMA prédit donc une concentration de CO2 très proche, très légèrement inférieure, par rapport au modèle simple polynomial. +L'intervalle de confiance des prédictions est légèrement plus large avec le modèle SARIMA. + +## Conclusion + +Dans cette analyse nous avons chargé un jeu de données correspondant à l'évolution de la concentration en CO2 atmosphérique provenant du site de Mauna Loa sur l'île d'Hawai. Pour construire le jeu de données à partir du fichier csv nous avons exclu des lignes correspondant à des informations autres que les données. Nous avons également géré des données manquantes codées par une valeur numérique. Enfin nous avons convertit une colonne en un format date reconnu par R. Cette dernière étape a demandé un peu de temps car les informations sur le format des colonnes représentant la date n'étaient pas précises. +Une fois la préparation du jeu de données réalisée nous avons représenté les données et constaté la présence d'une tendance et d'une saisonnalité fortes. Nous avons caractérisé les oscillations périodiques et la tendance en considérant les données comme une serie temporelle. Pour cela nous avons utilisé les fonction ts() et stl() du package stats. Ensuite, pour prédire la concentration en CO2 jsuqu'en 2025, nous avons construit un premier modèle simple de la tendance à l'aide d'un modèle linéaire correspondant à un polynôme de degré 2 du temps. Nous avons ensuite proposé un modèle stochastique SARIMA plus complexe mais plus adapté aux données qui permet de prédire la tendance et la saisonnalité. Pour cela nous avons utilisé les fonctions auto.arima() et forecast() du package forecast. Les deux approches donnent des résultats très similaires probablement parceque nous avons fait des prédictions à très court terme. + + diff --git a/module3/exo3/Module3Exo3Sujet1.pdf b/module3/exo3/Module3Exo3Sujet1.pdf new file mode 100644 index 0000000000000000000000000000000000000000..d4ad74034d76660306969a33a15e9e394dfcb372 Binary files /dev/null and b/module3/exo3/Module3Exo3Sujet1.pdf differ diff --git a/module3/exo3/exercice_fr.Rmd b/module3/exo3/exercice_fr.Rmd index 563588c9ff467d309feafe0577e9fd62c7ad3844..1f59b55e998c456e5880282da27faae020b56a06 100644 --- a/module3/exo3/exercice_fr.Rmd +++ b/module3/exo3/exercice_fr.Rmd @@ -173,7 +173,7 @@ 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 +Nous pouvons récupérer les différentes composantes de la série temporelle et représenter la saisonnalité seule. ```{r} data.trunc$seasonSTL<-as.numeric(decom.stl[1]$time.series[,1]) @@ -186,11 +186,19 @@ 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. +Nous pouvons 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. +La variance de la composante saisonnière n'est pas constante au cours du temps. Nous voyons sur le graphe précédent que cette variation intra-annuelle augmente légèrement au cours du temps. + +Pour mieux voir cet effet saisonnier nous pouvons le représenter en superposant toutes les années dans un graphe annuel avec la fonction ggseasonplot() du packages forecast. Pour mieux identifier où se situe les pics et les creux nous allons faire cette représentation à partir de la composante saisonnière obtenue avec la fonction stl(). + +```{r} +#install.packages("forecast") +library(forecast) +ggseasonplot(decom.stl[1]$time.series[,1]) +``` +Nous voyons qu'au sein d'une année le pic de concentration de CO2 a lieu vers le mois de mai et les valeurs les plus basses sont vers le mois de septembre octobre. Les données reflètent les variations de la concentration de CO2 atmosphérique dans l'hemisphère nord. Dans l'hemisphère nord, à partir du printemps, les plantes captent le CO2 atmosphérique expliquant une partie de ce cycle saisonnier (voir le document cité en introduction pour plus de détails). ### Caractérisation de la tendance @@ -227,21 +235,28 @@ legend("topleft",c("CO2 concentration","Column seasonally from file", ``` +Nous voyons qu'entre 1958 et 2023 la concentration de CO2 atmosphérique est passé de `r round(data.trunc[data.trunc$Date==as.Date("1958-03-15"),]$trendSTL,2)` ppm +à `r round(data.trunc[data.trunc$Date==as.Date("2023-05-15"),]$trendSTL,2)` ppm, soit une augmentation de `r +round(100*(data.trunc[data.trunc$Date==as.Date("2023-05-15"),]$trendSTL-data.trunc[data.trunc$Date==as.Date("1958-03-15"),]$trendSTL)/data.trunc[data.trunc$Date==as.Date("1958-03-15"),]$trendSTL,1)`% en 65 ans. + ## Prédiction de la concentration de CO2 jusqu'en 2025 ### Modélisatin de la tendance avec un modèle simple -Lorsque l'on zoom sur les données, un modèle linéaire semble convenir pour représenter la tendance d'augmentation du CO2 dans le temps. -Cependant en observant l'ensemble des données il apparaît clairement que la concentration de CO2 augmente de plus en plus vite. +Lorsque l'on zoom sur les données comme sur le graphe précédent, un modèle linéaire semble convenir pour représenter la tendance d'augmentation du CO2 dans le temps. Cependant en observant l'ensemble des données il apparaît clairement que la concentration de CO2 augmente de plus en plus vite. Un premier modèle simple pourrait être d'ajuster un modèle linéaire ayant une composante linéaire et une composante quadratique. -Pour simplifier nous allons ajuster ce modèle sur la tendance extraite précédemment. +Pour simplifier nous allons ajuster ce modèle sur la tendance extraite précédemment avec la fonction stl(). Pour ajuster le modèle nous allons convertir le champ date en numérique et utiliser cette variable comme variable explicative du modèle. Nous nommons cette variable "time". ```{r} data.trunc$time<-as.numeric(data.trunc$Date) head(data.trunc) +``` + +Nous ajustons le modèle et récupérons les valeurs ajustées et les résidus. +```{r} modeleSimple0<-lm(data=data.trunc, trendSTL ~ time + I((time)^2)) data.trunc$modeleSimple0Fit<-fitted(modeleSimple0) @@ -262,18 +277,17 @@ mtext("Résidus", side=3,line=1,cex=0.75,adj=0) abline(h=0) ``` -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. +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 Pour prédire la concentration de CO2 en 2025, nous allons utiliser le modèle précédent et l'appliquer jusqu'en 2025. -Au delà de mai 2023, nous n'avons plus de données. -Il s'agira d'extrapolations à l'aide du modèle simple. -Nous utilisons la fonction predict.lm qui permet d'obtenir l'intervalle de confiance des prédictions. +Au delà de mai 2023, nous n'avons plus de données. Il s'agira d'extrapolations à l'aide du modèle simple. +Nous utilisons la fonction predict.lm() qui permet d'obtenir l'intervalle de confiance des prédictions. Nous commençons par construire un data.frame avec une colonne "time" sur l'ensemble de la période de mars 1958 à décembre 2025. -Puis nous utilisons ce data.frame pour prédire les valeurs de CO2 à l'aide du modèle simple. +Puis nous utilisons ce data.frame pour prédire les valeurs de CO2 à l'aide du modèle simple et de la fonction predict.lm(). Nous créons un data.frame spécifiquement pour les données au-delà de mai 2023 qui correspond aux valeurs prédites extrapolées. ```{r} @@ -348,16 +362,15 @@ acf(resid(modelSarima)) pacf(resid(modelSarima)) ``` -Pour prédire la concentration de CO2 jusqu'en 2025 nous allons utiliser la fonction forecast() du package forecast. -Pour aller jusqu'à décembre 2025 nous allons prédire sur 31 valeur supplémentaires (de juin 2023 à décembre 2025). +Pour prédire la concentration de CO2 jusqu'en 2025 nous allons utiliser et modèle précédent et la fonction forecast() du package forecast. +Pour aller jusqu'à décembre 2025 nous allons prédire sur 31 valeurs supplémentaires (de juin 2023 à décembre 2025). ```{r} predictionsSARIMA<-data.frame(forecast(modelSarima,h=31)) tail(predictionsSARIMA) -names(predictionsSARIMA) ``` -Nous pouvons représenter les prédictions à l'aide de la fonction plot.forecast() du package forecast qui représente les prédictions avec les intervalles de confiance à 80 et 95%. +Nous pouvons représenter les prédictions à l'aide de la fonction plot.forecast() du package forecast qui représente les prédictions avec les intervalles de confiance à 80 et 95%. La fonction représente également les données qui ont servi à construire le modèle. Pour plus de lisibilité nous choisissons de ne représentées que les 50 dernières valeurs passées avec le paramètre include=50. ```{r} plot(forecast(modelSarima,h=31),include=50) @@ -379,3 +392,9 @@ predictionsSARIMA$Date <= as.Date("2025-12-31"),] Le modèle SARIMA prédit donc une concentration de CO2 très proche par rapport au modèle simple, très légèrement inférieure. L'intervalle de confiance des prédictions est légèrement plus large avec le modèle SARIMA. +## Conclusion + +Dans cette analyse nous avons chargé un jeu de données correspondant à l'évolution de la concentration en CO2 atmosphérique sur le site de Mauna Loa sur l'île d'Hawai. Pour construire le jeu de données à partir du fichier csv nous avons exclu des lignes correspondant à des informations autres que les données. Nous avons également géré des données manquantes codées par une valeur numérique. Enfin nous avons convertit une colonne en un format date reconnu par R. Cette dernière étape a demandé un peu de temps car les informations sur le format des colonnes représentant la date n'était pas précise. +Une fois la préparation du jeu de données réalisée nous avons représenté les données et constaté la présence d'une tendance et d'une saisonnalité fortes. Nous avons caractérisé les oscillations périodiques et la tendance en considérant les données comme une serie temporelle. Pour cela nous avons utilisé les fonction ts() et stl() du package stats. Ensuite, pour prédire la concentration en CO2 jsuqu'en 2025, nous avons construit un premier modèle simple de la tendance à l'aide d'un modèle linéaire correspondant à un polynôme de degré 2 du temps. Nous avons ensuite proposé un modèle stochastique SARIMA plus complexe mais plus adapté aux données qui permet de prédire la tendance et la saisonnalité. Les deux approches donnent des résultats très similaires probablement parceque nous avons fait des prédictions à très court terme. + +