diff --git a/module3/exo3/exercice_fr.Rmd b/module3/exo3/exercice_fr.Rmd index b97808914d6a6fb6529b56299236099867364cad..832bbb2149df12167e44ad365b294f16657bed83 100644 --- a/module3/exo3/exercice_fr.Rmd +++ b/module3/exo3/exercice_fr.Rmd @@ -18,9 +18,11 @@ Nous les récupérons sous forme d'un fichier en format CSV dont chaque ligne co En pratique, nous chargons les données grâce à une copie du fichier enregistrée localement. Un code télécharge directement et enregistre une copie au cas où ce fichier n'exste pas encore en local. ```{r} library(data.table) +library(timeSeries) library(parsedate) + if (!file.exists("monthly_in_situ_co2_mlo.csv")) download.file("https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/monthly/monthly_in_situ_co2_mlo.csv", "monthly_in_situ_co2_mlo.csv", method="auto") data_url <- fread(file = "monthly_in_situ_co2_mlo.csv",skip=54, sep=",",strip.white=TRUE, header = FALSE) @@ -48,7 +50,7 @@ plot(data_clean$DateExcel,data_clean$`CO2[ppm]`) -Bon on va les détecter et les virer y'en a pas non plus trop donc autant faire les choses bien +Bon on va les détecter et les modifier en les supprimant. ```{r} na_records = apply(data_clean, 1, function (x) any(is.na(x)+(x<0))) data_clean[na_records,] @@ -67,8 +69,76 @@ convert_iso_date = function(y,m) { data_clean2$dateStd<-as.Date(convert_iso_date(data_clean2$Yr,data_clean2$Mn)) head(data_clean2) + +plot(data_clean2$dateStd,data_clean2$`CO2[ppm]`) ``` On va pouvoir passer aux stats maintenant je pense, sauf si vous en avez déjà marre mais ça serait dommage d'arrêter en si bon chemin. Cette étude est vraiment pleine de suspens... -## Modèle statistique \ No newline at end of file +## Modèle statistique + +Regardons de plus près l'évolution de la série en focalisant sur des années, par exemple de 1980 à 1985: +```{r} + + +plot(data_clean2[Yr>1980 & Yr<1986,]$dateStd,data_clean2[Yr>1980 & Yr<1986,]$`CO2[ppm]`) +``` + +Nous observons bien une saisonalité de période annuelle. On ne va que s'intéresser dans la suite de cette étude aux données brutes c'est à dire le champ "*CO2(pm)*" car je ne sais pas trop à quoi correpond le reste des champs. Pour isoler la partie tendance, nous allons utiliser une technique simple de moyenne mobile. C'est un estimateur non-paramétrique de la tendance, au sens ou nous ne supposons pas de structure a-priori de cette tendance (par ex. linéaire ou polynomiale). On aurait aussi pu utiliser un lisasge par noyaux mais bon. + +```{r} + +MA<-data_clean2[, .(MA=mean(as.numeric(`CO2[ppm]`))),by=Yr] +MA[,Yr_date:=as.Date(convert_iso_date(Yr,"06"))] +MA$Yr<-as.numeric(MA$Yr) +tail(MA) + +plot(data_clean2$dateStd,data_clean2$`CO2[ppm]`,type='l') +lines(MA$Yr_date,MA$MA,col='red') + + +``` + + + +Faisons une régression sur la tendance. On pourrait enlever les premiers et derniers points car ce ne sont pas des donnés complètes. +```{r} + +plot(MA$Yr,log(MA$MA),col='red',type='l') +reg_lin_log<-lm(MA~log(Yr),data=MA) +summary(reg_lin_log) +plot(reg_lin_log) + +reg_lin<-lm(MA~Yr+I(Yr^2)+I(Yr^3),data=MA) +summary(reg_lin) +plot(reg_lin) +tail(predict(reg_lin)) +plot(MA$Yr,MA$MA,type='l') +lines(MA$Yr,predict(reg_lin),col='red') +``` + +Il est temps de passser à l'extrapolation: +```{r} + +dates_futur<-as.Date(convert_iso_date(seq(2021,2025),"06")) +futur<-data.table(Yr=seq(2021,2025)) + +pred<-predict(reg_lin,futur,interval = "prediction") +plot(MA$Yr,MA$MA,type='l',xlim=c(1960,2025),ylim=c(310,450)) +lines(seq(2021,2025),pred[,1],col='red') + +``` + +Sinon on aurait pu se concentrer sur la distribution des évolutions annuelles +```{r} +evol_ann<-diff(MA$MA)/MA$MA[-length(MA)] +plot(MA$Yr[-1],evol_ann) +plot(MA$Yr[-1],diff(MA$MA)) + + +reg_evol<-lm(evol_ann~MA$Yr[-1]) +summary(reg_evol) +plot(reg_evol) + + +``` \ No newline at end of file