## WinBUGS code for "Bayesian Modeling of Multivariate Spatial Binary Data with applications to Dental Caries"
## by D. Bandyopadhyay, B.J.Reich and E. Slate  
## Statistics in Medicine: Special Issue on Methodology in Oral Health Research


#The binary observations are y[sub,surf,tooth], sub=1-nsubs,
#surf = 1-5 (1=Occlusal), tooth=1:32

#b[j1,j2] controls the association within tooth between surfaces
#j1 and j2.  For each sub, b is symmetric with zero on the diagonal.

#c[j1,j2] controls the assocation between surfaces j1 and j2 on
#neigboring teeth on the same jaw.  c is symmetric.

#d[j] controls the relationship between observations on surface j
#for neighboring teeth on different jaws.  d is set to zero expect for surface 1.




model{
 for(sub in 1:nsubs){
  for(surf in 1:5){
   for(tooth in 1:32){ d1[sub,surf,tooth] <- equals(y[sub,surf,tooth],1)
                       d2[sub,surf,tooth] <- equals(y[sub,surf,tooth],0)

        y[sub,surf,tooth]~dbern(prob[sub,surf,tooth]) 
        y1[sub,surf,tooth]~dbern(prob[sub,surf,tooth]) 
      

       prob[sub,surf,tooth] <- max(0.00001,min(0.99999,prob1[sub,surf,tooth]))
       
}
  


#teeth on the ends:
   logit(prob1[sub,surf,1])<-
       mu[sub]+inprod(b[surf,],y[sub,,1])+
       inprod(c[surf,],y[sub,,2])+d[surf]*y[sub,surf,32]
       
      
   logit(prob1[sub,surf,16])<-
       mu[sub]+inprod(b[surf,],y[sub,,16]) +
       inprod(c[surf,],y[sub,,15])+d[surf]*y[sub,surf,17]
                         
      
   logit(prob1[sub,surf,17])<-
       mu[sub]+inprod(b[surf,],y[sub,,17]) +
       inprod(c[surf,],y[sub,,18])+d[surf]*y[sub,surf,16] 
       

   logit(prob1[sub,surf,32])<-
       mu[sub]+inprod(b[surf,],y[sub,,32]) +
       inprod(c[surf,],y[sub,,31])+d[surf]*y[sub,surf,1]                    
         

   #middle teeth:
   for(tooth in 2:15){
    logit(prob1[sub,surf,tooth])<- mu[sub] +
       inprod(b[surf,],y[sub,,tooth])+inprod(c[surf,],y[sub,,tooth-1])+
       inprod(c[surf,],y[sub,,tooth+1])+d[surf]*y[sub,surf,33-tooth]
       

    logit(prob1[sub,surf,tooth+16])<-mu[sub] +
       inprod(b[surf,],y[sub,,tooth+16])+inprod(c[surf,],y[sub,,tooth+15])+
       inprod(c[surf,],y[sub,,tooth+17])+d[surf]*y[sub,surf,16+tooth]
       
    }
  }

 

### subject-level random effect
  mu[sub]~dnorm(mnmu[sub],precsub)
  mnmu[sub]<-inprod(x[sub,],beta[])
 }


## Priors on precision parameters

 precsub <- pow(sigmasub, -2)
 sigmasub ~ dunif(0,100)


### Prior on fixed effects

beta[1] ~ dflat()
for(j in 2:6){beta[j]~dnorm(0,0.25)
odds[j] <- exp(beta[j])
}


### Priors on spatial parameters and odds

 b[1,1]~dnorm(0,10000)
 c[1,1]~ddexp(0,tau)
 d[1] ~ ddexp(0,tau)

 for(surf in 2:5){
 
    b[surf,surf]~dnorm(0,10000)
    c[surf,surf]~ddexp(0,tau)
    d[surf]~dnorm(0,10000)
  }

  for(j1 in 2:5){for(j2 in 1:(j1-1)){
    b[j1,j2]<-b[j2,j1]; b[j2,j1]~ddexp(0,tau)
    c[j1,j2]<-c[j2,j1]; c[j2,j1]~ddexp(0,tau)
  }}



for (j1 in 1:5) { oddsd[j1] <- exp(d[j1])
     for (j2 in 1:5) {
                 oddsb[j1,j2] <- exp(b[j1,j2])
                 oddsc[j1,j2] <- exp(c[j1,j2])
                     }
                 }

## Prior on double-exponential scale parameter

tau <- pow(tau2, 0.5)
tau2 ~ dgamma(1,1)


}