如何覆盖观察节点的值？

Question

是JAGS的新手，并使用this paper中引用的模型。看起来作者使用的JAGS的早期版本更接近BUGS，因为在某些时候，它出现在模型块中（第28-29行）：

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

JAGS抱怨以下错误：

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

检查JAGS手册，错误很明显。在A.4数据转换部分中，声明了BUGS允许在关系的左侧两次放置一个随机节点。作为解决方法，它提供了将数据转换包含在单独的data块中的功能。它仍然失败。

有人尝试复制这项工作并获得成功吗？有任何提示吗？

---

下面的完整JAGS模型。可疑分配是z[i,K] <- 0和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
}

Answer 1

我与论文作者取得了联系，并承认JAGS的后发行破坏了原始代码。用他的话说，解决方案是...

“但是有一个相对简单的解决方法，其中涉及使用一个额外的变量来有效表示潜在响应中的K-1差异，然后将其映射到Z变量。我将附加一个包含更新脚本的文件确实需要更改为算法创建初始值的方式。makeinits函数（在data.R文件中）也已更新。“

下面提供了新的模型代码：

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
}

具有新的makeinits功能：

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)]
}