#setwd("C:/Users/pgarnier/Documents/E/pesticide/evaluation-perso/pessac/Myriam")   

# ----------------------------------
# Chargement des fonctions et packages
# ----------------------------------

library(nlme)
library(ggplot2)
library(scatterplot3d)

# ----------------------------------
# Chargement jdd
# ----------------------------------

don=read.table("respi_vdt_PG_03_2019.csv", header=TRUE,sep=";", dec=".",skip=1)

# ----------------------------------
# Graphique donnees brutes
# ----------------------------------


summary(don)

# programme Michel pour box plot il faut l adapter
# rm( list=ls() )
# PourIFTt=read.table('PourIFTt.txt',header=TRUE,na.string="NA")
# PourIFTt
# PourIFTt$Trait<-factor(PourIFTt$Trait,levels=c("P","I", "SCV"))
# PourIFTt$Trait
# boxplot(IFTt~PourIFTt$Trait, data=PourIFTt) 

Pl=ggplot()+
  geom_line(data = don, aes(x = temps, y = respiration, color = traitement), size=0.5) +
  geom_point(data =don, aes(x = temps, y = respiration,color = traitement), size = 0.5) +
  scale_colour_manual(values=c("red","blue"))+
  scale_fill_manual(values=c("1","2"))+
  facet_wrap(~expe,scales="free_y")+
  theme(axis.text=element_text(size=8),
        axis.title=element_text(size=14,face="bold"))
Pl



# ----------------------------------
# Mise en forme jeu de donnees
# ----------------------------------

Exp<-c()
Time<-c()
Emission<-c()
Esp<-c()
rep<-c()
input=c()
type=c()
C_N=c()
temp=c()
densite=c()
Land_occupation=c()
cat_mo=c()

for (i in 1:length(unique(don$expe))) {
  
  Exp<-c(Exp,don$expe[don$expe==unique(don$expe)[i]&don$traitement=="CTRL"])
  Time<-c(Time,don$temps[don$expe==unique(don$expe)[i]&don$traitement=="CTRL"])
  Emission<-c(Emission,don$respiration[don$expe==unique(don$expe)[i]&don$traitement=="EW"]/don$respiration[don$expe==unique(don$expe)[i]&don$traitement=="CTRL"])
  Esp<-c(Esp,don$categorie[don$expe==unique(don$expe)[i]&don$traitement=="CTRL"])
  rep<-c(rep,don$Rep[don$expe==unique(don$expe)[i]&don$traitement=="CTRL"])
  input<-c(input,don$Mat_ajoutee[don$expe==unique(don$expe)[i]&don$traitement=="CTRL"])
  type<-c(type,don$type_exp[don$expe==unique(don$expe)[i]&don$traitement=="CTRL"])
  C_N<-c(C_N,don$C_N_Mat[don$expe==unique(don$expe)[i]&don$traitement=="CTRL"])
  temp=c(temp,don$temp[don$expe==unique(don$expe)[i]&don$traitement=="CTRL"])
  densite=c(densite,don$densite[don$expe==unique(don$expe)[i]&don$traitement=="CTRL"])
  Land_occupation<-c(Land_occupation,don$Land_occupation[don$expe==unique(don$expe)[i]&don$traitement=="CTRL"])
  cat_mo<-c(cat_mo,don$cat_MO[don$expe==unique(don$expe)[i]&don$traitement=="CTRL"])
}

Time2<-Time*Time
Time3<-Time*Time*Time
w<-1/rep

cat_mo<-as.factor(cat_mo)

TAB<-data.frame(Exp,Esp,Time,Time2,Time3,Emission,w,input,type,C_N,temp,densite,Land_occupation,cat_mo)
TAB<-TAB[is.na(TAB$Emission)==F,]
TAB$Emission_log=log(TAB$Emission)

#TAB2<-TAB[TAB$Esp==1|TAB$Esp==2,]

#*****************************************
#
# Il faut arranger l'histoire de groupe ci dessous parce que sur la figure les petits graphes sont bien par numéros
# de manip croissant mais ça ne colle plus avec le vrai numéro (la manip indiquée x n'est pas la xième); lignes 94 et 109

#*****************************************

