Problems converting Openbugs Mixed Treatment Comparison Code to JAGS

Steve Kay

  • Steve Kay

    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:

    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[,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.


  • Matt Denwood
    Matt Denwood

    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'



  • Martyn Plummer
    Martyn Plummer

    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
    Matt Denwood

    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 :)



  • Steve Kay

    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!)