Thanks for the reply Bob, but it still does not work, you see. I ran this
model, just with the main effects and it ran fine.
n=length(bi.bmi)
Lgen=2
Lrace=5
Lagegp=13
Lstra=15
Lpsu=2
bi.bmi.model=function(){
# likelihood
for(i in 1:n){
bi.bmi[i]~ dbern(p[i])
logit(p[i])<- a0 + a1[agegp[i]]+a2[gen[i]]+a3[race[i]]
+ g[stra[i]]+ u[psu[i],stra[i]]
}
# constraints for a1, a2, a3
a1[1]<-0.0
a2[1]<-0.0
a3[1]<-0.0
# priors
a0~ dnorm(0.0, 1.0E-4)
for(j in 2:Lagegp){a1[j]~ dnorm(0.0, 1.0E-4)}
for(j in 2:Lgen){ a2[j]~ dnorm(0.0, 1.0E-4)}
for(k in 2:Lrace){ a3[k]~ dnorm(0.0, 1.0E-4)}
for(l in 1:Lstra){
g[l]~dunif(0, 100)
}
for( m in 1:Lpsu){
for(l in 1:Lstra){
u[m,l]~ dnorm(0.0, tau.u)
}}
tau.u<-pow(sigma.u, -2)
sigma.u~ dunif(0.0,100)
}
library(BRugs)
writeModel(bi.bmi.model, con='bi.bmi.model.txt')
model.data=list( 'n','Lagegp', 'Lgen', 'Lrace', 'Lstra', 'Lpsu',
'bi.bmi','agegp', 'gen', 'race','stra', 'psu')
model.init=function(){
list( sigma.u=runif(1),
a0=rnorm(1), a1=c(NA, rep(0,12)),
a2=c(NA, rep(0, 1)),
a3=c(NA, rep(0, 4)),
g=rep(0,Lstra), u=matrix(rep(0, 30), nrow=2) )
}
model.parameters=c( 'a0', 'a1', 'a2', 'a3')
model.bugs=BRugsFit(modelFile='bi.bmi.model.txt',
data=model.data,
inits=model.init,
numChains=1,
para=model.parameters,
nBurnin=50, nIter=100)
This is just with the main effects, and this does not give me any problems,
and I also ran the following model with interaction term between gen and
race, and it also ran fine.
for (i in 1:n){
bi.bmi[i]~ dbern(p[i])
logit(p[i])<- a0 + a1[agegp[i]]+a2[gen[i]]+a3[race[i]]
+ a23[gen[i], race[i]]
+ gam[stra[i]]+ u[psu[i],stra[i]]
}
# constraints for a2, a3, a12 and a13
a1[1]<-0.0
a2[1]<-0.0
a3[1]<-0.0
a23[1,1]<-0.0
#gen x race
for(j in 2:Lrace){ a23[1,j]<-0.0}
for(k in 2:Lgen){ a23[k,1]<-0.0}
# priors
a0~ dnorm(0.0, 1.0E-4)
for(i in 2:Lagegp){a1[i]~dnorm(0.0, 1.0E-4)}
for(i in 2:Lgen){ a2[i]~ dnorm(0.0, 1.0E-4)}
for(i in 2:Lrace){ a3[i]~ dnorm(0.0, 1.0E-4)}
for(i in 2:Lgen){
for(j in 2:Lrace){
a23[i,j]~ dnorm(0.0, 1.0E-4)
}}
for(i in 1:Lstra){
gam[i]~dunif(0, 1000)
}
for( i in 1:Lpsu){
for(j in 1:Lstra){
u[i,j]~ dnorm(0.0, tau.u)
}}
tau.u<-pow(sigma.u, -2)
sigma.u~ dunif(0.0,100)
}
So, the error happens only when I try to plug in interaction with the agegp.
I still don't know how to correct it.
Thanks
