#data of Ashauer, 2012
#2,4,5-trichlorophenol
#choose directory session

#import data
data888 <- read.table("ashauer_trichlorophenol.csv", header=TRUE,
                      sep=";")
t<-c(data888[1:19,3])
Cobs<-c(data888[1:19,9])
Cw<-mean(data888[1:7,6])
C0<-Cobs[[1]]  
Cobsmet<-data888[1:19,11]
Cobsmet2<-data888[1:19,12]
plot(Cobs~t)
tc=0.979

require(rjags)
#model
mod1<- 
  "model
{
  
  for(i in 1:7) #accu 
  {
  Cpredald[i] <- ((kw*Cw)/(ke+km1+km2)) + (C0-((kw*Cw)/(ke+km1+km2))* (exp(-(ke+km1+km2)*t[i])))
  Cobsald[i] ~ dnorm(Cpredald[i],tau)   
  
  CpredMet1[i] <-  (((km1*(kw*Cw)/(ke+km1+km2))/((ke+km1+km2)-kemet1))*(exp(-(ke+km1+km2)*t[i])-exp(-kemet1*t[i]))) +
  ((((km1*(kw*Cw)/(ke+km1+km2))/kemet1)*(1-exp(-kemet1*t[i]))))
  CobsMet1[i] ~ dnorm(CpredMet1[i],tauMet1)  
  
  CpredMet2[i] <-  (((km2*(kw*Cw)/(ke+km1+km2))/((ke+km1+km2)-kemet2))*(exp(-(ke+km1+km2)*t[i])-exp(-kemet2*t[i]))) +
  ((((km2*(kw*Cw)/(ke+km1+km2))/kemet2)*(1-exp(-kemet2*t[i]))))
  CobsMet2[i] ~ dnorm(CpredMet2[i],tauMet2)  
  }
  
  for(i in 8:19) #depu 
  {
  Cpredald[i] <- ((kw*Cw)/(ke+km1+km2) * exp(-(ke+km1+km2)*(t[i]-tc)))+((C0-((kw*Cw)/(ke+km1+km2))*exp(-(ke+km1+km2)*t[i]))) 
  Cobsald[i] ~ dnorm(Cpredald[i],tau)  
  
  CpredMet1[i] <-  ((((km1*(kw*Cw)/(ke+km1+km2))/((ke+km1+km2)-kemet1))*(exp(-(ke+km1+km2)*t[i])-exp(-kemet1*t[i]))) +
  ((((km1*(kw*Cw)/(ke+km1+km2))/kemet1)*(1-exp(-kemet1*t[i])))))*exp(-kemet1*(t[i]-tc))
  CobsMet1[i] ~ dnorm(CpredMet1[i],tauMet1)  
  
  CpredMet2[i] <-  ((((km2*(kw*Cw)/(ke+km1+km2))/((ke+km1+km2)-kemet2))*(exp(-(ke+km1+km2)*t[i])-exp(-kemet2*t[i]))) +
  ((((km2*(kw*Cw)/(ke+km1+km2))/kemet2)*(1-exp(-kemet2*t[i])))))*exp(-kemet2*(t[i]-tc))
  CobsMet2[i] ~ dnorm(CpredMet2[i],tauMet2)   
  }
  
  #priors
  logkw~dunif(-10,5)                #dnorm(3.391641,0.1903153) #prior Ashauer 2010
  kw<-10^logkw
  logke ~ dunif(-10,2)              #dnorm(-0.02918839,1.124318) #prior Ashauer 2010
  ke<-10^logke
  logkm1 ~ dunif(-5,2)    
  km1<-10^logkm1 
  logkm2 ~ dunif(-5,2)    
  km2<-10^logkm2 
  logkemet1 ~ dunif(-5,2)     
  kemet1<-10^logkemet1 
  logkemet2 ~ dunif(-5,2)     
  kemet2<-10^logkemet2 
  sigma ~ dgamma(0.001,0.001)       
  tau <- 1/(sigma*sigma)    
  tauMet1<- 1/(sigmaMet1*sigmaMet1) 
  sigmaMet1 ~ dgamma(0.001,0.001) 
  tauMet2<- 1/(sigmaMet2*sigmaMet2) 
  sigmaMet2 ~ dgamma(0.001,0.001)
}
"

