1

Background
Trying to model volume of bikers in a rail trail which is less for a weekday as compared to a weekend. RailTrail from mosaicData contains data collected by the Pioneer Valley Planning Commission on the usage of a local rail-trail. For each of 90 days, they recorded the rail-trail volume (number of users) and whether it was a weekday (TRUE if yes and FALSE otherwise).

Model

Yi = trail volume (# of users) on day i
Xi = 1 for weekdays, 0 for weekends.

Likelihood

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

Priors

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

Code

Trying to implement this in R as follows:

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)
}"

Attempting to compile the model above using the following code:

# 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

However, I am getting the following error in my attempt to compile:

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

Question

Can you please help me with what is wrong here? I tested this in Ubuntu 20.04, MacOS Catalina as well as RStudio Cloud - same error. rjags.version() is 4.3.0.

hnovice
  • 43
  • 5
  • I don't follow what you are doing with the `b` parameter; do you want `m[i] <- a + b*X[i]` and to define the prior as `b ~ dnorm(0, 200^(-2))` – user20650 Jun 16 '20 at 22:07
  • `a` = typical weekend volume `a + b` = typical weekday volume `b` is the difference between weekday and weekend. b has 2 levels: `b[1]`, `b[2]`. `b[1]` is manually set to 0, and `b[2] ~ dnorm(0, 200^(-2))` – hnovice Jun 16 '20 at 22:14
  • or you could model them separately (which is knid of what you have done, I thnk); use `m[i] <- b[X[i]]` with prior ` for(i in 1:2) {b[i] ~ dnorm(0, 200^(-2))}` . For this second notation you will need `X = factor(RailTrail$weekday))` in `jags.model` – user20650 Jun 16 '20 at 22:17
  • ... re you comment. If oyu are including the intercept `a`, then `b` will have only one parameter so you only need to set the prior on `b`: `b ~ ... ` – user20650 Jun 16 '20 at 22:19
  • actually, I just ran you code but with `X = factor(RailTrail$weekday))` and it executes but you will notice that it is over-paramterised -- yo udont need `b[1] <- 0` if you include the intercept – user20650 Jun 16 '20 at 22:27
  • interesting. I checked that `RailTrail$weekday` explicitly before so did not think of including it in the model statement explicitly. Thanks for the tip! Works now. – hnovice Jun 16 '20 at 23:49

1 Answers1

3

as shared by @user20650:

The code works with using explicit X = factor(RailTrail$weekday)) in the compile statement, i.e., 

# 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)
)
hnovice
  • 43
  • 5