### WinBUGS code to implement the Full [0-1-augmented-Beta] model

model
{
Cte<-1000000

for(i in 1:n)
{
b[i]~dnorm(0,tau)
}


for(j in 1:N)
  {
   zeros[j]        <- 0
   zeros[j]        ~  dpois(zeros.means[j])
   zeros.means[j]  <- -lBetaInf[j]+Cte
   fdBeta[j]       <- exp(loggam(a1[j]+a2[j])-loggam(a1[j])-loggam(a2[j])+(a1[j]-1)*log(Y[j])+(a2[j]-1)*log(1-Y[j]))
   a1[j]           <- mu[j]*phi
   a2[j]           <- (1-mu[j])*phi
   logit(mu1[j])   <- Beta[1]  + Beta[2]   * gender[j] + Beta[3]  * age[j] + Beta[4]  * hba1cd[j] + Beta[5]  * smoker[j] + Beta[6]  * incisor [j]+ Beta[7]*premolar[j]+Beta[8] * molar[j]+ b[cluster[j]]
   mu[j]           <- max(0.00001,min(0.9999,mu1[j]))   
   logit(p0[j])    <- gamma[1]  + gamma[2] * gender[j] + gamma[3]  * age[j] + gamma[4]  * hba1cd[j] + gamma[5]  * smoker[j] + gamma[6]  * incisor [j]+ gamma[7] * premolar[j] + gamma[8] * molar[j] 
   logit(p1[j])    <- rho[1]    + rho[2]   * gender[j] + rho[3]    * age[j] + rho[4]    * hba1cd[j] + rho[5]    * smoker[j] + rho[6]    * incisor [j]+ rho[7]   * premolar[j] + rho[8]   * molar[j] 
   e[j]            <- equals(Y[j],0.0001)
   d[j]            <- equals(Y[j],0.9999)
   fdBetaInf1[j]   <- (e[j]*p0[j]+d[j]*p1[j]+(1-e[j])*(1-d[j])*fdBeta[j]*(1-p0[j]-p1[j]))*step(1-p0[j]-p1[j])
   fdBetaInf[j]    <- max(0.00000001,fdBetaInf1[j])
   lBetaInf[j]     <- log(fdBetaInf[j])
   }

#Priors on parameters

for (i in 1:8)
  {
   Beta[i]  ~ dnorm(0,0.001)
   gamma[i] ~ dnorm(0,0.001)
   rho[i] ~ dnorm(0,0.001)
  }

phi~dgamma(0.1,0.01)
tau <- pow(sigmabsd,-2)
sigmabsd~dunif(0,100)
sigma2b<-pow(sigmabsd,2)
}