#definition of data for the model
data=list(Cobsald=Cobs,CobsMet1=Cobsmet,CobsMet2=Cobsmet2,Cw=Cw,t=t,C0=C0, tc=tc)

inits=list(list(logkm1=0,logkm2=0,logkemet1=0,logkemet2=0,logkw=2,logke=0.372912,sigma=10,sigmaMet1=10,sigmaMet2=10),
           list(logkm1=1,logkm2=1,logkemet1=1,logkemet2=0,logkw=1,logke=0.372912,sigma=10,sigmaMet1=10,sigmaMet2=10),
           list(logkm1=0.5,logkm2=0.5,logkemet1=0.5,logkemet2=0,logkw=-1,logke=0.372912,sigma=10,sigmaMet1=10,sigmaMet2=10))
#model implementation
m1<- jags.model(file = textConnection(mod1), inits=inits,data=data, n.chains = 3)
update(m1,20000)
#dic<-dic.samples(m1, n.iter=301000, thin=85)
mcmc1<-coda.samples(m1, c("km1","km2","kemet1","kemet2","kw","ke","sigma","sigmaMet1","sigmaMet2"), n.iter= 2000000, thin=1000) #2000000 and 1000

#Rafetry and Lewis (1992) method
#mcmc1<-coda.samples(m1, c("km1","km2","kemet1","kemet2","kw","ke","sigma","sigmaMet1","sigmaMet2"), n.iter=5000, thin=1) 
#RL <- raftery.diag(mcmc1)
#resmatrixtot <- rbind(RL[[1]]$resmatrix,RL[[2]]$resmatrix,RL[[3]]$resmatrix )
#thin <- round(max(resmatrixtot[,"I"])+0.5)
#niter <- max(resmatrixtot[,"Nmin"])*thin
#mcmc1<-coda.samples(m1, c("km1","km2","kemet1","kemet2","kw","ke","sigma","sigmaMet1","sigmaMet2"), n.iter=niter, thin=thin)

#####################################################################
#Results
summary(mcmc1) #parameters  2.5, 50 and 97.5% quantiles
gelman.diag(mcmc1) #must be close to 1 
geweke.diag(mcmc1) #must be between -2 and 2
autocorr.plot(mcmc1[[1]]) #Visualization for autocorrelation
plot(mcmc1, density=F) #MCMC traces
mcmctot1<-as.data.frame(as.matrix(mcmc1)) #total of the 3 MCMC chains
mcmctotsample1<-mcmctot1[sample.int(nrow(mcmctot1), size=500),] #sample of 500 simulations

#####################################################################
#Prior and posterior distributions
#prior
data00=list(t=t,Cw=Cw,C0=C0,tc=tc) 
m00<- jags.model(file = textConnection(mod1), inits=inits,data=data00, n.chains = 3)
update(m00,50000)
mcmc00<-coda.samples(m00, c("kw","ke","sigma","kemet1","kemet2","km1","km2","kemet2","sigmaMet1","sigmaMet2"), n.iter=100000, thin=10)
mcmctot00<-as.data.frame(as.matrix(mcmc00)) 
#plot prior against posterior distributions
par(mfrow=c(2,2))
for (i in 1:ncol(mcmctot1))
{
  plot(density(mcmctot1[,i]), lwd=4,main=names(mcmctot1)[i])
  lines(density(mcmctot00[,i]), col="red", lwd=2)
}

######################################################################
#Predictions
vkw<-mcmctotsample1[,"kw"]
vke<-mcmctotsample1[,"ke"]
vkemet1<-mcmctotsample1[,"kemet1"]
vkemet2<-mcmctotsample1[,"kemet2"]
vkm1<-mcmctotsample1[,"km1"]
vkm2<-mcmctotsample1[,"km2"]
vsigma<-mcmctotsample1[,"sigma"]
vsigmaMet1<-mcmctotsample1[,"sigmaMet1"]
vsigmaMet2<-mcmctotsample1[,"sigmaMet2"]
vCpredald<-mcmctotsample1[,"Cpredald"]
vCpredMet1<-mcmctotsample1[,"CpredMet1"]
vCpredMet2<-mcmctotsample1[,"CpredMet2"]

