Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

Close

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