JAGS 中的逻辑回归

Question

我是贝叶斯分析的新手。我有一个带有二元响应变量的层次模型。只有一个预测变量（分类），它有 3 个级别：HLL、LHL 和 LLL。我通过对这些级别 all 进行虚拟编码来准备我的数据文件。我的机型规格如下：

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*HLL[j] 产生一个向量（长度为 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 由精度参数化，而不是标准偏差。这样写更清楚（IMO）：

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

西格玛 ~ dunif( 0 , 100 )

通常统一先验放在西格玛上，而不是精度上。

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

JAGS 中的逻辑回归

Logistic regression in JAGS

bayesian

jags

错误

仅拦截模型工作正常