Problems converting Openbugs Mixed Treatment Comparison Code to JAGS

Help
Steve Kay
2014-05-08
2014-05-08
  • 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

    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

    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

    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

     

Log in to post a comment.

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:

JavaScript is required for this form.





No, thanks