【问题标题】：Logistic regression in JAGSJAGS 中的逻辑回归
【发布时间】：2016-07-11 19:14:17
【问题描述】：

我是贝叶斯分析的新手。我有一个带有二元响应变量的层次模型。只有一个预测变量（分类），它有 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")

有什么想法吗？谢谢！

【问题讨论】：

标签： bayesian jags

【解决方案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*HLL[j] 产生一个向量（长度为 Ns），igit 无法处理。您可能打算在这一行中用 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 为您提供的错误消息更改为可能为您提供有关编码错误原因的更多线索的内容。

【讨论】：

非常感谢！为什么模型不会以不同的截距收敛？只是为了澄清一下：您那里的模型等效于随机效应，但不是随机截距，对吗？如果我希望截距变化而不是主效应本身，即具有随机截距但没有随机效应的模型怎么办？我已经使用glmer() 运行了此类模型的常客版本，并且该模型确实收敛了。最后，您能否详细说明您对 sigma 的先验的评论？
我的意思是说它可能不会像随机截距那样容易收敛——这取决于你的数据结构，它可能会或可能不会被识别。在任何情况下，通常最好从一个更简单的模型开始，然后逐渐增加复杂性，直到它“崩溃”。我已经编辑了我的答案以澄清这一点。对先验的讨论是一个很大的话题 - 请参阅stat.columbia.edu/~gelman/research/published/taumain.pdf 进行一些讨论（半柯西分布在 runjags 中实现）。总是值得比较不同的先验，看看它如何影响你的结果。
再次感谢，马特。最后一个问题：“加载 GLM 模块”到底是什么意思？我应该传递run.jags() 的额外参数，还是我遗漏了什么？
是的 - 参数应该是 module='glm'

【解决方案2】：

另外（作为上一个答案的扩展），请注意 JAGS 中标准偏差的参数化。 dnorm 由精度参数化，而不是标准偏差。像这样写它更清晰（IMO）：

b ~ dnorm( mu , 1/sigma^2 ) # 注意 1/sigma^2，不是 sigma

sigma ~ dunif(0, 100)

通常将统一先验放在 sigma 上，而不是放在精度上。

[想将其添加为评论而不是答案，但系统不允许。]

【讨论】：