Standard algorithms for doing hierarchical linear regression?

Are there standard algorithms (as opposed to programs) for doing hierarchical linear regression? Do people usually just do MCMC or are there more specialized, perhaps partially closed form, algorithms?

There's Harvey Goldstein's iterative generalized least-squares (IGLS) algorithm for one, and also it's minor modification, restricted iterative generalized least-squares (RIGLS), that gives unbiased estimates of the variance parameters.

These algorithms are still iterative, so not closed form, but they're computationally simpler than MCMC or maximum likelihood. You just iterate until the parameters converge.

• Goldstein H. Multilevel Mixed Linear-Model Analysis Using Iterative Generalized Least-Squares. Biometrika 1986; 73(1):43-56. doi: 10.1093/biomet/73.1.43

• Goldstein H. Restricted Unbiased Iterative Generalized Least-Squares Estimation. Biometrika 1989; 76(3):622-623. doi: 10.1093/biomet/76.3.622

• Fabulous! Exactly what I was looking for. Commented Feb 26, 2012 at 18:16

The lme4 package in R uses iteratively reweighted least squares (IRLS) and penalized iteratively reweighted least squares (PIRLS). See the PDF's here:

• Douglas Bates and Steven Walker have created a GitHub project whose goal is to use pure R code to implement the PIRLS algorithm above. github.com/lme4/lme4pureR. If you consider the base lmer() function in the lme4 package of R you would normally have to read through a whole bunch of C++ code understand the implementation of PIRLS in lmer() (which may be challenging for those of us not so well versed in C++ programming). Commented Mar 27, 2014 at 1:07

Another good source for "computing algorithms" for HLM's (again to the extent that you view them as similar specifications as LMM's) would be:

• McCulloch, C., Searle, S., Neuhaus, J. (2008). Generalized Linear and Mixed Models. 2nd Edition. Wiley. Chapter 14 - Computing.

Algorithms they list for computing LMM's include:

• EM algorithm
• Newton Raphson algorithm

Algorithms they list for GLMM's include:

• EM algorithm
• MCMC algorithms (as you mention)
• Stochastic approximation algorithms
• Simulated maximum likelihood

Other algorithms for GLMM's that they suggest include:

• Penalized quasi-likelihood methods
• Laplace approximations
• PQL/Laplace with bootstrap bias correction

If you consider the HLM to be a type of linear mixed model, you could consider the EM algorithm. Page 22-23 of the following course notes indicate how to implement the classic EM algorithm for mixed model:

http://www.stat.ucla.edu/~yuille/courses/stat153/emtutorial.pdf

###########################################################
#     Classical EM algorithm for Linear  Mixed Model      #
###########################################################
em.mixed <- function(y, x, z, beta, var0, var1,maxiter=2000,tolerance = 1e-0010)
{
time <-proc.time()
n <- nrow(y)
q1 <- nrow(z)
conv <- 1
L0 <- loglike(y, x, z, beta, var0, var1)
i<-0
cat("  Iter.       sigma0                 sigma1        Likelihood",fill=T)
repeat {
if(i>maxiter) {conv<-0
break}
V <- c(var1) * z %*% t(z) + c(var0) * diag(n)
Vinv <- solve(V)
xb <- x %*% beta
resid <- (y-xb)
temp1 <- Vinv %*% resid
s0 <- c(var0)^2 * t(temp1)%*%temp1 + c(var0) * n - c(var0)^2 * tr(Vinv)
s1 <- c(var1)^2 * t(temp1)%*%z%*%t(z)%*%temp1+ c(var1)*q1 -
c(var1)^2 *tr(t(z)%*%Vinv%*%z)
w <- xb + c(var0) * temp1
var0 <- s0/n
var1 <- s1/q1
beta <- ginverse( t(x) %*% x) %*% t(x)%*% w
L1 <- loglike(y, x, z, beta, var0, var1)
if(L1 < L0) { print("log-likelihood must increase, llikel <llikeO, break.")
conv <- 0
break
}
i <- i + 1
cat("  ", i,"  ",var0,"  ",var1,"  ",L1,fill=T)
if(abs(L1 - L0) < tolerance) {break}  #check for convergence
L0 <- L1
}
list(beta=beta, var0=var0,var1=var1,Loglikelihood=L0)
}

#########################################################
#  loglike calculates the LogLikelihood for Mixed Model #
#########################################################
loglike<- function(y, x, z, beta, var0, var1)
}
{
n<- nrow(y)
V <- c(var1) * z %*% t(z) + c(var0) * diag(n)
Vinv <- ginverse(V)
xb <- x %*% beta
resid <- (y-xb)
temp1 <- Vinv %*% resid
(-.5)*( log(det(V)) + t(resid) %*% temp1 )
}