######################################################################################################## ## The codes below use the zeroes trick available in WinBUGS to sample from (i) Simplex, (ii) Beta rectangular, ## and (iii) Beta densities. ######################################################################################################## ######################################################################################################## #(i) Simplex regression model ######################################################################################################## model { C <-1000 pi<-3.1416 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] <- -lDispInf[j]+C a[j] <- 2*pi*sigma2*pow((Y[j]*(1-Y[j])),3) fdDispersion[j] <- pow((a[j]),(-0.5))* exp((-0.5*pow((Y[j]-mu[j]),2))/(sigma2*Y[j]*(1-Y[j])*pow(mu[j]*(1-mu[j]),2))) 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]) <- psi[1] + psi[2] * gender[j] + psi[3] * age[j] + psi[4] * hba1cd[j] + psi[5] * smoker[j] + psi[6] * incisor [j] + psi[7] * premolar[j] + psi[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.01) d[j] <- equals (Y[j] , 0.99) fdDispInf1[j] <- (e[j] * p0[j] + d[j] * p1[j] + ( 1-e[j]) * (1-d[j]) * fdDispersion[j] * (1-p0[j]-p1[j]) )*step(1-p0[j]-p1[j]) fdDispInf[j] <- max(0.00000001,fdDispInf1[j]) lDispInf[j] <- log(fdDispInf[j]) } #Priors for the parameters in the model for (i in 1:8) { Beta[i] ~ dnorm(0,0.001) psi[i] ~ dnorm(0,0.001) rho[i] ~ dnorm(0,0.001) } tau <- pow(sigmabsd,-2) sigmabsd ~ dunif(0,100) sigma2b <- pow(sigmabsd,2) sigma ~ dunif(0,100) sigma2 <- pow(sigma,2) phi <- 1/sigma2 } ######################################################################################################## #(ii) Beta rectangular regression model ######################################################################################################## model { C <-1000 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]+C v[j] ~ dbern(theta[j]) u[j] <- v[j]+1 a1[j,1] <- mu[j]*phi a2[j,1] <- (1-mu[j])*phi a1[j,2] <- 1 a2[j,2] <- 1 logit(gamma1[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]] gamma[j] <- max(0.0001,min(0.9999,gamma1[j])) theta1[j] <- alpha * (1-abs(2*gamma[j]-1)) theta[j] <- max(0.0001,min(0.9999,theta1[j])) mu1[j] <- (gamma[j]-0.5*theta[j])/(1-theta[j]) mu[j] <- max(0.001,min(0.999,mu1[j])) fdBeta[j] <- exp( loggam ( a1 [ j, u[j] ] + a2 [ j,u[j] ] ) - loggam ( a1 [ j, u[j] ] )- loggam ( a2 [ j, u[j] ] ) + ( a1[ j, u[j] ] - 1 ) * log ( Y [j] ) + ( a2 [j, u[j] ] - 1 ) * log ( 1 - Y [j] ) ) d[j] <- equals(Y[j],0.01) e[j] <- equals(Y[j],0.99) logit(p0[j]) <- psi[1] + psi[2] * gender[j] + psi[3] * age[j] + psi[4] * hba1cd[j] + psi[5] * smoker[j] + psi[6] * incisor [j] + psi[7] * premolar[j] + psi[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] fdBetaInf1[j] <- (p0[j]*d[j]+(1-d[j])*(1-e[j])*fdBeta[j]*(1-p0[j]-p1[j])+p1[j]*e[j])*step(1-p0[j]-p1[j]) fdBetaInf[j] <- max(0.00000001,fdBetaInf1[j]) lBetaInf[j] <- log(fdBetaInf[j]) } #Priors for the parameters in the model for (i in 1:8) { Beta[i] ~ dnorm(0,0.001) psi[i] ~ dnorm(0,0.001) rho[i] ~ dnorm(0,0.001) } phi ~ dgamma(0.1,0.01) tau <- pow(sigmab,-2) sigmab ~ dunif(0,100) sigma2b <- pow(sigmab,2) alpha ~ dunif(0,1) } ######################################################################################################## #(iii) Beta regression model ######################################################################################################## model { C <-1000 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]+C 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.01) d[j] <- equals(Y[j],0.99) 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 for the parameters in the model 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) }