Category Archives: R

Multiple varying coefficients with multiple group-level predictors

Ever tried to set up a multilevel model for e.g. classical educational settings?
I do it with R and JAGS using the rjags package. If you want to use a sort of BUGS language on Mac OS X together with R you have to use JAGS.
My book recommendation for everything about Multilevel Models is  “Data Analysis Using Regression and Multilevel/Hierarchical Models (2007)” by Andrew Gelman.
I was able to set up the BUGS code of the book running in JAGS.
But be carefulI. There is an error in the BUGS code of Chapter 17.2 of the aforementioned book.
The error is in the file “17.2_Bugs_codes.bug” in the section “# Multiple varying coefficients with multiple group-level predictors” on line 131:
B.raw[j,1:K] ~ dnorm (B.raw.hat[j,], Tau.B.raw[,])
should be: 

B.raw[j,1:K] ~ dmnorm (B.raw.hat[j,], Tau.B.raw[,])
as on page 380, line 6 of the book.
A second error is on line 138:
G[k,l] <- xi[k] + G.raw[k,l]
should be: 

G[k,l] <- xi[k] * G.raw[k,l]
The complete corrected code looks like:
# Multiple varying coefficients with multiple group-level predictors
model {
  for (i in 1:n){
    y[i] ~ dnorm(y.hat[i], tau.y)
    y.hat[i] <- inprod(B[classid[i],],X[i,])
  }
  tau.y <- pow(sigma.y, -2)
  sigma.y ~ dunif(0, 100)

  for (k in 1:K){
    for (j in 1:J){
      B[j,k] <- xi[k] * B.raw[j,k]
    }
    xi[k] ~ dunif(0, 100)
  }  
  for (j in 1:J){
    B.raw[j,1:K] ~ dmnorm(B.raw.hat[j,], Tau.B.raw[,]) #this line is erroneous in 17.2_Bugs_codes.bug!
    for (k in 1:K){
      B.raw.hat[j,k] <- inprod(G.raw[k,], U[j,])
    }
  }
  for (k in 1:K){
    for (l in 1:L){
      G[k,l] <- xi[k] * G.raw[k,l] #this line is erroneous in 17.2_Bugs_codes.bug!
      G.raw[k,l] ~ dnorm(0, .0001)
    }
  }

  Tau.B.raw[1:K,1:K] ~ dwish(W[,], df)
  df <- K+1
  Sigma.B.raw[1:K,1:K] <- inverse(Tau.B.raw[,])
  for (k in 1:K){
    for (k.prime in 1:K){
      rho.B[k,k.prime] <- Sigma.B.raw[k,k.prime]/
        sqrt(Sigma.B.raw[k,k]*Sigma.B.raw[k.prime,k.prime])
    }
    sigma.B[k] <- abs(xi[k])*sqrt(Sigma.B.raw[k,k])
  }
}