Nested Random effect in JAGS/ WinBUGS

Question

I am interested in fitting the following nested random effect model in JAGS.

SAS code

proc nlmixed data=data1 qpoints=20;

    parms beta0=2 beta1=1 ;
    bounds vara >=0, varb_a >=0;
     eta = beta0+ beta1*t+ b2+b3;
    p = exp(eta)/(1+exp(eta));
    model TestResult ~ binary(p);
     random b2 ~ normal(0,vara) subject = HHcode;
     random b3 ~ normal(0,varb_a) subject = IDNo_N(HHcode);
    run;

My question: How to specify the random effect part?

I have repeated measurements on individuals. These individuals are further nested in the household. Note: The number of individuals per household vary!

Looking forward to hearing from you

Answer 1

Let's assume that we have two vectors which indicate which house and which individual a data point belongs to (these are things you will need to create, in R you can make these by changing a factor to numeric via as.numeric). So, if we have 10 data points from 2 houses and 5 individuals they would look like this.

house_vec = c(1,1,1,1,1,1,2,2,2,2) # 6 points for house 1, 4 for house 2

ind_vec = c(1,1,2,2,3,3,4,4,5,5) # everyone has two observations

N = 10 # number of data points

So, the above vectors tell us that there are 3 individuals in the first house (because the first 6 elements of house_vec are 1 and the first 6 elements of ind_vec range from 1 to 3) and the second house has 2 individuals (last 4 elements of house_vec are 2 and the last 4 elements of ind_vec are 4 and 5). With these vectors, we can do nested indexing in JAGS to create your random effect structure. Something like this would suffice. These vectors would be supplied in the data.list that you have to include with TestResult

for(i in 1:N){
mu_house[house_vec[i]] ~ dnorm(0, taua)
mu_ind[ind_vec[i]] ~ dnorm(mu_house[house_vec[i]], taub_a)
}

# priors
taua ~ dgamma(0.01, 0.01) # precision
sda <- 1 / sqrt(taua) # derived standard deviation
taub_a ~ dgamma(0.01, 0.01) # precision
sdb_a <- 1 / sqrt(taub_a) # derived standard deviation

You would only need to include mu_ind within the linear predictor, as it is informed by mu_house. So the rest of the model would look like.

for(i in 1:N){
logit(p[i]) <- beta0 + beta1 * t + mu_ind[ind_vec[i]]
TestResult[i] ~ dbern(p[i])
}

You would then need to set priors for beta0 and beta1