2

JAGS 的新手并使用本文中提到的模型。看起来作者使用了更接近 BUGS 的早期版本的 JAGS，因为在某些时候，这出现在模型块中（第 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
}
4

1 回答 1

2

我联系了论文作者并承认 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)]
}
于 2019-11-13T03:31:53.850 回答