Standard algorithms for doing hierarchical linear regression?

Question

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?

Answer 1

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

For more info on this and alternatives, see e.g.:

• Stephen W. Raudenbush, Anthony S. Bryk. Hierarchical linear models: applications and data analysis methods. (2nd edition) Sage, 2002.

Answer 2

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

http://rss.acs.unt.edu/Rdoc/library/lme4/doc/index.html

Answer 3

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:

• Numerical quadrature (GH quadrature)
• 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

Answer 4

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 )
}