有一个 Harvey Goldstein 的迭代广义最小二乘 (IGLS) 算法，它也是一个小的修改，限制迭代广义最小二乘 (RIGLS)，它给出了方差参数的无偏估计。

这些算法仍然是迭代的，所以不是封闭形式，但它们在计算上比 MCMC 或最大似然更简单。您只需迭代直到参数收敛。

• Goldstein H. 使用迭代广义最小二乘的多级混合线性模型分析。生物计量学1986；73（1）：43-56。doi: 10.1093/biomet/73.1.43

• Goldstein H. 受限无偏迭代广义最小二乘估计。生物计量学1989；76（3）：622-623。doi: 10.1093/biomet/76.3.622

有关此选项和替代方案的更多信息，请参见例如：

• 斯蒂芬·W·劳登布什、安东尼·S·布瑞克。分层线性模型：应用和数据分析方法。（第 2 版）圣人，2002 年。

R 中的 lme4 包使用迭代重加权最小二乘法 (IRLS) 和惩罚迭代重加权最小二乘法 (PIRLS)。在此处查看 PDF：

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

HLM 的“计算算法”的另一个很好的来源（在某种程度上，您将它们视为与 LMM 相似的规范）将是：

• McCulloch, C., Searle, S., Neuhaus, J. (2008)。广义线性和混合模型。第 2 版。威利。第 14 章 - 计算。

他们列出的用于计算 LMM 的算法包括：

• 电磁算法
• 牛顿拉夫森算法

他们为 GLMM 列出的算法包括：

• 数值求积（GH求积）
• 电磁算法
• MCMC算法（正如你提到的）
• 随机逼近算法
• 模拟最大似然

他们建议的其他 GLMM 算法包括：

• 惩罚准似然法
• 拉普拉斯近似
• 带自举偏差校正的 PQL/Laplace

如果您认为 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 )
}