我使用 R 中的 lme4 包来拟合混合效果模型:
library(lme4)
model <- lmer(y~x1*x2+(1|subect)+(0 + x1+x2+x1:x2|subject)
如果我想获得每个受试者的总效应,同时考虑固定和随机斜率,我可以使用函数 coef 来报告每个受试者的系数。
coef(model)
现在,我想计算每个受试者的估计值的置信区间,同时考虑固定效应和随机效应 CI。你有关于如何在 R 中做到这一点的提示吗?
谢谢!
更新 2024/25:另请参阅此答案。
这里我展示了三种获取这些值的方法:
lme4:::coef.merMod()
函数的修改/破解版本添加固定和随机效应分量的方差并计算正态 CI。这假设 (1) 这些分量的抽样分布是正态的,并且 (2) 固定效应和随机效应的方差是独立的。 (也不进行任何有限尺寸的修正)bootMer()
进行参数引导。这假设模型是正确的lmeresampler
进行非参数引导。这使得假设较弱,但很难实现具有(例如)交叉随机效应的情况(该包提供了一些替代方法,例如残差或狂野引导......)对非参数与参数引导程序进行了合理的讨论 这里
在这个(表现良好)示例中,所有三种方法都给出了合理地相似的结果...
下面的代码有点难看,可能可以清理(我主要避免使用 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 <- 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)
}