## This WinBUGS code below is for a nested random effects model with 4 covariates for $n$ subjects with $m$ locations recorded for each subject. The subject random effects 
## $i$ are independent Bridge random variables, and the location ("tooth") random effects $\epsilon_{ij}$ are multivariate bridge, with exponential spatial correlation matrix 
## $\Sigma$. The constants $m$, $n$, the matrix of distances between locations (DIST), and the mathematical constant $\pi$ need to be entered as data in addition to X and Y.



model{
for (i in 1:n){ for (j in 1:m){
Y[i,j]~ dbern(pi[i,j])
logit(pi[i,j])<- alpha+X1[i,j]*beta1+ X2[i,j]*beta2+X3[i,j]*beta3+ X4[i,j]*beta4+gam[i]+eps[i,j]

# get PHI(theta), P(z<theta)
ptheta[i,j]<-phi(theta[i,j]) #phi is N(0,1) cdf
eps[i,j]<- log( sin(bphi*pi*ptheta[i,j])/
sin(bphi*pi*(1-ptheta[i,j])))/bphi2

}#endj


#Draw the subject random effects from Bridge distribution
u[i]~dunif(0,1)
gam[i]<- log( sin(bphi1*pi*u[i])/sin(bphi1*pi*(1-u[i])))/(bphi1*bphi2)

# Draw the latent location random effects from a multivariate normal distribution

theta[i,1:m]~dmvnorm(0, V[,])

} #end i


#hyperpriors
range~dlnorm(-2,1)
beta1~dnorm(0,btau)
beta2~dnorm(0,btau)
beta3~dnorm(0,btau)
beta4~dnorm(0,btau)


# Define precision matrix for theta, e.g. exponential
# (could also use spatial.exp instead of dmvnorm above)
for(k in 1:m){ for(l in 1:m){
Sigma[k,l]<-exp(-DIST[k,l]/range)
}}

V[1:m,1:m]<-inverse(Sigma[,])

# Prior for precision in RE distributions is Gamma.
tau1~dgamma(1,.125)
tau2~dgamma(1,.125)

# define the bridge parameters phi... use 'bphi' since 'phi' is a function.
bphi1<-1/sqrt(3/(tau1*pi^2)+1)
bphi2<-1/sqrt(3/(tau2*pi^2)+1)

}