Data<-groupedData(Emission_log~Time|Exp, data=TAB)# a commenter et faire Data <- TAB 

#Data$input=as.factor(Data$input)
#levels(Data$input)=levels(don$Mat_ajoutee)

#Data$Esp=as.factor(Data$Esp)
#levels(Data$Esp)=levels(don$categorie)

# ----------------------------------
# Graphique log ratio
# ----------------------------------

Pl2=ggplot()+
  geom_line(data = TAB, aes(x = Time, y = Emission_log)) +
  geom_point( data =TAB, aes(x = Time, y = Emission_log)) + ylab("Emission log Ratio")+
  xlab("Days")+
  facet_wrap(~as.numeric(Exp),scales="free",ncol=6)+#on reforce Exp a redevenir un nombre et pas un facteur
  theme(axis.title.x=element_text(face = "bold",size=20),
        axis.title.y=element_text(face = "bold",size=20))
#+ ggsave(Pl2,filename = "joligraf.png",width=9.5,height=9)
Pl2


# ----------------------------------
# Comparaison modeles AIC/BIC
# ----------------------------------

#Data<-TAB

#modele mixte
Mod0<-lme(Emission_log~1,random=pdDiag(~1),data=Data,method="ML");summary(Mod0)[21:22]  #modele cst
Mod1<-lme(Emission_log~Time,random=pdDiag(~1+Time),data=Data,method="ML");summary(Mod1)[21:22]    #modele lineaire
Mod2<-lme(Emission_log~Time+Time2,random=pdDiag(~1+Time+Time2),data=Data,method="ML");summary(Mod2)[21:22]    #modele quadratique
Mod3<-lme(Emission_log~Time+Time2+Time3,random=pdDiag(~1+Time+Time2+Time3),data=Data,method="ML");summary(Mod3)[21:22]   #modele cubique

#modele effet fixe
Mod_fixe0=glm(log(Emission)~1,data=Data);AIC(Mod_fixe0);BIC(Mod_fixe0)
Mod_fixe1=glm(log(Emission)~Time,data=Data);AIC(Mod_fixe1);BIC(Mod_fixe1)
Mod_fixe2=glm(log(Emission)~Time+Time2,data=Data);AIC(Mod_fixe2);BIC(Mod_fixe2)
Mod_fixe3=glm(log(Emission)~Time+Time2+Time3,data=Data);AIC(Mod_fixe3);BIC(Mod_fixe3)
Mod_fixe3_dens=glm(log(Emission)~Time+Time2+Time3+densite,data=Data);AIC(Mod_fixe3_dens);BIC(Mod_fixe3_dens)
Mod_fixe3_input=glm(log(Emission)~Time+Time2+Time3+cat_mo,data=Data);AIC(Mod_fixe3_input);BIC(Mod_fixe3_input)
Mod_fixe3_esp=glm(log(Emission)~Time+Time2+Time3+Esp,data=Data);AIC(Mod_fixe3_esp);BIC(Mod_fixe3_esp)
Mod_fixe3_temp=glm(log(Emission)~Time+Time2+Time3+temp,data=Data);AIC(Mod_fixe3_temp);BIC(Mod_fixe3_temp)
Mod_fixe3_lo=glm(log(Emission)~Time+Time2+Time3+Land_occupation,data=Data);AIC(Mod_fixe3_lo);BIC(Mod_fixe3_lo)
Mod_fixe3_int=glm(log(Emission)~Time*densite+Time2+Time3,data=Data);AIC(Mod_fixe3_int);BIC(Mod_fixe3_int)
summary(Mod_fixe3_int)
# ----------------------------------
# Ajout covariables au modele mixte cubique
# ----------------------------------

#modele cubique
Mod3<-lme(Emission_log~Time+Time2+Time3,random=pdDiag(~1+Time+Time2+Time3),data=Data,method="ML");summary(Mod3)[21:22]

#modele cubique + densite
Mod_dens<-lme(Emission_log~Time+Time2+Time3+densite,random=pdDiag(~1+Time+Time2+Time3+densite),data=Data,method="ML");summary(Mod_dens)[21:22]
summary(Mod_dens)

#modele cubique + esp
Mod_esp<-lme(Emission_log~Time+Time2+Time3+Esp,random=pdDiag(~1+Time+Time2+Time3+Esp),data=Data,method="ML");summary(Mod_esp)[21:22]

