Algoritmos padrão para fazer regressão linear hierárquica?

9

Existem algoritmos padrão (em oposição a programas) para fazer regressão linear hierárquica? As pessoas geralmente fazem apenas o MCMC ou existem algoritmos mais especializados, talvez de forma parcialmente fechada?

John Salvatier
fonte

Respostas:

9

Há o algoritmo de mínimos quadrados generalizados iterativos (IGLS) de Harvey Goldstein para um, e também é uma modificação menor, os mínimos quadrados generalizados iterativos restritos (RIGLS), que fornece estimativas imparciais dos parâmetros de variação.

Esses algoritmos ainda são iterativos, portanto não são de forma fechada, mas são computacionalmente mais simples que o MCMC ou a probabilidade máxima. Você apenas itera até os parâmetros convergirem.

  • Goldstein H. Análise de modelo linear misto multinível usando mínimos quadrados generalizados iterativos. Biometrika 1986; 73 (1): 43-56. doi: 10.1093 / biomet / 73.1.43

  • Goldstein H. Estimativa de Mínimos Quadrados Generalizados Iterativos Restritos e Imparciais. Biometrika 1989; 76 (3): 622-623. doi: 10.1093 / biomet / 76.3.622

Para mais informações sobre isso e alternativas, consulte, por exemplo:

uma parada
fonte
Fabuloso! Exatamente o que eu estava procurando.
John Salvatier
4

O pacote lme4 em R usa mínimos quadrados ponderados iterativamente (IRLS) e mínimos quadrados ponderados iterativamente penalizados (PIRLS). Veja o PDF aqui:

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

Axl
fonte
11
Douglas Bates e Steven Walker criaram um projeto GitHub cujo objetivo é usar o código R puro para implementar o algoritmo PIRLS acima. github.com/lme4/lme4pureR . Se você considerar a lmer()função base no lme4pacote do R, normalmente precisará ler um monte de códigos C ++ para entender a implementação do PIRLS lmer()(o que pode ser um desafio para nós que não são tão versados ​​na programação C ++).
27414 Chris
1

Outra boa fonte para "algoritmos de computação" para HLMs (novamente na medida em que você os vê como especificações semelhantes às de LMMs) seria:

  • McCulloch, C., Searle, S., Neuhaus, J. (2008). Modelos Lineares e Mistos Generalizados. 2ª Edição. Wiley. Capítulo 14 - Computação.

Os algoritmos que eles listam para calcular os LMMs incluem:

  • Algoritmo EM
  • Algoritmo de Newton Raphson

Os algoritmos que eles listam para os GLMMs incluem:

  • Quadratura numérica (quadratura GH)
  • Algoritmo EM
  • Algoritmos MCMC (como você mencionou)
  • Algoritmos de aproximação estocástica
  • Probabilidade máxima simulada

Outros algoritmos para GLMMs que eles sugerem incluem:

  • Métodos penalizados de quase-probabilidade
  • Aproximações de Laplace
  • PQL / Laplace com correção de viés de bootstrap
Chris
fonte
0

Se você considerar o HLM um tipo de modelo misto linear, poderá considerar o algoritmo EM. As páginas 22-23 das seguintes notas de curso indicam como implementar o algoritmo EM clássico para modelo misto:

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