Statistical Modeling, Causal Inference, and Social Science

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)