如何在R中计算混合效应模型受试者总效应的置信区间？

你通过coef(model)获取的是固定效应+随机效应偏移的受试者水平系数（即BLUPs与固定效应的叠加值），以下是几种在R中计算该系数置信区间的实用方法：

方法1：借助merTools包的predictInterval #

merTools包专门针对混合模型提供了便捷的区间计算工具，predictInterval可直接生成受试者水平系数的置信区间：

library(merTools)
# 构造包含所有受试者的新数据集（按需设置x1、x2的取值，这里示例用均值）
newdat <- expand.grid(subject = unique(model@frame$subject),
                      x1 = mean(model@frame$x1),
                      x2 = mean(model@frame$x2))
# 计算95%置信区间，type="coef"指定针对受试者水平系数
ci_results <- predictInterval(model, newdata = newdat, type = "coef", level = 0.95)

输出结果会包含每个受试者各系数的置信上下限。

方法2：手动基于方差-协方差矩阵计算 #

利用lme4自带函数提取方差信息，结合正态近似手动计算区间：

# 提取固定效应值
fixef_vals <- fixef(model)
# 提取随机效应的BLUPs
ranef_vals <- ranef(model)$subject
# 合并得到每个受试者的总系数
coef_total <- t(fixef_vals + t(ranef_vals))

# 提取随机效应的方差-协方差矩阵
vc_ran <- VarCorr(model)$subject
# 提取固定效应的方差-协方差矩阵
vc_fix <- vcov(model)

# 计算标准误（固定效应方差+随机效应方差的平方根）
se_vals <- sqrt(diag(vc_fix) + diag(vc_ran))
# 计算95%置信区间（正态近似用1.96分位数）
ci_lower <- coef_total - 1.96 * se_vals
ci_upper <- coef_total + 1.96 * se_vals

注意：该方法基于正态分布假设，样本量较小时结果可能不够稳健。

方法3：用bootMer做bootstrap抽样 #

通过bootstrap重复抽样获取系数的经验分布，进而得到更可靠的置信区间：

# 定义bootstrap抽样时的提取函数：返回每个受试者的总系数
boot_coef_fun <- function(model) {
  cbind(fixef(model) + t(ranef(model)$subject))
}

# 执行bootstrap抽样，这里设置1000次（可根据计算资源调整）
boot_results <- bootMer(model, FUN = boot_coef_fun, nsim = 1000)

# 用分位数法计算95%置信区间
ci_bootstrap <- apply(boot_results$t, 2, quantile, probs = c(0.025, 0.975))

该方法无需分布假设，适合小样本或非正态数据，但计算耗时较长。

方法4：使用emmeans分析特定协变量组合的效应 #

如果需要针对x1、x2的特定取值组合计算受试者效应的置信区间，emmeans包可以快速实现：

library(emmeans)
# 指定x1、x2的取值，按受试者分组计算边际效应的置信区间
emm_obj <- emmeans(model, ~ x1*x2 | subject, at = list(x1 = c(0,1), x2 = c(0,1)))
# 输出置信区间
confint(emm_obj)

该方法适合分析不同协变量水平下的受试者特异性效应。

内容的提问来源于stack exchange，提问作者Raed Hamed