How to overwrite value of observed node?

Question

New to JAGS and working with a model referred in this paper. Looks like the authors used a previous version of JAGS that was closer to BUGS since at some point, this appears in the model block (lines 28-29):

  z[i,1:(K-1)] ~ dmnorm(etas[i,], omega)
  z[i,K] <- 0

JAGS complains with the following error:

Error in setParameters(init.values[[i]], i) : Error in node z[1,4]
Cannot overwrite value of observed node

Checking the JAGS manual, the error is evident. In section A.4 Data Transformations, it states that BUGS allows to put a stochastic node twice on the left hand side of a relation. As a workaround, it offers to include the data transformation on a separate data block. It still fails.

Has anybody tried to replicate this work and found success? Any hints?

Full JAGS model below. The suspect assignments are z[i,K] <- 0 and beta[j,K] <- 0

data {
   for (j in 1:(K-1)) {
      for (i in 1:N) {
         ones[i,j] <- 1
      }
      for (k in 1:K) {
         C[j,k] <- equals(j,k) - equals(k,K)
      }   
   }
   R <- (K-1) * C %*% t(C)
   lower <- -1e+3
   upper <-  1e+3
}
model { 
   for (i in 1:N) {
      for (k in 1:(K-1)) {
         bounds[i,k,1] <- equals(y[i,k],K)*lower + inprod(z[i,], equals(y[i,],y[i,k] + 1))
         bounds[i,k,2] <- equals(y[i,k],1)*upper + inprod(z[i,], equals(y[i,],y[i,k] - 1))
         ones[i,k] ~ dinterval(z[i,k], bounds[i,k,])
         etas[i,k] <- inprod(x[i,], beta[,k])
      }
      z[i,1:(K-1)] ~ dmnorm(etas[i,], omega)
      z[i,K] <- 0
   }
   for (j in 1:(P+1)) {
      for (k in 1:(K-1)) {
         beta[j,k] ~ dnorm(0,1e-3)
      }
      beta[j,K] <- 0
      for (k in 1:K) {   
          for (t in 1:K) {
             beta.identified[j,k,t] <- (beta[j,k] - beta[j,t])*sqrt(const)
          }
      }
   }
   omega ~ dwish(R, K-1)
   sigma <- inverse(omega)
   const <- pow(K/exp(logdet(sigma)), 1/K)
   sigma.identified <- sigma*const
}

score 3 · Answer 1 · answered Nov 13 '19 at 03:31

I reached out to the paper author and acknowledged that a posterior release of JAGS broke the original code. In his words, the solution is...

"But there is a relatively simple work-around that involves using an extra variable that effectively represents the K-1 differences in the latent responses and then mapping that to the Z variable. I'll attach a file that contains the updated scripts from that paper. This does require a change in the way that initial values are created for the algorithm. The makeinits function (in the data.R file) has been updated as well."

The new model code is provided below:

data {
   for (j in 1:(K-1)) {
      for (i in 1:N) {
         ones[i,j] <- 1
      }
      for (k in 1:K) {
         C[j,k] <- equals(j,k) - equals(k,K)
      }   
   }
   R <- (K-1) * C %*% t(C)
   lower <- -1e+3
   upper <-  1e+3
}
model { 
   for (i in 1:N) {
      for (k in 1:(K-1)) {
         bounds[i,k,1] <- equals(y[i,k],K)*lower + inprod(z[i,], equals(y[i,],y[i,k] + 1))
         bounds[i,k,2] <- equals(y[i,k],1)*upper + inprod(z[i,], equals(y[i,],y[i,k] - 1))
         ones[i,k] ~ dinterval(z[i,k], bounds[i,k,])
         etas[i,k] <- inprod(x[i,], beta[,k])
      }
      u[i,1:(K-1)] ~ dmnorm(etas[i,], omega)
      z[i,1:(K-1)] <- u[i,1:(K-1)]
      z[i,K] <- 0
   }
   for (j in 1:(P+1)) {
      for (k in 1:(K-1)) {
         beta[j,k] ~ dnorm(0,1e-3)
      }
      beta[j,K] <- 0
      for (k in 1:K) {   
          for (t in 1:K) {
             beta.identified[j,k,t] <- (beta[j,k] - beta[j,t])*sqrt(const)
          }
      }
   }
   omega ~ dwish(R, K-1)
   sigma <- inverse(omega)
   const <- pow(K/exp(logdet(sigma)), 1/K)
   sigma.identified <- sigma*const
}

with new makeinits function:

makeinits <- function (y) {
    n <- nrow(y)
    m <- ncol(y)
    z <- matrix(NA, n, m)
    for (i in 1:n) { 
        z[i,] <- sort(rnorm(m))[rank(-y[i,])]
        z[i,] <- z[i,] - z[i,m]
        }
    z[,-ncol(y)]
}

How to overwrite value of observed node?

1 Answers1