Стандартные алгоритмы для выполнения иерархической линейной регрессии?

9

Существуют ли стандартные алгоритмы (в отличие от программ) для выполнения иерархической линейной регрессии? Люди обычно просто делают MCMC или есть более специализированные, возможно частично закрытые формы, алгоритмы?

Джон Сальватье
источник

Ответы:

9

Существует алгоритм итеративного обобщенного наименьших квадратов (IGLS) Харви Голдштейна для одного, а также его незначительная модификация, ограниченные итеративные обобщенные наименьшие квадраты (RIGLS), который дает несмещенные оценки параметров дисперсии.

Эти алгоритмы все еще итеративны, поэтому не являются закрытой формой, но они вычислительно проще, чем MCMC или с максимальной вероятностью. Вы просто повторяете, пока параметры не сходятся.

  • Гольдштейн Х. Многоуровневый смешанный линейно-модельный анализ с использованием итерационных обобщенных наименьших квадратов. Биометрика 1986; 73 (1): 43-56. doi: 10.1093 / biomet / 73.1.43

  • Гольдштейн Х. Ограниченная несмещенная итерационная обобщенная оценка наименьших квадратов. Биометрика 1989; 76 (3): 622-623. doi: 10.1093 / biomet / 76.3.622

Для получения дополнительной информации об этом и альтернативах, см., Например:

универсальный
источник
Потрясающе! Именно то, что я искал.
Джон Сальватье
4

Пакет lme4 в R использует итеративно взвешенные наименьшие квадраты (IRLS) и штрафует итеративно повторно взвешенные наименьшие квадраты (PIRLS). Смотрите PDF здесь:

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

Axl
источник
1
Дуглас Бейтс и Стивен Уокер создали проект GitHub, целью которого является использование чистого кода R для реализации алгоритма PIRLS, описанного выше. github.com/lme4/lme4pureR . Если вы рассматриваете базовую lmer()функцию в lme4пакете R, вам, как правило, приходится читать весь набор кода C ++, чтобы понять реализацию PIRLS lmer()(что может быть непросто для тех из нас, кто не так хорошо разбирается в программировании на C ++).
Крис
1

Еще один хороший источник «вычислительных алгоритмов» для HLM (опять же, если рассматривать их как спецификации, аналогичные LMM):

  • McCulloch, C., Searle, S., Neuhaus, J. (2008). Обобщенные линейные и смешанные модели. 2-е издание. Wiley. Глава 14 - Компьютеры.

Алгоритмы, которые они перечисляют для вычисления LMM, включают:

  • ЭМ алгоритм
  • Алгоритм Ньютона Рафсона

Алгоритмы, которые они перечисляют для GLMM, включают:

  • Числовая квадратура (GH квадратура)
  • ЭМ алгоритм
  • Алгоритмы MCMC (как вы упоминаете)
  • Алгоритмы стохастической аппроксимации
  • Имитация максимального правдоподобия

Другие алгоритмы для GLMM, которые они предлагают, включают:

  • Оштрафованные квази-правдоподобные методы
  • Аппроксимации Лапласа
  • PQL / Laplace с коррекцией смещения при начальной загрузке
Крис
источник
0

Если вы рассматриваете HLM как тип линейной смешанной модели, вы можете рассмотреть EM-алгоритм. На стр. 22-23 следующих примечаний к курсу показано, как реализовать классический алгоритм EM для смешанной модели:

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 )
}
Крис
источник