Standardalgorithmen für die hierarchische lineare Regression?

9

Gibt es Standardalgorithmen (im Gegensatz zu Programmen) für die hierarchische lineare Regression? Machen die Leute normalerweise nur MCMC oder gibt es spezialisiertere, vielleicht teilweise geschlossene Algorithmen?

John Salvatier
quelle

Antworten:

9

Zum einen gibt es Harvey Goldsteins iterativen Algorithmus für generalisierte kleinste Quadrate (IGLS) und eine geringfügige Modifikation, eingeschränkte iterative generalisierte kleinste Quadrate (RIGLS), die unvoreingenommene Schätzungen der Varianzparameter liefert.

Diese Algorithmen sind immer noch iterativ, also keine geschlossene Form, aber rechnerisch einfacher als MCMC oder maximale Wahrscheinlichkeit. Sie iterieren einfach, bis die Parameter konvergieren.

  • Goldstein H. Mehrstufige gemischte lineare Modellanalyse unter Verwendung iterativer verallgemeinerter kleinster Quadrate. Biometrika 1986; 73 (1): 43-56. doi: 10.1093 / biomet / 73.1.43

  • Goldstein H. Restricted Unbias Iterative Generalized Least-Squares-Schätzung. Biometrika 1989; 76 (3): 622 & ndash; 623. doi: 10.1093 / biomet / 76.3.622

Weitere Informationen hierzu und zu Alternativen finden Sie zB:

ein Stop
quelle
Fabelhaft! Genau das, wonach ich gesucht habe.
John Salvatier
4

Das lme4-Paket in R verwendet iterativ neu gewichtete kleinste Quadrate (IRLS) und bestrafte iterativ neu gewichtete kleinste Quadrate (PIRLS). Sehen Sie die PDFs hier:

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

Axl
quelle
1
Douglas Bates und Steven Walker haben ein GitHub-Projekt erstellt, dessen Ziel es ist, reinen R-Code zur Implementierung des obigen PIRLS-Algorithmus zu verwenden. github.com/lme4/lme4pureR . Wenn Sie die lmer()Basisfunktion im lme4Paket von R berücksichtigen, müssen Sie normalerweise eine ganze Reihe von C ++ - Code durchlesen, um die Implementierung von PIRLS in zu verstehen lmer()(was für diejenigen von uns, die sich mit C ++ - Programmierung nicht so gut auskennen, eine Herausforderung sein kann).
Chris
1

Eine weitere gute Quelle für "Berechnungsalgorithmen" für HLMs (wiederum in dem Maße, in dem Sie sie als ähnliche Spezifikationen wie LMMs betrachten) wäre:

  • McCulloch, C., Searle, S., Neuhaus, J. (2008). Verallgemeinerte lineare und gemischte Modelle. 2. Auflage. Wiley. Kapitel 14 - Rechnen.

Zu den Algorithmen, die sie zur Berechnung von LMMs auflisten, gehören:

  • EM-Algorithmus
  • Newton Raphson-Algorithmus

Zu den Algorithmen, die sie für GLMMs auflisten, gehören:

  • Numerische Quadratur (GH-Quadratur)
  • EM-Algorithmus
  • MCMC-Algorithmen (wie Sie erwähnen)
  • Stochastische Approximationsalgorithmen
  • Simulierte maximale Wahrscheinlichkeit

Andere Algorithmen für GLMMs, die sie vorschlagen, umfassen:

  • Bestrafte Quasi-Likelihood-Methoden
  • Laplace-Annäherungen
  • PQL / Laplace mit Bootstrap-Bias-Korrektur
Chris
quelle
0

Wenn Sie das HLM als eine Art lineares gemischtes Modell betrachten, können Sie den EM-Algorithmus in Betracht ziehen. Seite 22-23 der folgenden Kursnotizen zeigen, wie der klassische EM-Algorithmus für gemischte Modelle implementiert wird:

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
quelle