Czy istnieją standardowe algorytmy (w przeciwieństwie do programów) do wykonywania hierarchicznej regresji liniowej? Czy ludzie zwykle po prostu wykonują MCMC, czy są bardziej wyspecjalizowane, być może częściowo zamknięte, algorytmy?
Czy istnieją standardowe algorytmy (w przeciwieństwie do programów) do wykonywania hierarchicznej regresji liniowej? Czy ludzie zwykle po prostu wykonują MCMC, czy są bardziej wyspecjalizowane, być może częściowo zamknięte, algorytmy?
Odpowiedzi:
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 .:
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:
lmer()
w lme4
pakiecie 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 ++).
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:
Wymienione przez nich algorytmy obliczania LMM obejmują:
Algorytmy, które wymieniają dla GLMM, obejmują:
Inne sugerowane przez nich algorytmy GLMM obejmują:
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 )
}