l<-500
vt<-seq(0,6.125,length.out = l)
niter= 2000000

qinf<-vector(length = l); qmed<-vector(length = l); qsup<-vector(length = l)
qinfm<-vector(length = l); qmedm<-vector(length = l); qsupm<-vector(length = l)
qinfm2<-vector(length = l); qmedm2<-vector(length = l); qsupm2<-vector(length = l)
niter=2000000

for (i in 1:l)
{ 
  if(vt[i]<1.00001)
  {vCpredE1<-(((vkw*Cw))/(vke+vkm1+vkm2)) + (C0-((vkw*Cw)/(vke+vkm1+vkm2))* (exp(-(vke+vkm1+vkm2)*vt[i])))
  vCpredE1obs <- rnorm(niter, vCpredE1, vsigma)}
  else {vCpredE1<-(((vkw*Cw)/(vke+vkm1+vkm2)) * exp(-(vke+vkm1+vkm2)*(vt[i]-tc)))+((C0-((vkw*Cw)/(vke+vkm1+vkm2)))*exp(-(vke+vkm1+vkm2)*vt[i])) 
  vCpredE1obs <- rnorm(niter, vCpredE1, vsigma)}
  
  if(vt[i]<1.00001)
  {vCpredMet1 <-  (((vkm1*(vkw*Cw)/(vke+vkm1+vkm2))/((vke+vkm1+vkm2)-vkemet1))*(exp(-(vke+vkm1+vkm2)*vt[i])-exp(-vkemet1*vt[i]))) +
    (((vkm1*(vkw*Cw)/(vke+vkm1+vkm2))/vkemet1)*(1-exp(-vkemet1*vt[i])))
  vCpredMet1obs <- rnorm(niter, vCpredMet1, vsigmaMet1)}
  else{vCpredMet1<- ((((vkm1*(vkw*Cw)/(vke+vkm1+vkm2))/((vke+vkm1+vkm2)-vkemet1))*(exp(-(vke+vkm1+vkm2)*(vt[i]))-exp(-vkemet1*(vt[i])))) +
                       (((vkm1*(vkw*Cw)/(vke+vkm1+vkm2))/vkemet1)*(1-exp(-vkemet1*(vt[i])))))*exp(-vkemet1*(vt[i]-tc))
  vCpredMet1obs <- rnorm(niter, vCpredMet1, vsigmaMet1)}
  
  if(vt[i]<1.00001)
  {vCpredMet2 <-  (((vkm2*(vkw*Cw)/(vke+vkm1+vkm2))/((vke+vkm1+vkm2)-vkemet2))*(exp(-(vke+vkm1+vkm2)*vt[i])-exp(-vkemet2*vt[i]))) +
    (((vkm2*(vkw*Cw)/(vke+vkm1+vkm2))/vkemet2)*(1-exp(-vkemet2*vt[i])))
  vCpredMet2obs <- rnorm(niter, vCpredMet2, vsigmaMet2)
  }
  else{vCpredMet2<- ((((vkm2*(vkw*Cw)/(vke+vkm1+vkm2))/((vke+vkm1+vkm2)-vkemet2))*(exp(-(vke+vkm1+vkm2)*(vt[i]))-exp(-vkemet2*(vt[i])))) +
                       (((vkm2*(vkw*Cw)/(vke+vkm1+vkm2))/vkemet2)*(1-exp(-vkemet2*(vt[i])))))*exp(-vkemet2*(vt[i]-tc))
  vCpredMet2obs <- rnorm(niter, vCpredMet2, vsigmaMet2)
  }
  
  
  qinf[i]<-quantile(vCpredE1obs, probs=0.025)
  qmed[i]<-quantile(vCpredE1obs, probs=0.5)
  qsup[i]<-quantile(vCpredE1obs, probs=0.975)
  
  qinfm[i]<-quantile(vCpredMet1obs, probs=0.025)
  qmedm[i]<-quantile(vCpredMet1obs, probs=0.5)
  qsupm[i]<-quantile(vCpredMet1obs, probs=0.975)
  
  qinfm2[i]<-quantile(vCpredMet2obs, probs=0.025)
  qmedm2[i]<-quantile(vCpredMet2obs, probs=0.5)
  qsupm2[i]<-quantile(vCpredMet2obs, probs=0.975)
}