#modele cubique + mo
Mod_mo<-lme(Emission_log~Time+Time2+Time3+cat_mo,random=pdDiag(~1+Time+Time2+Time3+cat_mo),data=Data,method="ML");summary(Mod_mo)[21:22]
summary(Mod_mo)

#modele cubique + lo
Mod_lo<-lme(Emission_log~Time+Time2+Time3+Land_occupation,random=pdDiag(~1+Time+Time2+Time3+Land_occupation),data=Data,method="ML");summary(Mod_lo)[21:22]


#modele cubique + temperature

Data_temp<-Data[is.na(Data$temp)==F,]
Mod_temp_0<-lme(Emission_log~Time+Time2+Time3,random=pdDiag(~1+Time+Time2+Time3),data=Data_temp,method="ML");summary(Mod_temp_0)[21:22]
Mod_temp<-lme(Emission_log~Time+Time2+Time3+temp,random=pdDiag(~1+Time+Time2+Time3+temp),data=Data_temp,method="ML");summary(Mod_temp)[21:22]
summary(Mod_temp)
#modele cubique + interraction densite
Mod_int<-lme(Emission_log~Time*densite+Time2+Time3,random=pdDiag(~1+Time*densite+Time2+Time3),data=Data,method="ML");summary(Mod_int)[21:22]



# ----------------------------------
# Ajustement modele final cubique + densite mixte
# ---------------------------------
#le model mixte avec densite s apelle mod et il n'est pas fixed il faut prendre REML a partir de ce niveau
Mod<-lme(Emission_log~Time+Time2+Time3+densite,random=pdDiag(~1+Time+Time2+Time3),data=Data,method="REML")
Mod_int<-lme(Emission_log~Time*densite+Time2+Time3,random=pdDiag(~1+Time*densite+Time2+Time3),data=Data,method="REML")

#Mod<-lme(Emission_log~Time+densite,random=pdDiag(~1+Time+Time2+Time3),data=Data,method="REML")
summary(Mod)
summary(Mod_int)
hist(residuals(Mod_int))
library(moments)
kurtosis(residuals(Mod_int))
skewness(residuals(Mod_int))

#Three responses corresponding to contrasted densities
par(mfrow=c(1,2))

TimeVec<-1:380
DensVec<-c(0.5,1.95,9,28)
Para_Mod_int<-fixef(Mod_int)
Pred1<-Para_Mod_int[1]+Para_Mod_int[2]*TimeVec+Para_Mod_int[3]*DensVec[1]+Para_Mod_int[4]*TimeVec^2+Para_Mod_int[5]*TimeVec^3+Para_Mod_int[6]*TimeVec*DensVec[1]
Pred2<-Para_Mod_int[1]+Para_Mod_int[2]*TimeVec+Para_Mod_int[3]*DensVec[2]+Para_Mod_int[4]*TimeVec^2+Para_Mod_int[5]*TimeVec^3+Para_Mod_int[6]*TimeVec*DensVec[2]
Pred3<-Para_Mod_int[1]+Para_Mod_int[2]*TimeVec+Para_Mod_int[3]*DensVec[3]+Para_Mod_int[4]*TimeVec^2+Para_Mod_int[5]*TimeVec^3+Para_Mod_int[6]*TimeVec*DensVec[3]
Pred4<-Para_Mod_int[1]+Para_Mod_int[2]*TimeVec+Para_Mod_int[3]*DensVec[4]+Para_Mod_int[4]*TimeVec^2+Para_Mod_int[5]*TimeVec^3+Para_Mod_int[6]*TimeVec*DensVec[4]
plot(TimeVec,Pred1, xlab="Time", ylab="log ratio",type="l", xlim=c(0,300), ylim=c(0,1.5))
lines(TimeVec, Pred2, lwd=2)
lines(TimeVec, Pred3, lwd=3)
lines(TimeVec, Pred4, lwd=4)

