# RMT Modélisation & Analyse de Données - Atelier Série temporelle 10 juillet 2014
# Lucie Michel (ACTA/INRA)

#dossier où ce trouve les fichiers
setwd("/Users/utilisateur/Documents/COURS_series_chrono")

#importation de la bdd blé
TAB_Yield <- read.table("BLE_DPT_surface_prod_rdt_R.txt", header = T)

head(TAB_Yield)

#Vecteurs avec l'ensemble des departements
Dpt_Nom <- levels(TAB_Yield$Dpt_Nom)

#Choix du departement
Dpt_i <- "AIN"
#Dpt_i <- "EURE"
#Dpt_i <- "CREUSE"

#BDDpour le departement voulu
TAB_Yield_i<-TAB_Yield[TAB_Yield$Dpt_Nom == Dpt_i,]

print(TAB_Yield_i)

#Annees 
Year <- TAB_Yield_i$Year
print(Year)

#Rendement
Yield <- TAB_Yield_i$Yield/10
print(Yield)

##############################
##    Modele statistique    ##
##############################

#Modele lineaire
lineaire <- glm(Yield~Year)
print(summary(lineaire))

#Rejutement des annees pour les modeles cubique et quadratique
Year_b <- Year - Year[1]
#Annes pour le modele cubique
Year2_b<-Year_b*Year_b

#Model quadratique
quadratiq <- glm(Yield~Year_b+Year2_b)
print(summary(quadratiq))

#Annees pour le modele cubique
Year3_b <- Year2_b*Year_b

#Model cubique
cubiq <- glm(Yield~Year_b+Year2_b+Year3_b)
print(summary(cubiq))

#Fonction pour le modele lineaire+plateau
LP<- function (t, Ymax, Tmax, P) {
  Y <- Ymax+P*(t-Tmax)
  Y[Y>Ymax] <- Ymax
  Y
}

# Liste pour l'AIN
list1<- list(Ymax= 7, P= 0.11, Tmax= 1998)
# Liste pour l'EURE
#list1<- list(Ymax= 8.5, P= 0.12, Tmax= 1998)
# Liste pour la CREUSE
#list1<- list(Ymax= 5, P= 0.07, Tmax=2000)

print(list1)

#Modele lineaire+plateau
lplateau <- nls(Yield ~LP(Year, Ymax, Tmax, P),start= list1, trace= T)
print(summary( lplateau))
print(AIC(lplateau))

########################################
##    Graphique rendement + modele    ##
########################################

#Graphique pour le modèle linéaire
plot(Year, Yield, xlab="", ylab="Yield (t.ha-1)", type ="l", lwd=1, xlim=c(1950,2011), 
     ylim=c(0,max(Yield)), cex.lab=1.2,col = "blue", col.lab = "blue")
lines(Year, predict(lineaire), lwd=2)
mtext("Year", side =1, line=2.5, at=1980, cex=2)
title(Dpt_i)

#Graphique pour le modèle quadratique
plot(Year, Yield, xlab="", ylab="Yield (t.ha-1)", type ="l", lwd=1, xlim=c(1950,2011), 
     ylim=c(0,max(Yield)), cex.lab=1.2,col = "blue", col.lab = "blue")
lines(Year, predict(quadratiq), lwd=2)
mtext("Year", side =1, line=2.5, at=1980, cex=2)
title(Dpt_i)


#Graphique pour le modèle cubique
plot(Year, Yield, xlab="", ylab="Yield (t.ha-1)", type ="l", lwd=1, xlim=c(1950,2011), 
     ylim=c(0,max(Yield)), cex.lab=1.2,col = "blue", col.lab = "blue")
lines(Year, predict(cubiq), lwd=2)
mtext("Year", side =1, line=2.5, at=1980, cex=2)
title(Dpt_i)

#Graphique pour le modèle linéaire + plateau
plot(Year, Yield, xlab="", ylab="Yield (t.ha-1)", type ="l", lwd=1, xlim=c(1950,2011), 
     ylim=c(0,max(Yield)), cex.lab=1.2,col = "blue", col.lab = "blue")
lines(Year, predict(lplateau), lwd=2)
mtext("Year", side =1, line=2.5, at=1980, cex=2)
title(Dpt_i)

########################
##    Calcul de R2    ##
########################

#R2 pour le modele lineaire
R_lineaire <-1-((sum((Yield-predict(lineaire))^2))/(sum((Yield-mean(Yield))^2)))
print(R_lineaire)

#R2 pour le modele quadratique
R_quadratiq <-1-((sum((Yield-predict(quadratiq))^2))/(sum((Yield-mean(Yield))^2)))
print(R_quadratiq)

#R2 pour le modele cubique
R_cubiq <-1-((sum((Yield-predict(cubiq))^2))/(sum((Yield-mean(Yield))^2)))
print(R_cubiq)

#R2 pour le modele lineaire + plateau
R_lplateau <-1-((sum((Yield-predict(lplateau))^2))/(sum((Yield-mean(Yield))^2)))
print(R_lplateau)

#################################
##    Graphique des résidus    ##
#################################

#Graphique des résidus pour le modele lineaire 
plot(Year,Yield-predict(lineaire))
abline(0,0)
title(Dpt_i)

#Graphique des résidus pour le modele quadratique
plot(Year,Yield-predict(quadratiq))
abline(0,0)
title(Dpt_i)

#Graphique des résidus pour le modele cubique
plot(Year,Yield-predict(cubiq))
abline(0,0)
title(Dpt_i)

#Graphique des résidus pour le modele  lineaire + plateau
plot(Year,Yield-predict(lplateau))
abline(0,0)
title(Dpt_i)