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?

score 9 · Accepted Answer · answered Feb 26 '12 at 14:22

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.

Fabulous! Exactly what I was looking for. – John Salvatier Feb 26 '12 at 18:16 — John Salvatier, Feb 26 '12 at 18:16

Axl · Answer 2 · 2012-03-22T23:29:42.933

4

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

edited Mar 22 '12 at 23:29

answered Mar 03 '12 at 00:56

Axl

151
3

1

Douglas Bates and Steven Walker have created a GitHub project whose goal is to use pure R code to implement the PIRLS algorithm above. https://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). – Chris Mar 27 '14 at 01:07

score 1 · Answer 3 · answered Mar 27 '14 at 01:17

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

score 0 · Answer 4 · answered Mar 27 '14 at 01:05

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

Standard algorithms for doing hierarchical linear regression?

4 Answers4