@@ -25,6 +25,8 @@ Dans ce document nous allons présenter et analyser des données d'évolution de
...
@@ -25,6 +25,8 @@ Dans ce document nous allons présenter et analyser des données d'évolution de
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 données sont disponibles sur le site suivant : [Scripps CO2 Program](https://scrippsco2.ucsd.edu/data/atmospheric_co2/primary_mlo_co2_record.html).
Nous avons téléchargé les données le 12 juin 2023.
Les premières 57 lignes du fichier donnent des informations sur les données.
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 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 23 à 37 indiquent comment citer les données.
...
@@ -107,20 +109,20 @@ Nous voyons clairement deux structures dans les données, une forte tendance à
...
@@ -107,20 +109,20 @@ Nous voyons clairement deux structures dans les données, une forte tendance à
### Transformation des données en série temporelle
### 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.
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.
Nous allons transformer les données en un objet ts (time series).
Pour cela nous devons utiliser une série de données sans valeurs manquantes.
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".
Plutôt que la colonne 5 nous allons donc utiliser la colonne 9 "CO2.1" décrite dans le fichier csv et 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.
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 :
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("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("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)
mtext("Data from Scripps CO2 Program", side=3,lin=0,cex=.75,adj=0)
points(data$Date,data$CO2,type="l",col=1)
points(data$Date,data$CO2,type="l",col=1,lwd=2)
legend("topleft",c("CO2 data with NA","CO2 data filled"),lty=1,col=c(1,2),bty="n")
legend("topleft",c("CO2 data with NA","CO2 data filled"),lty=1,col=c(1,"red"),lwd=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).
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).
...
@@ -210,7 +212,7 @@ Pour mieux voir cet effet saisonnier nous pouvons le représenter en superposant
...
@@ -210,7 +212,7 @@ Pour mieux voir cet effet saisonnier nous pouvons le représenter en superposant
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).
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).
...
@@ -260,16 +262,15 @@ legend("topleft",c("CO2 concentration","Column seasonally from file",
...
@@ -260,16 +262,15 @@ legend("topleft",c("CO2 concentration","Column seasonally from file",
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.
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).
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 simplifier nous allons ajuster ce modèle sur la tendance extraite précédemment avec la fonction stl(), variable nommée trendSTL dans le data frame.
Pour ajuster le modèle nous allons convertir le champ date en numérique et utiliser cette variable comme variable explicative du modèle.
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".
Nous nommons cette variable "time".
```{r}
```{r}
data.trunc$time<-as.numeric(data.trunc$Date)
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.
Nous ajustons le modèle et récupérons les valeurs ajustées et les résidus que nous représentons.
```{r}
```{r}
modeleSimple0<-lm(data=data.trunc, trendSTL ~ time + I((time)^2))
modeleSimple0<-lm(data=data.trunc, trendSTL ~ time + I((time)^2))
Nous pouvons estimer à l'aide du modèle simple la concentration moyenne attendue en CO2 pour l'année 2025.
Nous pouvons estimer à l'aide du modèle simple la concentration moyenne attendue en CO2 pour l'année 2025. Pour cela nous récupérons dans le data frame des extrapolations les données de l'année 2025 et nous calculons quelques statistiques pour l'année 2025.
**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.**
**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
### Modélisation à l'aide d'un modèle SARIMA
Pour aller un peu plus loin nous allons tenter d'utiliser un modèle stochastique de type SARIMA pour faire les prédictions.
Pour aller un peu plus loin nous allons tenter d'utiliser un modèle stochastique de type SARIMA ("Seasonal Autoregressive Integrated Moving Average") 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é.
En effet, nous avons vu précédemment que la série temporelle présentait à 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.
Nous pouvons également vérifier que les résidus obtenus avec la fonction stl() présentent une autocorrélation et une autocorrélation partielle à l'aide des fonctions acf et pcaf.
```{r}
```{r}
par(mfrow=c(1,2))
par(mfrow=c(1,2))
...
@@ -359,13 +362,13 @@ pacf(data.trunc$residSTL)
...
@@ -359,13 +362,13 @@ 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.
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. Même si cela rallonge le temps de calcul, nous fixons les paramètres stepwise et approximation à FALSE pour obtenir le meilleur modèle possible.
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.
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.
De la même façon que pour le modèle simple nous pouvons estimer la concentration de CO2 en 2025.
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.
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.
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.
Nous devons donc fixer la langue des dates en anglais pour récupérer la date dans le bon format puis la remettre en français.
```{r}
```{r}
Sys.setlocale("LC_TIME", "English")
Sys.setlocale("LC_TIME", "English")
...
@@ -411,7 +414,15 @@ L'intervalle de confiance des prédictions est légèrement plus large avec le m
...
@@ -411,7 +414,15 @@ L'intervalle de confiance des prédictions est légèrement plus large avec le m
## Conclusion
## 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.
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'Hawaii. 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.
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. L'avantage du modèle SARIMA est qu'il peut prédire également les variations saisonnières de la concentration de CO2 atmosphérique.
