使用分类变量进行 RJAGS 编译会引发索引超出范围错误答案

【问题标题】：RJAGS compilation with categorical variable throws index out of range error使用分类变量进行 RJAGS 编译会引发索引超出范围错误
【发布时间】：2020-06-16 21:57:07
【问题描述】：

背景
试图在铁路小道上模拟volume 的骑自行车的人，与周末相比，weekday 更少。来自mosaicData 的RailTrail 包含先锋谷规划委员会收集的有关当地铁路使用情况的数据。对于 90 天中的每一天，他们记录了轨道轨迹 volume（用户数量）以及它是否是 weekday（如果是，则为 TRUE，否则为 FALSE）。

型号

Yi = 第 i 天的跟踪量（用户数）
Xi = 1 表示工作日，0 表示周末。

可能性

易∼N(mi,s^2)
mi =a+bXi

先验

a ∼ N(400,100^2)
b ∼ N(0,200^2)
s ∼ Unif(0,200)

代码

尝试在 R 中按如下方式实现：

library(rjags)
library(mosaicData)

data(RailTrail)

# DEFINE the model    
rail_model_1 <- "model{
    # Likelihood model for Y[i]
    for(i in 1:length(Y)) {
      Y[i] ~ dnorm(m[i], s^(-2))
      m[i] <- a + b[X[i]]
    }

    # Prior models for a, b, s
    a ~ dnorm(400, 100^(-2))
    b[1] <- 0
    b[2] ~ dnorm(0, 200^(-2))
    s ~ dunif(0, 200)
}"

尝试使用以下代码编译上述模型：

# COMPILE the model
rail_jags_1 <- jags.model(
  textConnection(rail_model_1),
  data = list(Y = RailTrail$volume, X = RailTrail$weekday),
  inits = list(.RNG.name = "base::Wichmann-Hill", .RNG.seed = 10)
)

错误

但是，我在尝试编译时遇到以下错误：

Error in jags.model(textConnection(rail_model_1), data = list(Y = RailTrail$volume,  : 
  RUNTIME ERROR:
Compilation error on line 5.
Index out of range taking subset of  b

问题

你能帮我解决这里的问题吗？我在 Ubuntu 20.04、MacOS Catalina 以及 RStudio Cloud 中对此进行了测试——同样的错误。 rjags.version() 是 4.3.0。

【问题讨论】：

我不了解您对b 参数所做的事情；你想要m[i] <- a + b*X[i] 并将先验定义为b ~ dnorm(0, 200^(-2))
a = 典型周末成交量 a + b = 典型工作日成交量 b 是工作日和周末之间的差异。 b 有 2 个级别：b[1]、b[2]。 b[1]手动设置为0，b[2] ~ dnorm(0, 200^(-2))
或者你可以单独为它们建模（我认为这是你所做的事情）；将 m[i] <- b[X[i]] 与先前的 ` for(i in 1:2) {b[i] ~ dnorm(0, 200^(-2))}` 一起使用。对于第二个符号，您将需要 X = factor(RailTrail$weekday)) in jags.model
...你评论。如果oyu包含拦截a，那么b将只有一个参数，所以你只需要在b上设置先验：b ~ ...
实际上，我只是用X = factor(RailTrail$weekday)) 运行了你的代码，它会执行，但你会注意到它是过度参数化的——如果你包含拦截，你不需要b[1] <- 0

标签： r bayesian categorical-data rjags bayesglm

【解决方案1】：

@user20650分享：

代码在编译语句中使用显式X = factor(RailTrail$weekday))，即，

# COMPILE the model
rail_jags_1 <- jags.model(
  textConnection(rail_model_1),
  data = list(Y = RailTrail$volume, X = factor(RailTrail$weekday)),
  inits = list(.RNG.name = "base::Wichmann-Hill", .RNG.seed = 10)
)

【讨论】：