library(mvtnorm)
library(tmvtnorm)
library(mnormt)
library(BRugs)
library(R2WinBUGS)

folderSNI<-"C:/Users/pc/Documents/privado/SkewEllipticall/PaperNonLinear/paperCensurasT/Bayesian-NIPaperCens/ProgramasNI-Cens"
setwd(folderSNI)

iqr5 <- function (x) {(quantile(x, 0.025))}
iqr95 <- function (x) {(quantile(x, 0.975))}
 
################################################################################
## Linear Mixed Model Data Vaida
################################################################################
 
 
library(lmec)
data(UTIdata)
UTIdata <- subset(UTIdata, !is.na(RNA))
o <- order(UTIdata$Patid, UTIdata$Fup)
UTIdata <- UTIdata[o,]
cens = (UTIdata$RNAcens==1)+0
y = log10(UTIdata$RNA)
x = cbind((UTIdata$Fup==0)+0, (UTIdata$Fup==1)+0, (UTIdata$Fup==3)+0, (UTIdata$Fup==6)+0, (UTIdata$Fup==9)+0, (UTIdata$Fup==12)+0, (UTIdata$Fup==18)+0, (UTIdata$Fup==24)+0)
z = matrix(rep(1, length(y)), ncol=1)
cluster = as.numeric(UTIdata$Patid)

nj<-matrix(0,72,1)
for (j in 1:72) {
nj[j]=sum(cluster==j)
}

y1<-y
y.cen<- rep(10000,sum(nj))
y.cen[cens==1]=y[cens==1]
y[cens==1]=NA

 ### Caso Normal
datas <- list ("y","x","y.cen","cluster")
modelN<-file.path(folderSNI, "BugsLMMvaida.r")
initsN <- function(){list(beta=rep(2,8),sigmae.inv=1,sigmab.inv=1)}
parametersN<-c("beta","sigmae","sigmab","mu","b","w")

niter=50000  # numero de iteraçoes no GIBBS   30000, 50000, 100000, etc

## AJusta um modelo Normal



resu.simN <- bugs(datas, initsN, parametersN, modelN,
    n.chains=2, n.iter=niter, n.thin=2 , debug=FALSE, codaPkg=FALSE,
    program = "openbugs" ,
    working.directory=NULL,digits = 5,clearWD = TRUE)

print(resu.simN,5)
 
 
 ### Caso Student
  
  
datas <- list ("y","x","y.cen","cluster")
modelT<-file.path(folderSNI, "BugsLMM-tvaida.r")
initsT <- function(){list(beta=rep(2,8),sigmae.inv=1,sigmab.inv=1,nu=8,beta1=0.3)}
parametersT<-c("beta","sigmae","sigmab","nu","mu","b","w")

niter=50000  # numero de iteraçoes no GIBBS   30000, 50000, 100000, etc



resu.simT <- bugs(datas, initsT, parametersT, modelT,
    n.chains=2, n.iter=niter, n.thin=2 , debug=FALSE, codaPkg=FALSE,
    program = "openbugs" ,
    working.directory=NULL,digits = 5,clearWD = TRUE)
    
     print(resu.simT,5)
     

     
     ### Caso Slash
  
  
datas <- list ("y","x","y.cen","cluster")
modelSL<-file.path(folderSNI, "BugsLMM-SLvaida.r")
initsSL<- function(){list(beta=rep(2,8),sigmae.inv=1,sigmab.inv=1,nu=8)}
parametersSL<-c("beta","sigmae","sigmab","nu","mu","b","w")

niter=50000  # numero de iteraçoes no GIBBS   30000, 50000, 100000, etc



resu.simSL <- bugs(datas, initsSL, parametersSL, modelSL,
    n.chains=2, n.iter=niter, n.thin=2 , debug=FALSE, codaPkg=FALSE,
    program = "openbugs" ,
    working.directory=NULL,digits = 5,clearWD = TRUE)
    
     print(resu.simSL,5)
 
 
      ### Caso CN
  
  
datas <- list ("y","x","y.cen","cluster")
modelCN<-file.path(folderSNI, "BugsLMM-CNvaida.r")
initsCN<- function(){list(beta=rep(2,8),sigmae.inv=1,sigmab.inv=1,nu=0.1,gamma=0.1)}
parametersCN<-c("beta","sigmae","sigmab","nu","gamma","mu","b","w")





