dinterval()用于区间截断数据的密度函数？

Question

我是JAGS的新手，正在尝试理解JAGS中用于截断数据的dinterval()函数。我正在对只有每个数据点的上下限（而不是真实值）的粗略数据进行建模。以下是我认为它应该如何工作的简单示例：

每个点的上下限：

> head(lim)
        L        U
[1,] 14.98266 15.68029
[2,] 21.21827 21.91590
[3,] 18.34953 19.04716
[4,] 19.00186 19.69949
[5,] 15.39891 16.09654
[6,] 17.81705 18.51468

编写模型的函数（假设数据来自具有共同均值和方差的正态分布）：

playmodel <- function(){
           for (i in 1:50){
                is.censored[i] ~ dinterval(t[i], lim[i,])
                t[i] ~ dnorm(mu,tau)
               }
           mu ~ dnorm(0,.001)
           tau ~ dgamma(.01,.01)
          } 
          filename <- "toymod.bug"
          write.model(toymod,filename)

对于jags调用，有一些函数和赋值需要注意：

data <- list("lim"=lim)
inits <- list(mu=rnorm(1),tau=rgamma(1,.01,.01),t=as.vector(apply(lim,1,mean)))
#last part is to ensure the starting value is between the upper and lower limit
#each chain will start at the same place for t but this is just for this case
params <- c("mu","tau")

运行此模型：

playmodel.jags <- jags(data,inits, params, model.file="toymod.bug", n.chains=3,
                  n.iter=50000,n.burnin=30000, n.thin=1, DIC=TRUE, 
                  working.directory=NULL,refresh = 50000/50, progress.bar = "text")

当我运行这个程序时会发生什么？

1）我的mu估计值在0附近徘徊，而实际上应该是15。

2）如果DIC=TRUE，它将无法运行：

错误："Error in jags.samples(model, variable.names, n.iter, thin, type = "trace", : Failed to set trace monitor for node deviance"

我相信我做了一些愚蠢的事情，希望有人能帮助我找到正确的方向。

- scottyaz

---

此外，如果我使用 OpenBUGS 中的 I（Lower，Upper）函数对其进行建模，似乎也能正常工作。 - scottyaz

我不确定有多少人从事统计物理学。这可能是少数情况之一，此时在stats.stackexchange.com上进行跨帖可能是个好主意。或者，您可以联系软件包的作者，或查看软件包是否有邮件列表。 - Brandon Bertelsen

Answer 1

以下是Martyn Plummer的回复：

按照您的写法，您的模型没有任何观测结果。您可能已经注意到它非常快地运行。这是因为它从先验中进行前向采样。这就是为什么您的mu的后验均值与先验均值0相同。变量名“is.censored”适用于左截尾或右截尾数据，例如在生存分析中，但不适用于您的问题。因此，我将对其进行重命名为“y”。如果您有

y[j] ~ dinterval(t[j], lim[j,]) 

如果lim [j]有两列，则y [j]可以有三个可能的值。

y[j] = 0 if t[j] <= lim[j,1] 
y[j] = 1 if lim[j,1] < t[j] <= lim[j,2]
y[j] = 2 if lim[j,2] < t[j] 

为了对区间截尾数据进行建模，您需要将y[j]作为模型中的数据提供。在您的情况下，您知道t[j]总是在lim[j,1]和lim[j,2]之间，因此您的数据应该是这样的。

data <- list("lim"=lim, "y"=rep(1,nrow(lim))) 

DIC（Deviance Information Criterion）存在一些问题。因为你的模型没有任何结果数据，所以不能定义偏差。即使您提供了结果数据，您仍然无法获得希望得到的偏差统计数据（包括pD）。偏差将为零，“jags”函数将返回Gelman启发式方法来计算pD（我没有编写这个，所以请不要让我解释它），该值也将为零。实际上需要的似然是

p(lim[j,1] < t[j] <= lim[j,2] | mu, tau) 

但是JAGS为您提供

 p(y[j] | t[j]) 

这里的“focus”是错误的，DIC始终为1。我不知道WinBUGS在这种情况下会做什么。也许它对被审查变量有特殊规定。

- scottyaz