#Three responses corresponding to contrasted densities
TimeVec<-1:380
DensVec<-c(0.5,1.95,9, 28)
Para_Mod_int<-fixef(Mod)
Pred1<-Para_Mod_int[1]+Para_Mod_int[2]*TimeVec+Para_Mod_int[5]*DensVec[1]+Para_Mod_int[3]*TimeVec^2+Para_Mod_int[4]*TimeVec^3
Pred2<-Para_Mod_int[1]+Para_Mod_int[2]*TimeVec+Para_Mod_int[5]*DensVec[2]+Para_Mod_int[3]*TimeVec^2+Para_Mod_int[4]*TimeVec^3
Pred3<-Para_Mod_int[1]+Para_Mod_int[2]*TimeVec+Para_Mod_int[5]*DensVec[3]+Para_Mod_int[3]*TimeVec^2+Para_Mod_int[4]*TimeVec^3
Pred4<-Para_Mod_int[1]+Para_Mod_int[2]*TimeVec+Para_Mod_int[5]*DensVec[4]+Para_Mod_int[3]*TimeVec^2+Para_Mod_int[4]*TimeVec^3
plot(TimeVec,Pred1, xlab="Time", ylab="log ratio",type="l", xlim=c(0,300), ylim=c(0,1.5))
lines(TimeVec, Pred2, lwd=2)
lines(TimeVec, Pred3, lwd=3)
lines(TimeVec, Pred4, lwd=4)


pvalue_coef=summary(Mod)$tTable;pvalue_coef  #coefficients estimes dans colonne "value" et pvalue colonne "p-value"

# ----------------------------------
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# Pat a ecrit 182 a 189 Ajustement modele final cubique + densite  modele fixe
# ---------------------------------
#le model fixe
#fixemod <- lm(Emission_log~Time+Time2+Time3+densite+Time:densite,data=Data)
fixemod <- lm(Emission_log~Time+Time2+Time3+densite,data=Data)
#fixemod <- lm(Emission_log~Time+densite,data=Data)
summary(fixemod)
summary(fixemod)$coefficients[,c("Estimate","Pr(>|t|)")]


# ----------------------------------
# Graphiques modele avec donnees
# ---------------------------------

#donnees avec prediction pour chaque exp
res=as.data.frame(Data$Time)
names(res)="Time"
res$Exp=Data$Exp
res$ratio=Data$Emission_log
res$type=as.factor(rep("observed",dim(res)[1]))


res4=as.data.frame(Data$Time)
names(res4)="Time"
res4$Exp=Data$Exp
res4$ratio=predict(Mod_int,Data)
res4$type=as.factor(rep("Fitted cubic model inter time density (random effects)",dim(res4)[1]))
res5=rbind(res4,res)

res6=as.data.frame(Data$Time)
names(res6)="Time"
res6$Exp=Data$Exp
res6$ratio=predict(Mod_fixe3_int,Data)
res6$type=as.factor(rep("Fitted cubic model inter time density (fixed effects)",dim(res4)[1]))
res7=rbind(res5,res6)

#Ajout d'une variable pour avoir les graphiques par ordre croissant de numéro de Exp
NameExp<-as.numeric(as.character(res7$Exp))
res7<-cbind(res7,NameExp)

res7$type=factor(res7$type,c("observed","Fitted cubic model inter time density (random effects)","Fitted cubic model inter time density (fixed effects)"))

res7

#  geom_point()+

p1<- ggplot(res7,aes(x=Time,y=ratio,linetype=type,colour=type)) + 
  geom_line()+
  scale_linetype_manual(values=c("blank","longdash","dashed"))+
  scale_color_manual(values=c('blue','green', 'red'))+
  geom_point(aes(size=type))+
  scale_shape_manual(values=c(16, 0, 0))+
  scale_size_manual(values=c(1,0,0))+
  ylab("Emission log Ratio")+
  xlab("Days")+
  facet_wrap(~NameExp,scales="free")+
  theme(legend.justification = "center",
        legend.position = "bottom", 
        legend.direction="vertical", 
        legend.background = element_blank(),
        panel.grid.major = element_line(colour = "grey"),
        panel.grid.minor = element_blank(),
        panel.background = element_rect(fill = "white"),
        axis.title.x=element_text(face = "bold",size=15),
        axis.title.y=element_text(face = "bold",size=15),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank(),
        axis.text=element_text(color = "black", size = 20),
        panel.border = element_rect(fill = NA, colour = "grey50"),
        axis.text.y=element_text(face = "bold",size=5),
        axis.text.x=element_text(face = "bold",size=5),
        legend.text=element_text(color = "black", size = 10),
        plot.title = element_text(vjust=1.5,face = "bold",size=20))
