How to specify Bayesian mixed effects model in BUGS

Question

I posted this earlier in the week then retracted the question when I found a good source, not wanting to waste people's time. I haven't made much progress I'm afraid. In trying to be a good citizen here I will make the problem as clear as possible. I suspect there will be few takers.

I have a dataframe in R I wish to analyse in BUGS or R. It is in long format. It consists of multiple observations on 120 individuals, with a total of 885 rows. I am examining the occurrence of a categorical outcome - but that's not really relevant here. The question is about something deeper.

The model I have been using up to here is

  mymodel<-gee(Category ~ Predictor 1 + Predictor 2..family=binomial(link="logit"),
  data=mydata, 
   id=Person)

with a marginal model essentially accounting for the clustering of patients. I then examined

 mymodel<-gee(Category ~ Predictor 1 + Predictor 2.. , family=binomial(link="logit"), 
  corstr = "AR-M", 
  data=mydata, id=Person)

in order to account for the time ordering of the observations on the individual people.

This didn't change much.

Then I tried to model them using the following set of MCMCPack commands:

 mymodel<-MCMCglmm(category~  Predictor1 + Predictor2..,  
 data=mydata, family=binomial(link="logit"))

An examination of the output was thrilling, showing statistical significance for many predictors. I hailed myself as a newly converted bayesian, until I realised I hadn't accounted for repeated measures within patients.

I understand that I have to account for that. I understand that this may mean fitting a hyperprior for each individual - is that right? What form will this take in BUGS?

Here's a basic log reg model: (kudos to Kruschke, J., Indiana)

    model {
  for( i in 1 : nData ) {
    y[i] ~ dbern( mu[i] )
    mu[i] <- 1/(1+exp(-( b0 + inprod( b[] , x[i,] ))))
  }
  b0 ~ dnorm( 0 , 1.0E-12 )
   for ( j in 1 : nPredictors ) {
    b[j] ~ dnorm( 0 , 1.0E-12 )
  }
}

However, no hyperprior here for the individual. Here's my best attempt so far at a within-individual design, accounting for repeated measures within people:

Here's Jackman's model for JAGS

1 model{
2 ## loop over data for likelihood
3 for(i in 1:n){
4  y[i] ~ dbern( mu[i] )
    mu[i] <- 1/(1+exp(-( b0 + inprod( b[] , x[i,] ))))
6 }
7 sigma ˜ dunif(0,20) ## prior on standard deviation
8 tau <- pow(sigma,-2) ## convert to precision
9
10 ## hierarchical model for each state’s intercept & slope
11 for(p in 1:50){
12 beta[p,1:2] ˜ dmnorm(mu[1:2],Tau[,]) ## bivariate normal
13 }
14
15 ## means, hyper-parameters
16 for(q in 1:2){
17 mu[q] ˜ dnorm(0,.0016)

}

Here's my bastard-child model for BUGS

1 model{
2 ## loop over data for likelihood
3 for(i in 1:n){
4 mu.y[i] <- alpha + beta[s[i],1] + beta[s[i],2]*(j[i]-jbar)
5 demVote[i] ˜ dnorm(mu.y[i],tau)
6 }
7 sigma ˜ dunif(0,20) ## prior on standard deviation
8 tau <- pow(sigma,-2) ## convert to precision
9
10 ## hierarchical model for each state’s intercept & slope
11 for(p in 1:120){
12 beta[p,1:2] ˜ dmnorm(mu[1:2],Tau[,]) ## bivariate normal
13 }
14
15 ## means, hyper-parameters
16 for(q in 1:2){
17 mu[q] ˜ dnorm(0,.0016)
  }

Can somebody let me know if I'm heading in the right direction. My understanding of this is growing, but slowly. Please be gentle. I'm a medic, not a statistic! I have used R quite a bit, but I'm new to BUGS, and new to Bayes.

Thanks,

R

Answer 1

You are (were) almost there. Just a few comments - you don't have to make the prior for the beta[,1:2] parameters a joint MV normal; you can make the prior such that beta[i,1] and beta[i,2] are independent, which simplifies things (for example, no prior covariance need be specified.) Note that doing so doesn't mean they will be independent in the posterior.

Other comments: Since you have a constant term - alpha - in the regression, the components beta[,1] should have zero mean in the prior. Also, you don't have a prior for alpha in the code.

Here's a model with hierarchical intercept and slope terms; I've tried to keep to your priors and notation where possible, given the changes:

model {
  for(i in 1:n){
    mu.y[i] <- alpha + beta0[s[i]] + beta1[s[i]]*(j[i]-jbar)
    demVote[i] ~ dnorm(mu.y[i],tau)
  }

  alpha ~ dnorm(0, 0.001) ## prior on alpha; parameters just made up for illustration
  sigma ~ dunif(0,20) ## prior on standard deviation
  tau <- pow(sigma,-2) ## convert to precision

  ## hierarchical model for each state’s intercept & slope
  for (p in 1:120) {
     beta0[p] ~ dnorm(0, tau0)
     beta1[p] ~ dnorm(mu1, tau1)
  }

  ## Priors on hierarchical components; parameters just made up for illustration
  mu1 ~ dnorm(0, 0.001) 
  sigma0 ~ dunif(0,20)
  sigma1 ~ dunif(0,20)
  tau0 <- pow(sigma0,-2)
  tau1 <- pow(sigma1,-2)
}

A very useful resource for hierarchical models, including some "tricks" to speed up convergence, is Gelman and Hill.

(A little late with the answer, but may be helpful to some future questioner.)