Algoritma standar untuk melakukan regresi linier hirarkis?

9

Apakah ada algoritma standar (sebagai lawan program) untuk melakukan regresi linier hirarkis? Apakah orang biasanya hanya melakukan MCMC atau ada lebih khusus, mungkin sebagian bentuk tertutup, algoritma?

John Salvatier
sumber

Jawaban:

9

Ada Harvey Goldstein iteratif generalisasi-kuadrat terkecil (IGLS) algoritma untuk satu, dan juga itu modifikasi kecil, dibatasi iteratif umum-kuadrat terkecil (RIGLS), yang memberikan perkiraan yang tidak bias dari parameter varians.

Algoritma ini masih berulang, jadi bukan bentuk tertutup, tetapi secara komputasi lebih sederhana daripada MCMC atau kemungkinan maksimum. Anda hanya mengulanginya sampai parameter bertemu.

Untuk info lebih lanjut tentang ini dan alternatifnya, lihat misalnya:

onestop
sumber
Menakjubkan! Persis apa yang saya cari.
John Salvatier
4

Paket lme4 dalam R menggunakan kuadrat iteratif reweighted terkecil (IRLS) dan menghukum iteratif kuadrat terkecil reweighted (PIRLS). Lihat PDF di sini:

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

Axl
sumber
1
Douglas Bates dan Steven Walker telah menciptakan proyek GitHub yang tujuannya adalah menggunakan kode R murni untuk mengimplementasikan algoritma PIRLS di atas. github.com/lme4/lme4pureR . Jika Anda mempertimbangkan lmer()fungsi dasar dalam lme4paket R, Anda biasanya harus membaca seluruh kode C ++ untuk memahami implementasi PIRLS di lmer()(yang mungkin menantang bagi kita yang tidak begitu berpengalaman dalam pemrograman C ++).
Chris
1

Sumber lain yang baik untuk "algoritma komputasi" untuk HLM (sekali lagi sejauh Anda melihatnya sebagai spesifikasi yang mirip dengan LMM) adalah:

  • McCulloch, C., Searle, S., Neuhaus, J. (2008). Model Linear dan Campuran Umum. Edisi ke-2. Wiley. Bab 14 - Komputasi.

Algoritma yang mereka daftarkan untuk komputasi LMM meliputi:

  • Algoritma EM
  • Algoritma Newton Raphson

Algoritma yang mereka daftarkan untuk GLMM meliputi:

  • Quadrature numerik (GH quadrature)
  • Algoritma EM
  • Algoritma MCMC (seperti yang Anda sebutkan)
  • Algoritma perkiraan stokastik
  • Simulasi kemungkinan maksimum

Algoritma lain untuk GLMM yang mereka sarankan meliputi:

  • Metode kuasi-kemungkinan hukuman
  • Laplace aproksimasi
  • PQL / Laplace dengan koreksi bias bootstrap
Chris
sumber
0

Jika Anda menganggap HLM sebagai jenis model campuran linier, Anda dapat mempertimbangkan algoritma EM. Halaman 22-23 dari catatan kursus berikut menunjukkan bagaimana menerapkan algoritma EM klasik untuk model campuran:

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
sumber