p1

#predit vs observe
Data$pred=predict(Mod,Data)
plot(Data$pred~Data$Emission_log,xlab="observed",ylab="fitted cubic model",xlim=c(-0.5,1.5),ylim=c(-0.5,2))
abline(a=0,b=1)

#residus en fonction du temps
plot(residuals(Mod)~Data$Time,xlab="Days",ylab="Residuals cubic model",ylim=c(-0.3,0.3))
abline(a=0,b=0)

#coefficients par experience
plot(coef(Mod),xlab="coefficients fitted cubic model + density",nrow=1)

#coefficients par experience tries
plot(coef(Mod)[order(coefficients(Mod)[,1]),],xlab="Coefficients model cubic + density",cex=1)

# ----------------------------------
# Graphique avec IC pour faible et forte densite
# ---------------------------------

coef_fixe=Mod$coefficients$fixed

densite=150  #rentrer valeur densite

plot(coef_fixe[1]+coef_fixe[2]*(0:400)+coef_fixe[3]*(0:400)^2+
       coef_fixe[4]*(0:400)^3+coef_fixe[5]*densite,
     ylim=c(-6,6),ylab="predicted values",xlab="days")



lines(coef_fixe[1]+coef_fixe[2]*(0:400)+coef_fixe[3]*(0:400)^2+
        coef_fixe[4]*(0:400)^3+coef_fixe[5]*densite+
       0.2643036^2+((0.005762125)^2)*(0:400)^2+((1.380068e-06)^2)*(0:400)^4+((4.536256e-13)^2)*(0:400)^6,
      col="red")

lines(coef_fixe[1]+coef_fixe[2]*(0:400)+coef_fixe[3]*(0:400)^2+
        coef_fixe[4]*(0:400)^3+coef_fixe[5]*densite-
        ((0.2643036^2+((0.005762125)^2)*(0:400)^2+((1.380068e-06)^2)*(0:400)^4+((4.536256e-13)^2)*(0:400)^6)),
      col="red")





# ----------------------------------
# Calcul des 2 pentes
# ----------------------------------


Model<-function(X, A,B,D,S)
{ 
  Y<-rep(NA,length(X))
  C<-A+S*(B-D)
  Y[X<=S]<-A[X<=S]+B[X<=S]*X[X<=S]
  Y[X>S]<-C[X>S]+D[X>S]*X[X>S]
  print(Y)
  return(Y)
}


summary(Data)


library(nlme)
Fit<-nlme(Emission_log~Model(Time,A,B,D,S),start=c(A=0.03,B=0.01,S=50,D=0),random=pdDiag(A+B+S+D~1),fixed=list(A+B+S+D~1),
          data=Data)
summary(Fit)

plot(coefficients(Fit))

plot(Fit)


Model<-function(X, B,D,S)
{ 
  Y<-rep(NA,length(X))
  C<-S*(B-D)
  Y[X<=S]<-B[X<=S]*X[X<=S]
  Y[X>S]<-C[X>S]+D[X>S]*X[X>S]
  print(Y)
  return(Y)
}


summary(Data)


library(nlme)
Fit<-nlme(Emission_log~Model(Time,B,D,S),start=c(B=0.01,S=50,D=0),random=pdDiag(B+S+D~1),fixed=list(B+S+D~1),
          data=Data)
summary(Fit)

plot(coefficients(Fit))

plot(Fit)

# ----------------------------------
# graphique 3D temps + densite
# ----------------------------------

Forte=Data[which(Data[,"densite"]>20),]
scatterplot3d(Forte$Time,Forte$densite,log(Forte$Emission),angle=-150,
              color=as.numeric(as.factor(as.character(Forte$Exp))),zlim=c(0,2))

Faible=Data[which(Data[,"densite"]<20),]
scatterplot3d(Faible$Time,Faible$densite,log(Faible$Emission),angle=-150,
              color=as.numeric(as.factor(as.character(Faible$Exp))),zlim=c(0,2))