Standardowe algorytmy do wykonywania hierarchicznej regresji liniowej?


Odpowiedzi:


9

Istnieje jeden iteracyjny algorytm uogólnionego najmniejszego kwadratu Harveya Goldsteina (IGLS), a także jego niewielka modyfikacja, ograniczona iteracyjna uogólniona uogólniona najmniejszych kwadratów (RIGLS), który daje obiektywne oszacowania parametrów wariancji.

Algorytmy te są wciąż iteracyjne, więc nie są zamknięte, ale są obliczeniowo prostsze niż MCMC lub maksymalne prawdopodobieństwo. Po prostu iterujesz, aż parametry się zbiegną.

  • Goldstein H. Wielopoziomowa mieszana analiza modelu liniowego z wykorzystaniem iteracyjnych uogólnionych najmniejszych kwadratów. Biometrika 1986; 73 (1): 43–56. doi: 10.1093 / biomet / 73.1.43

  • Goldstein H. Ograniczał bezstronną iteracyjną uogólnioną estymację metodą najmniejszych kwadratów. Biometrika 1989; 76 (3): 622–623. doi: 10.1093 / biomet / 76.3.622

Aby uzyskać więcej informacji o tym i alternatywach, zobacz np .:


Fantastyczny! Dokładnie tego szukałem.
John Salvatier

4

Pakiet lme4 w R używa iteracyjnie ponownie ważonych najmniejszych kwadratów (IRLS) i karany iteracyjnie ponownie ważonych najmniejszych kwadratów (PIRLS). Zobacz pliki PDF tutaj:

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


1
Douglas Bates i Steven Walker stworzyli projekt GitHub, którego celem jest użycie czystego kodu R do implementacji powyższego algorytmu PIRLS. github.com/lme4/lme4pureR . Jeśli weźmiesz pod uwagę funkcję podstawową lmer()w lme4pakiecie R, zwykle będziesz musiał przeczytać całą masę kodu C ++, zrozumieć implementację PIRLS w lmer()(co może być trudne dla tych z nas, którzy nie są tak dobrze zaznajomieni z programowaniem w C ++).
Chris

1

Innym dobrym źródłem „algorytmów obliczeniowych” dla HLM (ponownie w zakresie, w którym postrzegasz je jako podobne specyfikacje jak LMM), byłoby:

  • McCulloch, C., Searle, S., Neuhaus, J. (2008). Uogólnione modele liniowe i mieszane. 2. edycja. Wiley. Rozdział 14 - Informatyka.

Wymienione przez nich algorytmy obliczania LMM obejmują:

  • Algorytm EM
  • Algorytm Newtona Raphsona

Algorytmy, które wymieniają dla GLMM, obejmują:

  • Kwadratura numeryczna (kwadratura GH)
  • Algorytm EM
  • Algorytmy MCMC (jak wspomniałeś)
  • Stochastyczne algorytmy aproksymacyjne
  • Symulowane maksymalne prawdopodobieństwo

Inne sugerowane przez nich algorytmy GLMM obejmują:

  • Karane metody quasi-prawdopodobieństwa
  • Przybliżenia Laplace'a
  • PQL / Laplace z korekcją błędu początkowego

0

Jeśli uważasz HLM za rodzaj liniowego modelu mieszanego, możesz rozważyć algorytm EM. Strony 22-23 poniższych notatek kursowych wskazują, jak zaimplementować klasyczny algorytm EM dla modelu mieszanego:

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 )
}
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.