Steve Kay
2014-05-08
Hi,
NICE produce some great docs on Mixed Treatment Comparisons - search "NICE TSD" in google. The code provided is all in Winbugs. I'm trying to convert to JAGS to take advantage of running parallel chains on multi-processors using runjags.
Here's my R jags code for example 5 in NICE's TSD2:
library(runjags)
jags.MTC <- "model{
for(i in 1:ns){ # LOOP THROUGH STUDIES
mu[i] ~ dnorm(0,.0001) # vague priors for all trial baselines
for (k in 1:na[i]) { # LOOP THROUGH ARMS
vari[i,k] <- pow(se[i,k],2) # calculate variances
prec[i,k] <- 1/vari[i,k] # set precisions
y[i,k] ~ dnorm(theta[i,k],prec[i,k]) # normal likelihood
theta[i,k] <- mu[i] + d[t[i,k]] - d[t[i,1]] # model for linear predictor
}
}
d[1]<-0 # treatment effect is zero for reference treatment
for (k in 2:nt){ d[k] ~ dnorm(0,.0001) } # vague priors for treatment effects
} "
# Initial Values
inits1 <- list("d"=c(NA, 0,0,0,0), "mu"=c(0, 0, 0, 0, 0, 0, 0))
inits2 <- list("d"=c(NA, -1,-3,-1,1), "mu"=c(-3, -3, -3, -3, -3, -3, -3))
inits3 <- list("d"=c(NA, 2,2,2,2), "mu"=c(-3, 5, -1, -3, 7, -3, -4))
inits <- list(inits1, inits2, inits3)
t<-matrix(nrow=7,ncol=3)
y<-matrix(nrow=7,ncol=3)
se<-matrix(nrow=7,ncol=3)
t[,1] <- c(1,1,1,3,3,4,4)
t[,2] <- c(3,2,2,4,4,5,5)
t[,3] <- c(NA,NA,4,NA,NA,NA,NA)
y[,1] <- c(-1.22,-0.7,-0.3,-0.24,-0.73,-2.2,-1.8)
y[,2] <- c(-1.53,-2.4,-2.6,-0.59,-0.18,-2.5,-2.1)
y[,3] <- c(NA,NA,-1.2,NA,NA,NA,NA)
se[,1] <- c(0.504,0.282,0.505,0.265,0.335,0.197,0.2)
se[,2] <- c(0.439,0.258,0.51,0.354,0.442,0.19,0.25)
se[,3] <- c(NA,NA,0.478,NA,NA,NA,NA)
na <- c(2,2,3,2,2,2,2)
datalist <- list(nt=5,t=t,y=y,se=se,ns=7,na=na)
results <- run.jags(model=jags.MTC, n.chains=3, inits=inits,sample = 20000,burnin = 10000,
method="rjparallel", data=datalist,monitor=c("d"))
It compiles then sends out the following message
"Starting the rjags simulations using a socket cluster with 3 nodes on host ‘localhost’
Error in checkForRemoteErrors(val) :
3 nodes produced errors; first error: RUNTIME ERROR:
Compilation error on line 11.
d[1] is a logical node and cannot be observed"
As my programming seems in line with the advice on page 8 of the JAGS 3.4 manual as regards setting parameters to 0 and using NA in the initial values I'm not sure what is it I'm doing wrong. Any help much appreciated.
Steve
Matt Denwood
2014-05-08
Hi Steve
There appears to be an issue with the initial values being passed through using the 'rjparallel' method of runjags - your code works fine for me using the other JAGS interface methods. I'll take a closer look hopefully tomorrow and get back to you with what's going on, but in the meantime try using method='parallel' or method='snow'
Cheers,
Matt
Martyn Plummer
2014-05-08
Hi Matt. Where are you these days?
It's not the initial values but the replication of the data set. The data() member function for jags model objects should return a list of observed nodtes that can be read back into a new model. But in this case it adds an unwanted component for d
, i.e.
> newdata <- m$data() > newdata$d [1] 0 NA NA NA NA
When you try to use this as data for a new model it fails because d[1]
is already defined in the model as a constant node.
I knew about this problem, and it is already fixed in the development version. Unfortunately it can't be fixed in JAGS 3.x.y.
Steve, a workaround is to remove the definition of d[1] <- 0
from the model and include it instead as data, i.e.
d <- c(0, NA, NA, NA, NA) datalist <- list(nt=5,t=t,y=y,se=se,ns=7,na=na,d=d )
Matt Denwood
2014-05-08
Hi Martyn
Thanks for the tip - I would have been on a wild goose chase trying to track that one down I think!
I've just moved to University of Copenhagen (as of last week). It's early days yet, but I am optimistic that this will give me more time to spend on my core research activities - as well as playing around with R code of course :)
Matt
Steve Kay
2014-05-08
Brilliant - many thanks to both of you for your prompt replies (I've spent many hours trying different ways to solve this - would have got to your solution Martyn in another 100 hours!)
Best,
Steve