# Cross-classification BUGS

#### Lazar

##### Phineas Packard
Hi All,

I have a hierarchical Bayes model that looks like:

Code:
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:

Code:
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

##### Cookie Scientist
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

##### Phineas Packard
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

##### Phineas Packard
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

##### Phineas Packard
Ok so this works fine:

Code:
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:
Code:
icc.1 <- sigma.a/(sigma.y + sigma.a )
So to extend that to the model above I assume:
Code:
icc.1 <- sigma.a1/(sigma.y + sigma.a1 + sigma.a2)
Would work fine?

#### Jake

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