qinf[which(qinf < 0)] <- 0
qmed[which(qinf < 0)] <- 0
qsup[which(qinf < 0)] <- 0
qinfm[which(qinfm < 0)] <- 0
qmedm[which(qinfm < 0)] <- 0
qsupm[which(qinfm < 0)] <- 0
qinfm2[which(qinfm2 < 0)] <- 0
qmedm2[which(qinfm2 < 0)] <- 0
qsupm2[which(qinfm2 < 0)] <- 0

#Plot Predictions
library(Hmisc)
par(mfrow=c(1,1))
plot(Cobs~t, pch=16, ylim=c(0,17000),col="black",xlab="Time (d)", ylab="Contaminant concentration in organism (nmol g-1 ww)") #donnes observees
polygon(c(vt,rev(vt)),c(qinfSm2,rev(qsupSm2)),col="#fde4c4",border=NA, lwd=2)
polygon(c(vt,rev(vt)),c(qinfSm,rev(qsupSm)),col="#fde4f3",border=NA, lwd=2)
polygon(c(vt,rev(vt)),c(qinfS,rev(qsupS)),col="#edc3c3",border=NA, lwd=2)
polygon(c(vt,rev(vt)),c(qinf,rev(qsup)),col="#c1d2d2fe",border=NA, lwd=2)
polygon(c(vt,rev(vt)),c(qinfm,rev(qsupm)),col="#99ca968a",border=NA, lwd=2)
polygon(c(vt,rev(vt)),c(qinfm2,rev(qsupm2)),col="#99ca968a",border=NA, lwd=2)

points(Cobs~t, pch=16, col="darkgreen")
points(Cobsmet~t, pch=16, col="chartreuse3")
points(Cobsmet2~t, pch=16, col="#089043")
lines(vt, qmed, col="darkgreen", lwd=2, lty=1)
lines(vt, qmedm, col="chartreuse3", lwd=2, lty=1)
lines(vt, qmedm2, col="#089043", lwd=2, lty=1)
abline(v=tc,lty=2)
lines(vt, qmedS, col="#b40e3e", lwd=2, lty=2)
lines(vt, qmedSm, col="#e09c9c", lwd=2, lty=2)
lines(vt, qmedSm2, col="orange", lwd=2, lty=2)


#####################################################################
#Comparison of parameter: bayesian and orignal inference methods

#50% mean
kw=mean(vkw)
kwS=1389
ke=mean(vke)
keS=0
km1=mean(vkm1)
kmS=14.52
kemet1=mean(vkemet1)
kemetS=0.822
km2=mean(vkm2)
kmS2=2.35
kemet2=qmean(vkemet2)
kemetS2=1.78

#2.5%
qkw=quantile(vkw, probs=0.025)
qkwS=0
qke=quantile(vke, probs=0.025)
qkeS=0
qkm=quantile(vkm1, probs=0.025)
qkmS=1.1
qkemet=quantile(vkemet1, probs=0.025)
qkemetS=0.42
qkm2=quantile(vkm2, probs=0.025)
qkmS2=0
qkemet2=quantile(vkemet2, probs=0.025)
qkemetS2=0

#97.5%
skw=quantile(vkw, probs=0.975)
skwS=4073
ske= quantile(vke, probs=0.975)
skeS=36
skm= quantile(vkm1, probs=0.975)
skmS=28.0
skemet=quantile(vkemet1, probs=0.975)
skemetS=1.23
skm2= quantile(vkm2, probs=0.975)
skmS2=9.7
skemet2=quantile(vkemet2, probs=0.975)
skemetS2=9

