Компиляция RJAGS с категориальной переменной выдает ошибку индекса вне допустимого диапазона

Предпосылки
Попытка моделироватьvolume байкеров на рельсовом пути, что меньше для weekday по сравнению с выходными. RailTrail от mosaicDataсодержит данные, собранные Комиссией по планированию Pioneer Valley по использованию местной железной дороги. За каждые 90 дней записывались рельсовые пути.volume (количество пользователей) и было ли это weekday (ИСТИНА, если да, и ЛОЖЬ в противном случае).

Модель

Yi = объем трейлов (количество пользователей) в день i
Xi = 1 для будних дней, 0 для выходных.

Вероятность

  • Yi ∼ N(mi,s^ ​​2)
  • mi = a + bXi

Приоры

  • а ∼ N(400,100^2)
  • б ∼ 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.

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)
)
Другие вопросы по тегам