Algoritmos estándar para hacer regresión lineal jerárquica?

¿Existen algoritmos estándar (a diferencia de los programas) para hacer una regresión lineal jerárquica? ¿Las personas generalmente solo hacen MCMC o hay algoritmos más especializados, quizás de forma parcialmente cerrada?

Existe el algoritmo iterativo de mínimos cuadrados generalizados (IGLS) de Harvey Goldstein para uno, y también es una modificación menor, mínimos cuadrados generalizados iterativos restringidos (RIGLS), que proporciona estimaciones imparciales de los parámetros de varianza.

Estos algoritmos aún son iterativos, por lo que no son de forma cerrada, pero son computacionalmente más simples que MCMC o de máxima probabilidad. Simplemente itera hasta que los parámetros converjan.

• Goldstein H. Análisis de modelo lineal mixto multinivel utilizando mínimos cuadrados generalizados iterativos. Biometrika 1986; 73 (1): 43-56. doi: 10.1093 / biomet / 73.1.43

• Goldstein H. Estimación de mínimos cuadrados iterativa restringida imparcial imparcial. Biometrika 1989; 76 (3): 622-623. doi: 10.1093 / biomet / 76.3.622

Para obtener más información sobre esto y alternativas, consulte, por ejemplo:

• Stephen W. Raudenbush, Anthony S. Bryk. Modelos lineales jerárquicos: aplicaciones y métodos de análisis de datos . (2a edición) Sage, 2002.

El paquete lme4 en R usa mínimos cuadrados repesados ​​iterativamente (IRLS) y mínimos cuadrados repesados ​​iterativamente penalizados (PIRLS). Vea los PDF aquí:

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

Otra buena fuente de "algoritmos informáticos" para HLM (de nuevo en la medida en que los vea como especificaciones similares a las de LMM) sería:

• McCulloch, C., Searle, S., Neuhaus, J. (2008). Modelos lineales y mixtos generalizados. 2da edición. Wiley Capítulo 14 - Computación.

Los algoritmos que enumeran para computar LMM incluyen:

• Algoritmo EM
• Algoritmo de Newton Raphson

Los algoritmos que enumeran para GLMM incluyen:

• Cuadratura numérica (cuadratura GH)
• Algoritmo EM
• Algoritmos MCMC (como mencionas)
• Algoritmos de aproximación estocástica
• Máxima probabilidad simulada

Otros algoritmos para GLMM que sugieren incluyen:

• Métodos de cuasi-verosimilitud penalizados
• Aproximaciones de Laplace
• PQL / Laplace con corrección de sesgo bootstrap

Si considera que el HLM es un tipo de modelo mixto lineal, podría considerar el algoritmo EM. Las páginas 22-23 de las siguientes notas del curso indican cómo implementar el algoritmo EM clásico para el modelo mixto:

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