OpenBUGS: Hierarchical Model unable to compile

Question

I've implement a model in OpenBUGs from the following data:

The data consists of 10 independent studies - 5 trials of 2 different medications (Med = 0 or 1) then a total of infected participants for that trial (infected) - and a Total of participants in the trial.

bugsData(list(y = df$Infected, n = 10, J = 5 , t = df$Total, x = df$Med) , file="Bugs_Data.txt")

My original model can written as:

Y[i]|mu[i],t[i] ~ Bin(mu[i], t[i]), i = 1,.....10
logit(mu[i]) = Beta0 + Beta1*x[i]

Which I'm able to implement in OpenBUGS like so:

model {
## LIKELIHOOD
  for (i in 1:n) {
    y[i] ~ dbin(mu[i], t[i])
    logit(mu[i]) <- Beta0 + Beta1*(x[i])
  }
  ## NORMAL PRIORS 
  Beta0 ~ dnorm(0, 0.0001)
  Beta1 ~ dnorm(0, 0.0001)
}

However I then would like to implement a hierarchical model, in which the s(i) denotes the trial from column trial where observation i was taken:

Y[i]|mu[i],t[i] ~ Bin(mu[i], t[i]), i = 1,.....10
logit(mu[i]) = Beta0,s(i) + Beta1*x[i]
Beta0,j ~ Normal(mu_bo, sigma2_b0), j = 1,....5

I've attempted a model to the best of my ability with no success as of yet.

model {
## LIKELIHOOD
  for (j in 1:J) {
    for (i in 1:n) {
      y[i] ~ dbin(mu[i], t[i])
      logit(mu[i]) <- Beta0[j] + Beta1*(x[i])
    Beta0[j] ~ dnorm(mu_b0, tau_b0)
    }
  }
  
  ## PRIORS 
  mu_b0 ~ dnorm(0, 0.0001)
  tau_b0 ~ dgamma(0.0001, 0.0001)
  sigma2_b0 <- 1 / tau_b0
  Beta1 ~ dnorm(0, 0.0001)
}

The current model isn't able to compile with my data as is.

Answer 1

You need to supply a vector of length n which indexes what group data point i belongs to.

model {
## LIKELIHOOD
for (i in 1:n) {
  y[i] ~ dbin(mu[i], t[i])
  logit(mu[i]) <- Beta0[group_vec[i]] + Beta1*(x[i])
}
## PRIORS 
mu_b0 ~ dnorm(0, 0.0001)
tau_b0 ~ dgamma(0.0001, 0.0001)
sigma2_b0 <- 1 / tau_b0
Beta1 ~ dnorm(0, 0.0001)
for(j in 1:J){
  Beta0[j] ~ dnorm(mu_b0, tau_b0)
}
}

group_vec is of length n, and it takes an integer of which group data point i belongs to. So, for example, if we had 21 datapoints and 3 groups.

groups <- factor(rep(c("A", "B", "C"), each = 7))

then we could construct group_vec as

group_vec <- as.numeric(groups)
[1] 1 1 1 1 1 1 1 2 2 2 2 2 2 2 3 3 3 3 3 3 3

This is called nested indexing, and it's a really common way to specify hierarchical models, you just supply it along with your data. In your case, you already have trial as an integer. So you could either make group_vec be the trial column in your data, or rename group_vec in the model to be Trial, whichever you prefer.