如何覆盖被观察节点的值？

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
}

- Paco Cruz

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

- Paco Cruz