在这里，我展示了三种获取这些值的方法：

• 使用lme4:::coef.merMod()函数的修改/破解版本，将固定效应分量和随机效应分量的方差相加，并计算正态CI。这假定（1）这些分量的抽样分布呈正态分布，以及（2）固定效应和随机效应的方差是独立的。（也不做任何有限大小的修正）
• 使用bootMer()的参数引导。这是假设模型是正确的
• 使用lmeresampler的非参数 Bootstrap 。这使得假设较弱，但很难实现具有（例如）交叉随机效应的情况（该软件包提供了一些替代方法，例如残差或Wild bootstrapping...）

有一个合理的讨论非参数与参数引导here
所有三种方法在这个（表现良好的）例子中给予了 * 合理 * 相似的结果......
下面的代码有点难看，可能可以清理一下（我主要避免tidyverse，这实际上可能是有用的...）

加载包并适配基本模型

library(lme4)
library(ggplot2)
library(lmeresampler)
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)

参数化引导

f <- function(x) unlist(coef(x)$Subject)
b <- bootMer(fm1, FUN = f, seed = 101, nsim = 1000,
             use.u = TRUE,  ## important - *don't* resample RE values
             parallel = "multicore", ncpus = 5)
## rearrange
bootmer_res <- t(apply(as.data.frame(b), 2, quantile, c(0.025, 0.975)))
bootmer_res <- setNames(as.data.frame(bootmer_res), c("lwr", "upr"))

方差和

cc <- coef2.merMod(fm1)  ## defined below
coef_res <- data.frame(unlist(cc$values$Subject + qnorm(0.025)*cc$sd$Subject),
                       unlist(cc$values$Subject + qnorm(0.975)*cc$sd$Subject))
coef_res <- setNames(as.data.frame(coef_res), c("lwr", "upr"))

非参数引导

set.seed(101)
dd <- bootstrap(fm1, .f = f, type = "case", B = 1000, 
     ## resample only within Subjects
     resample = c(FALSE, TRUE))
lmeresamp_res <- t(apply(dd$replicates, 2, quantile, c(0.025, 0.975)))
lmeresamp_res <- setNames(as.data.frame(lmeresamp_res), c("lwr", "upr"))

合并估计和绘图

all_res <- rbind(transform(bootmer_res, nm = "bootmer"),
                 transform(coef_res, nm = "coef2"),
                 transform(lmeresamp_res, nm = "lmeresamp"))
## better to do this by manipulating existing metadata (e.g. rownames)
nsubj <- 18
all_res$var <- rep(rep(c("(Intercept)", "Days"), each = nsubj), 3)
all_res$subj <- rep(seq(nsubj), 6)
ggplot(all_res, aes(xmin = lwr, xmax = upr, y = subj, colour = nm)) +
    facet_wrap(~var, scale = "free") +
    geom_linerange(position = position_dodge(width = 0.5))
                  
ggsave("coef_SDs.png")

coef2.merMod（）函数

coef2.merMod <- function (object, ...)  {
    fef <- data.frame(rbind(fixef(object)), check.names = FALSE)
    fef_var <- data.frame(rbind(diag(vcov(object))), check.names = FALSE) ## *
    ref <- ranef(object, condVar = TRUE)  ## *
    refnames <- unlist(lapply(ref, colnames))
    nmiss <- length(missnames <- setdiff(refnames, names(fef)))
    if (nmiss > 0) {
        fillvars <- setNames(data.frame(rbind(rep(0, nmiss))), 
            missnames)
        fef <- cbind(fillvars, fef)
        fef_var <- cbind(fillvars, ref_var)
    }
    dfun <- function(x) lapply(ref, function(y) x[rep.int(1L, nrow(y)), , drop = FALSE])
    val <- dfun(fef)
    val_var <- dfun(fef_var)
    for (i in seq_along(val)) {
        refi <- ref[[i]]
        refi_var <- t(apply(attr(refi, "postVar"), 3, diag))
        colnames(refi_var) <- colnames(refi)
        row.names(val[[i]]) <- row.names(val_var[[i]]) <- row.names(refi) ## *
        nmsi <- colnames(refi)
        if (!all(nmsi %in% names(fef))) 
            stop("unable to align random and fixed effects")
        ## 
        for (nm in nmsi) {
            val[[i]][[nm]] <- val[[i]][[nm]] + refi[, nm]
            val_var[[i]][[nm]] <- val_var[[i]][[nm]] + refi_var[, nm]
        }
        val_var[[i]] <- as.data.frame(lapply(val_var[[i]], sqrt),
                                      check.names = FALSE)
    }
    list(values = val, sd = val_var)
}