resu.simCN <- bugs(datas, initsCN, parametersCN, modelCN,
    n.chains=2, n.iter=niter, n.thin=2 , debug=FALSE, codaPkg=FALSE,
    program = "openbugs" ,
    working.directory=NULL,digits = 5,clearWD = TRUE)
    
     print(resu.simCN,5)    
        
     
      
     
###Criterios

critT<-criterios(cens,y1,x,z,nj,espac=100,resu.simT,type="T")
critN<-criterios(cens,y1,x,z,nj,espac=100,resu.simN,type="Normal")

critSL<-criterios(cens,y1,x,z,nj,espac=100,resu.simSL,type="SL")
critCN<-criterios(cens,y1,x,z,nj,espac=100,resu.simCN,type="CN")

        critN$IL[41]
###Figuras Influence

par(mfrow = c(1, 2))
plot(critN$IL, ylab="K-L divergence", xlab="Index",type="h",main="Normal",ylim=range(0,1.5))
identify(critN$IL,n=4)
plot(critT$IL, ylab="K-L divergence", xlab="Index",type="h",main="Student-t",ylim=range(0,1.5))
plot(critSL$IL, ylab="K-L divergence", xlab="Index",type="h",main="Slash",ylim=range(0,1.5))
plot(critCN$IL, ylab="K-L divergence", xlab="Index",type="h",main="Contaminated Normal",ylim=range(0,1.5))
#identify(critT$IL,n=2)

### Residuals
# plot(a,b, main = substitute("Student "*{italic(t)}))
par(mfrow = c(1, 2))

qqnorm(resu.simN$mean$b)
qqline(resu.simN$mean$b)

qqnorm(resu.simT$mean$b)
qqline(resu.simT$mean$b)


qqnorm(y1-resu.simN$mean$mu)
qqline(y1-resu.simN$mean$mu)

qqnorm(y1[cens==0]-resu.simT$mean$mu[cens==0])
qqline(y1[cens==0]-resu.simT$mean$mu[cens==0])

####

#Graficos de 4 individuals

 ##Graficos
 
par(mfrow = c(2, 2))
j=10
plot(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))], y1[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],ylim=c(4.5, 6.5), main="ID = 10", ylab="Log10 (RNA) ", xlab="days",type="p")
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simN$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=1,col=1)
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simT$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=2,col=2)
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simSL$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=3,col=3)
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simCN$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=3,col=4)
points(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],cens[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],col="red",pch=19)

#abline(h=2,lty=5)
j=52
y22=y1[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
x22=UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
cens[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]

plot(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))], y22,ylim=c(2, 4), main="ID = 17", , ylab="Log10 (RNA) ", xlab="days",type="p")
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simN$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=1,col=1)
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simT$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=2,col=2)
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simSL$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=3,col=3)
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simCN$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=3,col=4)
points(x22[cens[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]==1],y22[cens[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]==1],col="red",pch=25)
#leg=c("Normal fit","Student-t fit","Slash fit")
#legend(locator(1),leg,lty=1:3,bty="n",col=1:3)
leg=c("Normal","Student-t","Slash", "Cont. Normal")
legend(locator(1),leg,lty=1:4,bty="n",col=1:4)

#abline(h=2,lty=5)


j=30

plot(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))], y1[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],ylim=c(2, 4), main="ID = 30", , ylab="Log10 (RNA) ", xlab="days",type="p")
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simN$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=1,col=1)
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simT$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=2,col=2)
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simSL$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=3,col=3)
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simCN$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=3,col=4)
#leg=c("Normal fit","Student-t fit","Slash fit")

#legend(locator(1),leg,lty=1:3,bty="n",col=1:3)
#abline(h=2,lty=5)

j=70

plot(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))], y1[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],ylim=c(3,5), main="ID = 70", , ylab="Log10 (RNA) ", xlab="days",type="p")
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simN$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=1,col=1)
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simT$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=2,col=2)
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simSL$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=3,col=3)
lines(UTIdata$Fup[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],resu.simCN$mean$mu[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))],lty=3,col=4)



