# Trap 66 in WinBUGS in a hierarchical Bayesian modeling

I want to analyze a multilevel multidimensional model in WinBUGS. the model is as below (N=2362 students responding to K=45 items of a test, students are nested within J=116 schools):

model{
#responses
for(i in 1:N){
for(j in 1:K){
logit(p[i,j])<- a1[j]*th[i,1]+a2[j]*th[i,2]-b[j]
y[i,j]~dbern(p[i,j] )
}
th[i,1:2]~dmnorm(mu[sc[i],1:2],tau.p[1:2,1:2])
}
#school level
for(j in 1:J){
mu[j,1:2]~dmnorm(m[j,1:2],tau.s[1:2,1:2])
}

#priors
for(j in 1:J){
m[j,1:2]~dmnorm(m0[1:2],cov[1:2,1:2])
}

tau.p[1:2,1:2]~dwish(cov[1:2,1:2],2)
tau.s[1:2,1:2]~dwish(cov[1:2,1:2],2)
sigma.p[1:2,1:2]<-inverse(tau.p[,])
sigma.s[1:2,1:2]<-inverse(tau.s[,])
s2p<-sum(sigma.p[,])
s2s<-sum(sigma.s[,])
rho<-(s2s)/(s2s+s2p)

a1~dlnorm(0,4)
a2<-0
b~dnorm(0,1)
for(s in 2:K) {
a1[s]~dlnorm(0,4)
a2[s]~dlnorm(0,4)
b[s]~dnorm(0,1)
}
}


I've set these functions as initial values:

ini<-function(){
list(tau.p=matrix(rgamma(4,100,100),2,2),
tau.s=matrix(rgamma(4,100,100),2,2),
th=rmvnorm(N,mean=c(0,0),sigma=diag(2)),
m=rmvnorm(J,mean=c(0,0),sigma=diag(2)),
mu=rmvnorm(J,mean=c(0,0),sigma=diag(2)),
a1=rlnorm(K,0, 0.4),
a2=c(NA,rlnorm(K-1,0, 0.4)),
b=rnorm(45,0,0.5))
}


I use rube package in R to check and run my analysis and everything looks fine. When I run the model I receive "Trap 66 (postcondition violated)" or "undefined real result". I think the problem is from the initials but I have no idea how to solve it.

Any idea?

• any chance you can provide a slimmed down data set, otherwise it is difficult for one to replicate your error and help you debug? Oct 10, 2013 at 16:07
• Just as an FYI you shouldn't use "t" as a variable name since it's the name of an internal R function ?t. Oct 10, 2013 at 22:58
• @gjabel it's a huge data set.I'm not sure how I can provide a slim data.Do you the problem due to data?
– Amin
Oct 11, 2013 at 16:25
– user37573
Jan 20, 2014 at 20:06
• I'd suggest trying JAGS. The error messages are far more transparent and you can use essentially the same code (maybe a small tweak here and there). Even if you are wedded to WinBUGS for some reason, you might get some clue as to what is wrong.
– guy
Jan 20, 2014 at 20:51

• How can I narrow down priors of variances while I'm using Wishart Dist.? The cov[1:2,1:2] is an identity matrix and I set df of Wishart Dist. as 2 because of the dimensions of th. I completely agree with you about trickiness of BUGS models.