Harold Doran writes,

I know that you have a new book coming out on multilevel modeling and I see that you have some code using lmer on your website. I don’t know if it is helpful or not, but below is some code I wrote (basically tweaking the code in your book) to learn how to estimate a multilevel model using the gibbs sampler. This is equivalent to the following model in lmer

library(mlmRev)
fm <- lmer(math ~ 1 +(1|sch), star) It is nothing too fancy, but it was helpful for me as a practice set and maybe it might be for others. I am often desirous for a "Rosetta Stone" so I can learn to move between different ways for estimating models.

Here’s the code:

theta.update <- function() { theta.hat <- (mu/tau^2 + (n_j/sigma.y^2)*y_j)/(1/tau^2 + n_j/sigma.y^2) V.theta <- 1/(1/tau^2 + n_j/sigma.y^2) rnorm(J, theta.hat, sqrt(V.theta)) } mu.update <- function(){ rnorm(1,mean(theta), tau/sqrt(J)) } tau.update <- function(){ sqrt(sum((theta-mu)^2)/rchisq(1, J-1)) } sigma.update <- function(){ n * (1/n * (sum((star$math - rep(theta.update(), n_j))^2, na.rm=T))) / rchisq(1,n) } n.chains <- 5 n.iter <- 1000 J <- length(levels(star$sch)) # number of schools n_j <- with(star, tapply(id, sch, length)) # number of students within each school n <- sum(n_j) # total number of students sims <- array(NA, c(n.iter, n.chains, J+3)) dimnames(sims) <- list(NULL, NULL, c(paste ('theta[', 1:length(levels(star$sch)), ']', sep=''), 'mu', 'tau', 'sigma')) y_j <- as.numeric(with(star, tapply(math, sch, mean, na.rm=T))) sigma.y <- sd(y_j) for(m in 1:n.chains){ mu <- rnorm(1, mean(y_j), sd(y_j)) tau <- runif(1,0, sd(y_j)) for(t in 1:n.iter){ theta <- theta.update() mu <- mu.update() sigma2 <- sqrt(sigma.update()) tau <- tau.update() sims[t,m,] <- c(theta, mu, tau, sigma2) } } R2WinBUGS:::monitor(sims, n.chains=5)