### Graficos de Nu

#graficos Densidades
par(mfrow = c(2, 2))
dT<-attach.bugs(resu.simT,overwrite = FALSE)
hist(dT$nu,xlab=expression(nu), main="(a) Student-t",ylab="Density",xlim=range(1,6),probability=TRUE,col="grey")
lines(density(dT$nu,bw=0.1))
abline(v=iqr5(dT$nu),lty=3)
abline(v=iqr95(dT$nu),lty=3)

#### BI

#hist(resu.simT$mean$b[cens==0],xlab=expression(nu), main="Student-t",ylab="Density",xlim=range(-5,5),probability=TRUE,col="grey")
#lines(density(resu.simT$mean$bi[cens==0],bw=0.1))
#abline(v=iqr5(dT$nu),lty=3)
#abline(v=iqr95(dT$nu),lty=3)
########
  
leg=c(paste("[", round(iqr5(dT$nu),2),";",round(iqr95(dT$nu),2),"]"))
legend(locator(1),leg,bty="n")

expression(paste("BS(", alpha,"=1,",beta,"=1)"))

#dSL<-attach.bugs(resu.simSL,overwrite = FALSE)
hist(dSL$nu,,xlab=expression(nu),main="(b) Slash",ylab="Density",xlim=range(0.8,1.5),probability=TRUE,col="grey")
lines(density(dSL$nu,bw=0.015))
abline(v=iqr5(dSL$nu),lty=3)
abline(v=iqr95(dSL$nu),lty=3)
leg=c(paste("[", round(iqr5(dSL$nu),2),";",round(iqr95(dSL$nu),2),"]"))
legend(locator(1),leg,bty="n")


#dCN<-attach.bugs(resu.simCN,overwrite = FALSE)                               
hist(dCN$nu,,xlab=expression(nu),main="(c) Contaminated Normal",ylab="Density",probability=TRUE,col="grey")
lines(density(dCN$nu,bw=0.01))
abline(v=iqr5(dCN$nu),lty=3)
abline(v=iqr95(dCN$nu),lty=3)
leg=c(paste("[", round(iqr5(dCN$nu),2),";",round(iqr95(dCN$nu),2),"]"))
legend(locator(1),leg,bty="n")


hist(dCN$gamma,,xlab=expression(rho),main="(d) Contaminated Normal",ylab="Density",probability=TRUE, col="grey")
lines(density(dCN$gamma,bw=0.01))
abline(v=iqr5(dCN$gamma),lty=3)
abline(v=iqr95(dCN$gamma),lty=3)              
leg=c(paste("[", round(iqr5(dCN$gamma),2),";",round(iqr95(dCN$gamma),2),"]"))
legend(locator(1),leg,bty="n")


### Figura Pesos

#par(mfrow = c(2, 1))
#dT<-attach.bugs(resu.simT,overwrite = FALSE)
M5<-matrix(0,72,1)
M95<-matrix(0,72,1)
for(j in 1:72){
M5[j]<-iqr5(dT$w[,j])
M95[j]<-iqr95(dT$w[,j])
}
aa=resu.simT$mean$w
M95[72]<-M95[72]+0.1
M5[72]<-M5[72]+0.1
aa[72]<-aa[72]+0.1
plot(aa,type="p",ylim=c(-1,6),xlab="Index",ylab="u",main="Student-t")
points(M5,type="l")
points(M95,type="l")
identify(M5,n=4)
abline(h=M5[20],lty=3)
points(x=c(20,35,42,57),y=resu.simT$mean$w[c(20,35,42,57)],col="red",pch=19)


#dSL<-attach.bugs(resu.simSL,overwrite = FALSE)
M5SL<-matrix(0,72,1)
M95SL<-matrix(0,72,1)
for(j in 1:72){
M5SL[j]<-iqr5(dSL$w[,j])
M95SL[j]<-iqr95(dSL$w[,j])
}

plot(resu.simSL$mean$w,type="p",ylim=c(-1,2),xlab="Index",ylab="u",main="Slash")
points(M5SL,type="l")
points(M95SL,type="l")
identify(M5SL,n=3)
abline(h=M5SL[20],lty=3)