Dear BUGS users:
I am having trouble fitting an IRT model in Linux with JAGS and wonder
if this is the appropriate forum to seek help. Suggestions will be
much appreciated.
I have used winBUGS (via R2winBUGS) for some time with success, but am
now confronting a model with thousands of parameters and a fairly
large dataset. Updating takes very long --- 5k iterations with 3
chains can take up to four hours, and I need many many more. Attempts
to use snow for parallel computation in the quad core processor have
failed beyond my understanding, so I abandoned that route.
I am now attempting to fit in Linux to compare speed using R's
system.time, and for this I migrated from BUGS to JAGS. My first Linux
trial run with 20 iterations appeared very promising --- the model
reports to run much much faster:
> Compiling model graph
> Resolving undeclared variables
> Allocating nodes
> Graph Size: 838547
>
> Initializing model
>
> user system elapsed
> 97.310 0.440 97.806
The equivalent triad in winBugs was about 241, .32, and 242, a
substantial improvement --- apparently. I then attempted a more
serious 10k burnin with 10k iterations. To my astonishment, the
reported system.time was virtually identical and the posterior sample
looks unmodified compared to the priors. There is, evidently, a
problem that is going unreported.
Can anyone offer clues about what may be going on and how I can
attempt to fix it? The model and key commands appear below for
reference.
Many thanks!
Eric Magar
Ciencia PolĂtica, ITAM
Mexico City
+52 55 5628 4079
@emagar
###########
## MODEL ##
###########
# rjags version
cat("
model {
for (j in 1:J){
for (i in 1:I){
v.star[j,i] ~ dnorm(mu[j,i],1)T(lo.v[j,i],hi.v[j,i]); ##
truncated normal sampling (BUGS uses I)
mu[j,i] <- beta[i]*x[j] - alpha[i] + delta[i]*y[j]
}
}
## priors ################
for (j in 1:J){
x[j] ~ dnorm(0, 1)
y[j] ~ dnorm(0, 1)
}
for(i in 1:(west-1)){
alpha[i] ~ dnorm( 0, 1)
beta[i] ~ dnorm( 0, 1)
delta[i] ~ dnorm( 0, 1)
}
alpha[west] ~ dnorm( 0, 1)
beta[west] ~ dnorm(-4, 20) ## informative
delta[west] ~ dnorm(-4, 20) ## informative
for(i in (west+1):(east-1)){
alpha[i] ~ dnorm( 0, 1)
beta[i] ~ dnorm( 0, 1)
delta[i] ~ dnorm( 0, 1)
}
alpha[east] ~ dnorm( 0, 1)
beta[east] ~ dnorm(-4, 20) ## informative
delta[east] ~ dnorm( 4, 20) ## informative
for(i in (east+1):(north-1)){
alpha[i] ~ dnorm( 0, 1)
beta[i] ~ dnorm( 0, 1)
delta[i] ~ dnorm( 0, 1)
}
alpha[north] ~ dnorm( 0, 1)
beta[north] ~ dnorm( 4, 20) ## informative
delta[north] ~ dnorm( 4, 20) ## informative
for(i in (north+1):I){
alpha[i] ~ dnorm( 0, 1)
beta[i] ~ dnorm( 0, 1)
delta[i] ~ dnorm( 0, 1)
}
}"
, file="model.jag")
#end model##############
##
dip.parameters <- c("delta","beta", "alpha", "x", "y")
dip.inits<-function() {
list(v.star=vstar, x = rnorm(J), y = rnorm(J), alpha = rnorm(I),
delta = rnorm(I), beta = rnorm(I))
}
dip.data <- list ("J"=J, "I"=I, "lo.v"=lo.v, "hi.v"=hi.v, "west"=west,
"east"=east, "north"=north)
##
#test ride to see program works
myWrap <- function(){
modelInit <- jags.model (
"model.jag",
inits=dip.inits,
data=dip.data,
n.chains=3,
n.adapt=10)
dip.60 <- coda.samples (modelInit,
dip.parameters,
n.iter=10,
thin=2)
return(dip.60)
}
print(system.time(myWrap()))
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|