## R2WinBUGS code from 
## Mustvari T, Bandyopadhyay D, Lesaffre E and  Declerck D. (2013). "A multilevel model for spatially correlated binary data in the presence of misclassification: 
## An application in oral health research",  Statistics in Medicine, 32(30), 5241-5259
## Written by T. Mustvari



#### Setting up the directory, and reading datasets

wd<-"D:/PhD/New CAR"
setwd(wd);
datos<-read.csv("cent_2gamma.csv",header=TRUE,na.string="");
attach(datos);


### Generating initial values from a simple GLM fit

y<-ES_status;

xcen1<-(xcen-mean(xcen))/sd(xcen)

ycen1<-(ycen-mean(ycen))/sd(ycen)
AGE1<-(AGE-mean(AGE))/sd(AGE)

G1<-ifelse(GENDER=='girl',1,0);

T1<-ifelse(Type==1,1,0);

TT2<-ifelse(TType==2,1,0);
TT3<-ifelse(TType==3,1,0);
TT1<-ifelse(TType==1,1,0);

model1<-glm(ES_status~ as.factor(GENDER) + AGE1 + as.factor(Type) + as.factor(jaw) + as.factor(quad) + xcen1 + ycen1, 
            family=binomial(link="logit"),data=datos);

##print(confint(model1, level = 0.95),2)

######
beta<-model1$coefficients;names(beta)<-NULL;
beta<-c(-6.542,-0.274, 0.006,-0.989, 1.355, 3.644,-0.658, 0.071, 0.168, 0.081, 0.119)


########
### Creating some other variables for later input in the WinBUGS program

t1<-unique(tid);
tid2<-rep(0,length(t1));
for (i in 1:length(t1)){tid2<-ifelse(tid==t1[i],i,tid2)}

id1<-unique(idnr2);
idd2<-rep(0,length(id1));
for (i in 1:length(id1)){idd2<-ifelse(idnr2==id1[i],i,idd2)}

ex1<-unique(EXAMINER);
ex2<-rep(0,length(ex1));
for (i in 1:length(ex1)){ex2<-ifelse(EXAMINER==ex1[i],i,ex2)}

N<-nrow(datos);
n1<-length(unique(idd2));
n2<-length(unique(tid2));
n3<-length(unique(ex2));


#######################

bm<-rnorm(n1,0,1);
bt<-rnorm(n2,0,1);
be<-rnorm(n3,0,1);



#######################
### The dataset to be inputted

data<-list("N","n1","n2","n3","y","G1","T1","surface1","surface2","surface3","surface4","surfnext","num","num2","jaw","quad","idd2","tid2","ex2","xcen","ycen","AGE");

### Initial values
inits <- function(){list(beta=beta,bm=bm,sigma.m=1, bt=bt,sigma.t=1,be=be,sigma.e=1,gamma1=1, gamma2=1)}

### Parameters
parameters  <- c("beta","gamma1","gamma2","sigma2.m","sigma2.t","sigma2.e");


########################
### Calling the model

cat("model
         {
          # Modeling at the observation level 
          for(i in 1:N)
           { 
            y[i]~dbern(p[i])


       mu[i] <-beta[1] + beta[2]*G1[i] + beta[3]*(AGE[i]-mean(AGE[]))/sd(AGE[]) + beta[4]*T1[i] + beta[5]*jaw[i] + beta[6]*quad[i] + 
               beta[7]*(xcen[i]-mean(xcen[]))/sd(xcen[]) + beta[8]*(ycen[i]-mean(ycen[]))/sd(ycen[]) +  
               bm[idd2[i]] + bt[tid2[i]] + be[ex2[i]] 
             
       B1[i] <- gamma1*equals(num[i],1)*(surface1[i] + surface2[i] + surface3[i])
       B2[i] <- gamma1*equals(num[i],2)*(surface1[i] + surface2[i] + surface3[i])
       B3[i] <- gamma1*equals(num[i],3)*(surface1[i] + surface2[i])
       B4[i] <- gamma1*equals(num[i],4)*(surface1[i] + surface2[i] + surface3[i])
       B5[i] <- gamma1*equals(num[i],5)*(surface1[i] + surface2[i] + surface3[i] + surface4[i])
       B6[i] <- gamma2*equals(num2[i],1)*(surfnext[i])


       logit(p[i])<-mu[i] + B1[i] +  B2[i] + B3[i] + B4[i] + B5[i] + B6[i]

       #a[i]<-exp(mu[i])/(1+exp(mu[i]))

       #p[i] <- eta[i]*t11 +  (1-eta[i])*t10


           }  

        # Level: Mouth
        for(m1 in 1:n1)
          {
             bm[m1]~dnorm(0,tau.m)
          }
          

        # Level: Tooth
        for(t1 in 1:n2)
          {
             bt[t1]~dnorm(0,tau.t)
          }
          

        # Level: Examiner
        for(e1 in 1:n3)
          {
             be[e1]~dnorm(0,tau.e)
          }
          

       ### Priors for fixed effects
       for (l in 1:8) { beta[l] ~ dnorm(0.0, 1.0E-6) }


       ### Other priors
       sigma.m ~ dunif(0.0,100)
       sigma2.m <- pow(sigma.m,2)
       tau.m<-1/sigma2.m
       
       sigma.t ~ dunif(0.0,100)
       sigma2.t <- pow(sigma.t,2)
       tau.t<-1/sigma2.t

       sigma.e ~ dunif(0.0,100)
       sigma2.e <- pow(sigma.e,2)
       tau.e<-1/sigma2.e

 
       gamma1 ~ dunif(0,3)
       gamma2 ~ dunif(0,3)}", file="model2f.txt")


library(R2WinBUGS);
library(R2jags);
library(BRugs);

## Calling WinBUGS inside R
model1.sim <- bugs(data,inits,parameters,"model2f.txt",n.chains=3,n.iter=15000,
                   n.burnin=10,n.thin=1, bugs.directory="c:/Program Files/WinBUGS14/",debug=TRUE,
                   working.directory=wd, clearWD=FALSE);
a<-proc.time()


## You may also call JAGS
model2.sim <- jags2(data,inits,parameters,"model2f.txt",n.chains=3,n.iter=10000,n.burnin=100, working.directory=wd);
c<-proc.time()


## Working with the diagnostics
jagsfit.mcmc <- as.mcmc.list(model2.sim)
traceplot(jagsfit.mcmc)
xyplot(jagsfit.mcmc)
densityplot(jagsfit.mcmc)


#### ROC Curves ######
a<-model2.sim$mean$pred
dd<-data.frame(cbind(y,a))
write.foreign(dd, "c:/PhD/mydata.txt", "c:/PhD/mydata.sas",   package="SAS")