library(Hmisc)
par(mfrow=c(1,1))
plot(c(kw,kwS,ke,keS,km,kmS,kemet,kemetS,km2,kmS2,kemet2,kemetS2)~c(0:11),xlim=c(0,11), ylim=c(0,100), xaxt="n", pch=16,
     ylab="model parameter value",xlab="Model parameter")
axis(1, at =c(0:11), labels = c("kw Bayesian", "kw Ashauer","ke Bayesian", "ke Ashauer","km Bayesian", "km Ashauer","kem Bayesian", "kem Ashauer","km2 Bayesian", "km2 Ashauer","kem2 Bayesian", "kem2 Ashauer"), tick = TRUE, col = "black",
     lwd = 1, cex.axis=0.6, las=3)
errbar(c(0:11), c(kw,kwS,ke,keS,km,kmS,kemet,kemetS,km2,kmS2,kemet2,kemetS2), yplus=c(skw,skwS,ske,skeS,skm,skmS,skemet,skemetS,skm2,skmS2,skemet2,skemetS2), yminus=c(qkw,qkwS,qke,qkeS,qkm,qkmS,qkemet,qkemetS,qkm2,qkmS2,qkemet2,qkemetS2), add=TRUE, 
       lty=1, lwd=1,col=c("red","black","red","black","red","black","red","black","red","black","red","black"))
title(main="Ashauer's data - G. pulex exposed to 2,4,5-trichlorophenol")
abline(v=1.5, lty=2)
abline(v=3.5, lty=2)
abline(v=5.5, lty=2)
abline(v=7.5, lty=2)
abline(v=9.5, lty=2)

#####################################################################
#PPC
vx<-Cobs
vy<-qmed[c(1,18,18,31,32,80,81,100,101,125,126,164,164,246,246,409,409,500,500)]
vyinf<-qinf[c(1,18,18,31,32,80,81,100,101,125,126,164,164,246,246,409,409,500,500)]
vysup<-qsup[c(1,18,18,31,32,80,81,100,101,125,126,164,164,246,246,409,409,500,500)]
vxM<-Cobsmet
vyM<-qmedm[c(1,18,18,31,32,80,81,100,101,125,126,164,164,246,246,409,409,500,500)]
vyMinf<-qinfm[c(1,18,18,31,32,80,81,100,101,125,126,164,164,246,246,409,409,500,500)]
vyMsup<-qsupm[c(1,18,18,31,32,80,81,100,101,125,126,164,164,246,246,409,409,500,500)]
vxM2<-Cobsmet2
vyM2<-qmedm2[c(1,18,18,31,32,80,81,100,101,125,126,164,164,246,246,409,409,500,500)]
vyMinf2<-qinfm2[c(1,18,18,31,32,80,81,100,101,125,126,164,164,246,246,409,409,500,500)]
vyMsup2<-qsupm2[c(1,18,18,31,32,80,81,100,101,125,126,164,164,246,246,409,409,500,500)]


plot(vy~vx, col="darkgreen", pch=19, xlab="Observed", ylab="Predicted",ylim=c(0,18000), xlim=c(0,15000))
errbar(vx, vy, yplus=vysup, yminus=vyinf, add=TRUE, lty=1, lwd=1.5, errbar.col="darkgreen", cap=0)
errbar(vxM, vyM, yplus=vyMsup, yminus=vyMinf, add=TRUE, 
       lty=1, lwd=1.5,errbar.col="chartreuse3", cap=0)
errbar(vxM2, vyM2, yplus=vyMsup2, yminus=vyMinf2, add=TRUE, 
       lty=1, lwd=1.5,errbar.col="purple", cap=0)

points(vyM~vxM, col="chartreuse3", pch=19)
points(vyM2~vxM2, col="purple", pch=19)
abline(0,1)


#####################################################################
#Correlations of parameters 
require(GGally)
ggscatmat(mcmctotsample1)