JAGS中的Logistic回归

Question

我是贝叶斯分析的新手。我有一个具有二进制响应变量的层次模型。只有一个预测因子(范畴)，它有3个等级: HLL、LHL和LLL。我准备了我的数据文件通过虚拟编码所有这些级别。我的模型规范如下：

cat("

model{
  for(i in 1:Ny){
    y[i] ~ dbern(mu[s[i]])
  }
for(j in 1:Ns){
mu[j] <- ilogit(b0[j] + b1*HLL[j] + b2*LHL[j] + b3*LLL[j])

b0[j] ~ dnorm(mu0, sigma0)
b1[j] ~ dnorm(mu1, sigma1)
b2[j] ~ dnorm(mu2, sigma2)
b3[j] ~ dnorm(mu3, sigma3)

  }

mu0 ~ dnorm(0, 0.001)
sigma0 ~ dunif(0, 100)

mu1 ~ dnorm(0, 0.001)
sigma1 ~ dunif(0, 100)

mu2 ~ dnorm(0, 0.001)
sigma2 ~ dunif(0, 100)

mu3 ~ dnorm(0, 0.001)
sigma3 ~ dunif(0, 100)
}
", fill = TRUE, file = "generalModel.txt")

错误

基本上，我想得到HLL和LHL的估计值(使用LLL作为参考。职等)。这个模型不能运行，我也不知道为什么。以下是错误消息：

Calling 3 simulations using the parallel method...
Following the progress of chain 1 (the program will wait for all chains to finish
before continuing):
Welcome to JAGS 4.2.0 on Sun Jul 10 00:10:00 2016
JAGS is free software and comes with ABSOLUTELY NO WARRANTY
Loading module: basemod: ok
Loading module: bugs: ok
. . Reading data file data.txt
. Compiling model graph
   Resolving undeclared variables
   Allocating nodes
RUNTIME ERROR:
Invalid vector argument to ilogit
Deleting model
. Reading parameter file inits1.txt
Can't set RNG name. No model!
Can't set initial values. No model!
. Can't initialize. No model!
. Adaptation skipped: model is not in adaptive mode.
. Updating 1000
-------------------------------------------------| 1000
Can't update. No model!

. Can't set monitor. No model!
. Can't set monitor. No model!
. Can't set monitor. No model!
. Can't set monitor. No model!
. Updating 10000
-------------------------------------------------| 10000
Can't update. No model!

. No model
. Can't dump CODA output. No model!
. Can't dump samplers. No model!
. Updating 0
Can't update. No model!
Can't update. No model!
. Deleting model
. 
All chains have finished
Note: the model did not require adaptation
Error in runjags.readin(directory = startinfo$directory, silent.jags = silent.jags,  : 
  All the simulations appear to have crashed - check the model output in failed.jags() for clues
In addition: Warning messages:
1: No initial values were provided - JAGS will use the same initial values for all chains 
2: You attempted to start parallel chains without setting different PRNG for each chain, which is not recommended.  Different .RNG.name values have been added to each set of initial values. 
Note: Either one or more simulation(s) failed, or there was an error in processing
the results.  You may be able to retrieve any successful simulations using:
results.jags("/private/var/folders/nv/gznh75k93cv1wp35q1hvkkg00000gn/T/RtmpYRkQYd/runjagsfiles7c8d79109b99",
recover.chains=TRUE)
See the help file for that function for possible options.
To remove failed simulation folders use cleanup.jags() - this will be run
automatically when the runjags package is unloaded

纯拦截模型运行良好

我成功地运行了仅拦截的模型，这些模型与上面的相同。在本例中，我运行了一个模型，例如，只使用LLL，而另一个模型只使用HLL。然后，我绘制了两个后验的差值，结果与glmer()模型中的HLL估计值基本一致，其中LLL是参考值。关卡。

cat("

model{
  for(i in 1:Ny){
    y[i] ~ dbern(mu[s[i]])
  }
  for(j in 1:Ns){
    mu[j] <- ilogit(b0[j])
    b0[j] ~ dnorm(mu0, sigma)
  }

  mu0 ~ dnorm(0, 0.001)
  sigma ~ dunif(0, 100)
}
", fill = TRUE, file = "modelA.txt")

有什么想法吗？谢谢！

Answer 1

输出的重要部分是这个错误：

RUNTIME ERROR:
Invalid vector argument to ilogit

这来自您的模型的以下部分：

mu[j] <- ilogit(b0[j] + b1*HLL[j] + b2*LHL[j] + b3*LLL[j])

参数b1、b2和b3由j在下面的循环中索引，因此b1*HLLj生成一个向量(长度为Ns)，ilogit无法处理该向量。您可能打算在这一行中用j索引b1等，即：

mu[j] <- ilogit(b0 + b1[j]*HLL[j] + b2[j]*LHL[j] + b3[j]*LLL[j])

请注意，我已经删除了b0上的索引，以便有固定的截获，否则模型可能不会那么快地收敛，但是一旦模型的其余部分能够正常工作，您可以在以后再添加这个索引。这和你的意图之间可能还有其他的区别，但是如果没有更多关于你想要做的事情的信息，我很难说出来。例如，您可能不想用j索引b0、b1、b2和b3，而是希望在循环之外定义单个b0、b1、b2和b3 --但这并不是分层的，您说这是您想要的。

还有几点：

( 1)我会非常小心处理以下类型的前科：

sigma0 ~ dunif(0, 100)

对于sigma0来说，这可能是一个信息丰富的先导，并且将倾向于将参数拉向更高的值。

2)确保加载GLM模块--这将允许块更新参数，从而提高收敛性。

( 3)这是一个较次要的问题，但更习惯于这样写：

logit(mu[j]) <- ....

超过：

mu[j] <- ilogit(...)

除了让这一点更清楚地说明您正在编写GLM之外，它还将更改JAGS提供给您的错误消息，从而为您提供更多关于编码错误原因的线索。

Answer 2

此外(在扩展前面的答案)，注意参数化的标准偏差在JAGS。dnorm是由精度而不是标准差来参数化的。(国际海事组织)用这样的方式写它更清楚：

B~d范数( mu，1/σ^2)#注意1/西格玛^2，而不是西格玛

sigma ~ dunif( 0,100 )

通常，均匀的先验是在西格玛，而不是在精度上。

我想把它作为评论而不是回答，但是系统不允许我这样做。