如何为带嵌套随机效应的GLMM固定效应系数生成95%非参数Bootstrap置信区间？

问题背景 #

我正在用带Poisson误差族的广义线性混合模型（GLMM），检验三个二分变量（AGE=ADULT/JUVENILE、SEX=MALE/FEMALE、MEDICATION=NEW/OLD）及AGE与MEDICATION的交互项对计数型响应变量COUNT的显著影响。

数据存在嵌套依赖性：来自22个站点（SITE有33个水平）、21个年份（YEAR为分类变量），且站点采样年份不完全覆盖；同时数据稀疏，每个站点每年的COUNT观测值较少。

构建的Poisson GLMM代码如下：

model <- glmer(data = mydata,
               family = poisson(link = "log"),
               formula = COUNT ~ SEX + SEX:MEDICATION + AGE + AGE:SEX + MEDICATION + AGE:MEDICATION + (1|SITE/YEAR), 
               offset = log(COUNT.SAMPLE.SIZE),
               nAGQ = 0)

尝试用glmmboot包的bootstrap_model生成固定效应系数的95%非参数Bootstrap置信区间时，工具提示“Performing case resampling (no random effects)”，无法识别嵌套随机效应；尝试设置块重采样参数均报错，重构模型后问题依旧。希望获得可行的解决方案，附基于lme4::grouseticks的可复现代码。

解决方案 #

方法1：使用lmeresampler包（推荐） #

lmeresampler支持混合模型的块重采样，能适配嵌套随机效应的依赖性结构。以下是基于grouseticks的完整可复现代码：

# 加载所需包
library(tidyverse)
library(lme4)
library(lmeresampler)

# 加载数据集
data("grouseticks")

# 添加模拟的SEX、AGE、MEDICATION变量
set.seed(1)
grouseticks$SEX <- factor(sample(c("MALE", "FEMALE"), nrow(grouseticks), replace = TRUE))
set.seed(2)
grouseticks$AGE <- factor(sample(c("ADULT", "JUVENILE"), nrow(grouseticks), replace = TRUE))
set.seed(3)
grouseticks$MEDICATION <- factor(sample(c("OLD", "NEW"), nrow(grouseticks), replace = TRUE))

# 计算每个LOCATION-YEAR的样本量（用于偏移项）
grouseticks <- grouseticks %>%
  group_by(LOCATION, YEAR) %>%
  mutate(SAMPLE.SIZE = n()) %>%
  ungroup()

# 构建Poisson GLMM（修正原代码重复SEX项的问题）
model <- glmer(data = grouseticks,
               family = poisson(link = "log"),
               formula = TICKS ~ SEX + AGE + MEDICATION + AGE:MEDICATION + (1|LOCATION/YEAR),
               offset = log(SAMPLE.SIZE),
               nAGQ = 0)

# 执行块重采样Bootstrap：按LOCATION-YEAR组合保留嵌套结构
set.seed(123)
boot_results <- case_bootstrap(model,
                               resample = "block",
                               block_var = c("LOCATION", "YEAR"),
                               B = 1000,
                               type = "perc")

# 查看固定效应的95%置信区间
confint(boot_results)

说明

• block_var = c("LOCATION", "YEAR")指定按站点-年份组合块重采样，保留数据的嵌套依赖性；
• type = "perc"生成百分位数置信区间，也可替换为"bca"（偏差校正加速区间，更稳健但计算稍慢）；
• 输出结果包含所有固定效应系数的95%置信区间。

方法2：使用boot包手动实现（灵活自定义） #

如果需要更精细的控制，可通过boot包手动编写重采样逻辑：

library(boot)

# 定义Bootstrap函数：输入块索引，返回固定效应系数
boot_fun <- function(data, indices) {
  # 提取重采样后的数据集
  resampled_data <- data[indices, ]
  # 尝试拟合模型，捕获拟合失败的情况
  fit <- tryCatch(glmer(data = resampled_data,
                        family = poisson(link = "log"),
                        formula = TICKS ~ SEX + AGE + MEDICATION + AGE:MEDICATION + (1|LOCATION/YEAR),
                        offset = log(SAMPLE.SIZE),
                        nAGQ = 0),
                  error = function(e) NA)
  # 拟合成功则返回系数，否则返回NA
  if (!is.na(fit)) return(fixef(fit)) else return(rep(NA, length(fixef(model))))
}

# 定义站点-年份块
grouseticks$block <- paste(grouseticks$LOCATION, grouseticks$YEAR, sep = "_")
block_ids <- unique(grouseticks$block)
obs_block <- match(grouseticks$block, block_ids)

# 执行Bootstrap
set.seed(123)
boot_obj <- boot(data = grouseticks,
                 statistic = boot_fun,
                 R = 1000,
                 strata = obs_block)

# 计算95%百分位数置信区间
lapply(1:length(fixef(model)), function(i) {
  boot.ci(boot_obj, type = "perc", index = i)
})

说明

• 手动定义站点-年份块，按块分层重采样，确保保留原数据的嵌套结构；
• tryCatch处理数据稀疏导致的模型拟合失败问题，避免中断Bootstrap流程。

关于glmmboot的说明 #

glmmboot的bootstrap_model对嵌套随机效应的支持有限，无法很好适配这类数据的依赖性结构，因此更推荐使用上述两种方法实现块重采样Bootstrap。

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