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[1]~dlnorm(0,4)
a2[1]<-0
b[1]~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?