Rewriting WinBUGS code into JAGS code

2014-08-04
2014-08-08
  • Inga Schwabe

    Inga Schwabe - 2014-08-04

    Dear all,
    I am trying to rewrite code that was written in WinBUGS to JAGS code.
    The code I want to use is the following:

    model{
        for( i in 1 : nmz ){
            g2mz[i, 1:N] ~ dmnorm(mu[1:N],tau.g[1:N, 1:N]) 
            g1mz[i, 1:N] ~ dmnorm(g2mz[i, 1:N],tau.g[1:N, 1:N])
            esmz[i] <- g1mz[i, 2]
            eemz[i] <- g1mz[i, 1]
    
            for(j in 1:2){
                ermz[i, j] <- g1mz[i, 3] + esmz[i] * emz[i, j] 
                emz[i, j] ~ dnorm(eemz[i],tau.e)
                rmz[i, j] ~ dnorm(ermz[i, j],tau.r)
                smz[i, j] ~ dnorm(esinmz[i],tau.s)
                esinmz[i] <- beta[1] + beta[2] * esmz[i]
            }   
        }
    
        for (k in 1:N){ 
            for (l in 1:N){
                sigma2.g[k,l] <- inverse(tau.g[,],k,l)*2 # sigma2.g is genetic covariance matrix 
            } 
        }
    
    #Priors
    tau.g[1:N, 1:N] ~ dwish(omega.g[,], N) 
    mu[1:N] ~ dmnorm(mean[],precis[,])
    tau.e ~ dgamma(0.001,0.001) 
    tau.s ~ dgamma(0.001,0.001)
    tau.r ~ dgamma(0.001,0.001)
    beta[1:2] ~ dmnorm(mbeta[],prbeta[ , ])
    }
    

    The code is adapted from the article 'Markov chain Monte Carlo Approaches to Analysis of Genetic and Environmental Components of Human Development and G x E Interaction' by Eaves and Erkanli (2003, Behavior Genetics). However, when I try to run the code, I am getting the following error message:
    Compilation error on line 3.
    Unable to resolve node mu[1:3]
    This may be due to an undefined ancestor node or a directed cycle in the graph
    While searching the internet for any helpful resources, I found out that in JAGS, the multivariate normal distribution has to be specified as dmnorm(mu, Omega) (http://www.jrnold.me/blog/jagsopenbugs-to-stan-distributions.html), so I changed the code to the following:

    model{
        for( i in 1 : nmz ){
            g2mz[i, 1:N] ~ dmnorm(mu,tau.g)  ## CHANGED CODE !
            g1mz[i, 1:N] ~ dmnorm(g2mz, tau.g) ## CHANGED CODE !
            esmz[i] <- g1mz[i, 2]
            eemz[i] <- g1mz[i, 1]
    
            for(j in 1:2){
                ermz[i, j] <- g1mz[i, 3] + esmz[i] * emz[i, j] 
                emz[i, j] ~ dnorm(eemz[i],tau.e)
                rmz[i, j] ~ dnorm(ermz[i, j],tau.r)
                smz[i, j] ~ dnorm(esinmz[i],tau.s)
                esinmz[i] <- beta[1] + beta[2] * esmz[i]
            }   
        }
    
        for (k in 1:N){ 
            for (l in 1:N){
                sigma2.g[k,l] <- inverse(tau.g[,],k,l)*2 # sigma2.g is genetic covariance matrix 
            } 
        }
    
    #Priors
    tau.g[1:N, 1:N] ~ dwish(omega.g[,], N) 
    mu[1:N] ~ dmnorm(0, .1)  # ## CHANGED CODE !
    
    tau.e ~ dgamma(0.001,0.001) 
    tau.s ~ dgamma(0.001,0.001)
    tau.r ~ dgamma(0.001,0.001)
    beta[1:2] ~ dmnorm(mbeta[],prbeta[ , ])
    }
    

    However, I am still getting the same error. I have no experience in using multivariate distributions in JAGS, could anyone point me to the error?
    I believe that the error is somewhere in the specification of the distribution. However, if it is helpful to have any example R code to run the model, please let me know!
    Many thanks,
    Inga Schwabe

     
  • Martyn Plummer

    Martyn Plummer - 2014-08-05

    Did you supply mean (a vector of length 3) and precis (a 3x3 matrix) as data when you run your model?

     
    • Inga Schwabe

      Inga Schwabe - 2014-08-06

      Dear Martyn,

      Thanks a lot for our quick response.

      Yes, I have, but your reaction pointed me to another problem, I did not supply omega.g as data when I was running the model. I am sorry for the trivial error!

      However, when I am supplying data for omega.g I am getting the following error:

      Error in jags.model(model.file, jagsdata, inits, n.chains = 1, quiet = FALSE) :
      RUNTIME ERROR:
      Incorrect number of parameters in function inverse

      Is the dimension of the omega.g data incorrect?

      This is part of the code I am using in R to specify my data:

      mu <- c(0,0,0)
      precis <- matrix(c(1, 1/10, 1/10, 1/10, 1, 1/10, 1/10, 1/10, 1), 3,3)
      omega.g <- diag(3) * var_a #var_a =0.5, 3 = N (number of variables, as in JAGS code)
      
      jagsdata <- list(Ymz, Emz, Smz, nMZ, N, mu, precis, omega.g)
      names(jagsdata)<- c("rmz", "emz", "smz", "nmz", "N", "mu", "precis", "omega.g") 
      model.file <- "eaves.bug"
      inits <- NULL
      

      I changed the specification of beta[1:2] ~ dnorm(mbeta[], prbeta[,]) to beta_0 ~ dnorm(0,.1) and beta_1(0,.1), so I do not have to provide any data for that prior. To provide any misunderstandings, this is the current version of the syntax I am using:

      model
      {
          for( i in 1 : nmz ){
              g2mz[i, 1:N] ~ dmnorm(mu[1:N],tau.g[1:N, 1:N]) 
              g1mz[i, 1:N] ~ dmnorm(g2mz[i, 1:N],tau.g[1:N, 1:N])
              esmz[i] <- g1mz[i, 2]
              eemz[i] <- g1mz[i, 1]
      
              for(j in 1:2){
                  ermz[i, j] <- g1mz[i, 3] + esmz[i] * emz[i, j] 
                  emz[i, j] ~ dnorm(eemz[i],tau.e)
                  rmz[i, j] ~ dnorm(ermz[i, j],tau.r)
                  smz[i, j] ~ dnorm(esinmz[i],tau.s)
                  esinmz[i] <- beta_1 + beta_2 * esmz[i]
              }   
          }
      
          for (k in 1:N){ 
              for (l in 1:N){
                  sigma2.g[k,l] <- inverse(tau.g[,],k,l)*2 # sigma2.g is genetic covariance matrix 
              } 
          }
      
      #Priors
      tau.g[1:N, 1:N] ~ dwish(omega.g[,], N) 
      mu[1:N] ~ dmnorm(mean[],precis[,])
      tau.e ~ dgamma(0.001,0.001) 
      tau.s ~ dgamma(0.001,0.001)
      tau.r ~ dgamma(0.001,0.001)
      beta0 ~ dnorm(0, .1)
      beta1 ~ dnorm(0, .1)
      
      }
      

      Many thanks.

       
  • Martyn Plummer

    Martyn Plummer - 2014-08-06

    OK I see the problem. Early versions of WinBUGS did not have a proper matrix inverse function. It would only return a scalar, and you needed to specify which element of the inverse matrix you wanted with second and third arguments to inverse(). To get the whole inverse matrix you needed to construct it element by element like this:

    for (k in 1:N){ 
       for (l in 1:N){
          sigma2.g[k,l] <- inverse(tau.g[,],k,l)*2 # sigma2.g is genetic covariance matrix 
       } 
    }
    

    JAGS has a proper inverse function that returns the whole inverse matrix, so you can replace those nested loops with

    sigma2.g <- inverse(tau.g) * 2 # sigma2.g is genetic covariance
    

    This will also work in OpenBUGS. Your problem stems from the fact that you are adapting some quite old BUGS code that is now obsolete.

     
    • Inga Schwabe

      Inga Schwabe - 2014-08-07

      Thanks for your quick response.

      I changed the code to the following (using the inverse function):

      model
      {
          for( i in 1 : nmz ){
              g2mz[i, 1:N] ~ dmnorm(mu[1:N],tau.g[1:N, 1:N]) 
              g1mz[i, 1:N] ~ dmnorm(g2mz[i, 1:N],tau.g[1:N, 1:N])
              esmz[i] <- g1mz[i, 2]
              eemz[i] <- g1mz[i, 1]
      
              for(j in 1:2){
                  ermz[i, j] <- g1mz[i, 3] + esmz[i] * emz[i, j] 
                  emz[i, j] ~ dnorm(eemz[i],tau.e)
                  rmz[i, j] ~ dnorm(ermz[i, j],tau.r)
                  smz[i, j] ~ dnorm(esinmz[i],tau.s)
                  esinmz[i] <- beta_1 + beta_2 * esmz[i]
              }   
          }
      
      sigma2.g <- inverse(tau.g) * 2 # sigma2.g is genetic covariance
      
      #Priors
      tau.g[1:N, 1:N] ~ dwish(omega.g[,], N) 
      mu[1:N] ~ dmnorm(mean[],precis[,])
      tau.e ~ dgamma(0.001,0.001) 
      tau.s ~ dgamma(0.001,0.001)
      tau.r ~ dgamma(0.001,0.001)
      beta0 ~ dnorm(0, .1)
      beta1 ~ dnorm(0, .1)
      
      }
      

      However, now I am getting the following error message:
      Error in jags.model(model.file, jagsdata, inits, n.chains = 1, quiet = FALSE) :
      RUNTIME ERROR:
      Compilation error on line 4.
      Unable to resolve node mu[1:3]
      This may be due to an undefined ancestor node or a directed cycle in the graph

      although I passed to the JAGS model:

      mu <- c(0,0,0)
      precis <- matrix(c(1, 1/10, 1/10, 1/10, 1, 1/10, 1/10, 1/10, 1), 3,3)
      

      Thank you so much for your help.

       
  • Martyn Plummer

    Martyn Plummer - 2014-08-07

    Surely you want to pass mean there and not mu, i.e.

    mean <- c(0,0,0)
    precis <- matrix(c(1, 1/10, 1/10, 1/10, 1, 1/10, 1/10, 1/10, 1), 3,3)
    
     
    • Inga Schwabe

      Inga Schwabe - 2014-08-08

      Of course.. Sorry to waste your time on such a stupid mistake! I guess I just have seen and adjusted the code too often to see such a simple error.

      Thanks a lot for your time and effort, I really appreciate it.

      Kind regards,
      Inga

       

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