Cross-classification BUGS

[Lazar]

Hi All,

I have a hierarchical Bayes model that looks like:

model{
  #Overaching model
  for (i in 1:n){
    y[i] ~ dnorm (yHat[i], tau.y)
    yHat[i] <- a[id[i]] + b.group*group[i] + b.wave*wave[i] + b.int*wave[i]*group[i]
  }
  #Fixed Effects
  b.group ~ dnorm(0, .0001)
  b.wave ~ dnorm(0, .0001)
  b.int ~ dnorm(0, .0001)
  
  #For missing y
  tau.y <- pow(sigma.y, -2)
  sigma.y ~ dunif (0, 100)
  
  #Loop through id
  for (j in 1:S2) {
      a[j] ~dnorm(mu.a, tau.a)
  }
  #Random Effects
  mu.a ~ dnorm(0, .0001)
  tau.a <- pow(sigma.a, -2)
  sigma.a ~ dunif (0, 100)

  #Effect of group at particular points in time
  b1 <- b.group + b.int*1
  b4 <- b.group + b.int*4
  b12 <- b.group + b.int*12
  b24 <- b.group + b.int*24

}

The model is basically a repeated measures RCT with observations at baseline, week 1,4,12, and 24. So the data is nested under participant.

I also have another nesting called tid. What I want to do is to be able to alter the code above to account for this additional nesting. Something like the lme4 type model:

lmer(y ~ group*wave + (1|id) + (1|tid), data=temp)

Any help on how to modify the bugs/JAGS script above to be able to do this would be greatly appreciated.

 

[Jake]

I'm not really familiar with BUGS syntax but can't you just add another intercept term to yHat (e.g., so there's an a.id and an a.tid) and another loop to set the distribution of the a.tids?

 

[Lazar]

I tried that but got nowhere BUT that could be because of a stray curly bracket I just found. Working on fixing that now.

 

[Lazar]

All works. The problem was a funny one. The ids I had were not consecutive numbers and so when I looped over them I kept getting subscript out of bounds. Did as.numeric(factor(x)) and all worked fine.

 

[Lazar]

Ok so this works fine:

model{
for (i in 1:n){
y[i] ~ dnorm (yHat[i], tau.y)
yHat[i] <- a1[id[i]] + a2[tid[i]] + b.group*group[i] + b.wave*wave[i] + b.int*wave[i]*group[i] +
           b.site*sid[i]
}

b.group ~ dnorm(0, .0001)
b.wave ~ dnorm(0, .0001)
b.int ~ dnorm(0, .0001)
b.site ~ dnorm(0, .0001)

tau.y <- pow(sigma.y, -2)
sigma.y ~ dunif (0, 100)

for (j in 1:S2) {
  a1[j] ~dnorm(mu.a1, tau.a1)
}

mu.a1 ~ dnorm(0, .0001)
tau.a1 <- pow(sigma.a1, -2)
sigma.a1 ~ dunif (0, 100)

for (k in 1:S1) {
  a2[k] ~dnorm(mu.a2, tau.a2)
}

mu.a2 ~ dnorm(0, .0001)
tau.a2 <- pow(sigma.a2, -2)
sigma.a2 ~ dunif (0, 100)
  
a <- mu.a1 + mu.a2
b1 <- b.group + b.int*1
b4 <- b.group + b.int*4
b12 <- b.group + b.int*12
b24 <- b.group + b.int*24

}

I would now like to estimate ICCs in this model directly. For a simple two level model (see op) I assume the code would be:

icc.1 <- sigma.a/(sigma.y + sigma.a )

So to extend that to the model above I assume:

icc.1 <- sigma.a1/(sigma.y + sigma.a1 + sigma.a2)

Would work fine?

 

[Jake]

Yes, that would be one of the ICCs one could compute from this model.