0

我绞尽脑汁想办法解决这个问题，但我束手无策！首先，必要的背景：用耙子对湖泊中的水生植物进行采样。你把耙子扔进湖里，然后把它拉回你的船上，然后你就会弄清楚它的尖叉上有什么植物。在我们的例子中，我们测量存在/不存在以及“丰度”，但是以序数/间隔审查的方式->如果物种 X 根本没有在耙子上被注意到，则为 0，如果覆盖 < 25，则为 1耙齿的百分比，如果覆盖范围在 25% 到 75% 之间，则为 2，如果覆盖 > 75%，则为 3。然而，当一个物种的丰度较低时，很容易完全错过它，所以 0 是粗略的——它们可能并不代表真正的缺失，而这正是我们的模型试图探索的问题。

所以，这里真的有三层——一个真正的、完全潜在的丰度，我们根本没有直接观察到，一个部分潜在的“真正的存在/不存在”，因为我们知道真正的存在在哪里，但不知道真正的缺席在哪里，然后我们有我们观察到的存在/不存在数据。更有趣的是，我们认为一些环境变量可能会影响真实丰度和真实发生率，但会有所不同，然后其他变量可能会影响可检测性，而我们正试图梳理这些过程。

所以，无论如何，我的实际模型比我在下面粘贴的要大得多，也更复杂，但这里有一种功能性（但可能在学术上毫无价值）的训练版本，它复制了我得到的错误。

#data setup
N = 1500 #Number of cases
obs = sample(c(0,1,2,3), N, 
             replace=T, prob=c(0.7, 0.2, 0.075, 0.025)) #Our observed, interval-censored data.
X1 = rnorm(N) #Some covariate that probably affects both occurrance and abundance but maybe in different ways.
abundances = rep(NA, times = N) #Abundance is a latent variable we don't directly observe. From elsewhere, I know the values here need to be NAs so the model will know to impute them
occur = rep(1, times = N) #Occurance is a degraded form of our abundance data.

#d will be the initials for the abundance data, since this is apparently needed to jumpstart the imputation.
d = vector()
for(o in 1:N) {
  if (obs[o]==0) { d[o] = 0.025; occur[o] = 0 }
  if (obs[o]==1) { d[o] = 0.15 }
  if (obs[o]==2) { d[o] = 0.5 }
  if (obs[o]==3) { d[o] = 0.875 }
}

#Data
test.data = list("N" = N, 
                 "obs" = obs,
                 "X1" = X1, 
                 "abund" = abundances,
                 "lim" = c(0.05, 0.25, 0.75, 0.9999),
                 "occur" = occur)
#Inits
inits = list(abund = d)

cat("model
    {
    for (i in 1:N) {
    
obs[i] ~ dinterval(abund[i], lim)
abund[i] ~ dbeta(theta[i], rho[i]) T(0.0001, 0.9999)
theta[i] <- mu[i] * epsilon
rho[i] <- epsilon * (1-mu[i])

logit(mu[i]) <- alpha1 + X.beta1 * X1[i]

occur[i] ~ dbern(phi[i])
logit(phi[i]) <- alpha2 + X.beta2 * X1[i]

    }
    
#Priors
epsilon ~ dnorm(5, 0.1) T(0.01, 10)
alpha1 ~ dnorm(0, 0.01)
X.beta1 ~ dnorm(0, 0.01)
alpha2 ~ dnorm(0, 0.01)
X.beta2 ~ dnorm(0, 0.01)
}
 ", file = "training.txt")

test.run = jags.model(file = "training.txt", inits = inits, data=test.data, n.chains = 3) 

params = c("epsilon", 
                "alpha1",
                "alpha2",
                "X.beta1", 
                "X.beta2")

run1 = run.jags("training.txt", data = test.data, n.chains=3, burnin = 1000, sample = 5000, adapt = 4000, thin = 2,
                monitor = c(params), method="parallel", modules = 'glm')

最后，我得到了这个错误，每当我尝试做这样的远程操作时，我总是会得到这个错误：

图信息：观察到的随机节点：3000 未观察到的随机节点：1505 总图大小：19519。读取参数文件 inits1.txt。节点obs 1中的初始化模型 错误节点 与父节点不一致

我已经阅读了所有我能找到的包含此错误的帖子，包括this、this、this和this。我可以从我的研究和测试中推测出该错误可能是由于以下原因之一而发生的。

  1. 我对潜在丰度变量的缩写不知何故是不够的。听起来这需要非常有用的初始值才能工作。
  2. 我的一个或多个先验是允许不允许的值或它们太宽泛，这会以某种方式造成问题。这可能是一个特别严重的问题，因为我使用的 beta 分布对没有 0 和 1 之外的值有很强的要求。
  3. 我错误地使用了 dinterval() 函数，这似乎是因为包含它的行总是会触发错误。
  4. 我的模型不知何故指定错误。

但我看不出哪里出错了——我已经为 1 和 2 尝试了许多不同的选项，据我从文档中可以看出（参见第 55-56 页），我正在正确使用 dinterval . 我错过了什么？？

如果它是相关的，从我收集的内容来看， dinterval() 的想法是 ~ 左侧的变量是第一个参数中给出的变量的区间删失版本（这里是丰度）。然后，第二个参数（此处为 lim）是一个“断点”向量，它指示丰度数据最终位于哪个区间。因此，在这里，如果您低于最低 lim，则最终观察到的丰度代码为 0 （此处为 0.05），如果您在 lim 的前两个值之间，则为 1 等。这就像丰度变量被推过由 lim 变量创建的“分箱筛”以产生分箱输出变量，我们观察到的丰度.

任何指导将是最受欢迎的！

4

1 回答 1

1

我已经使用 JAGS 4.3.0 和 rjags 4-12 运行了您的示例。对我来说，带有 rjags 的版本运行正常。带有 runjags 的版本不起作用，因为您没有提供初始值。这很容易通过添加参数来解决

inits=list(inits, inits, inits)

调用 run.jags()。

您已经正确理解了 dinterval 的用途。这是一个“可观察函数”，它通过可能性对其参数施加约束。使用 dinterval 时，您必须始终提供满足第一次迭代约束的初始值。据我所见，您的初始值确实满足约束，这通过我可以运行您的示例（使用初始值）这一事实得到验证。

于 2022-02-21T10:45:49.587 回答