非线性混合效应模型（nlme）预测值95%置信区间计算咨询

我尝试采用Ben Bolker分享的Bootstrap方法，为非线性混合效应模型（nlme）生成预测值的95%置信区间，目标是针对预测变量（示例中为year）的每个水平模拟置信区间。

以下是使用FlexParamCurve包的penguin.data数据集的可复现代码，改编自Dr. Bolker的示例：

library(FlexParamCurve) # 同时加载所需的nlme包
library(ggplot2)

set.seed(1234)

## 构建模型
fm2 <- nlme(weight ~ SSlogis(ckage, Asym, R0, lrc),
            data = penguin.data,
            fixed= list(Asym ~ year,
                        R0 ~ year,
                        lrc ~ year),
            random = Asym ~ 1,
            start = c(Asym = 1000, 0,
                      R0 = 21, 0,
                      lrc = 1, 0),
            control = list(maxIter = 100),
            na.action = na.pass)

# 创建用于预测的模拟x（ckage）值
xvals.peng <- with(penguin.data,seq(min(ckage),max(ckage),length.out=100))

nresamp <- 100

## 工具函数：计算置信区间
get_CI <- function(y,pref = "") {
  r1 <- t(apply(y , 1 , quantile , c(0.025 , 0.975)))
  setNames(as.data.frame(r1) , paste0(pref , c("lwr" , "upr"))) 
}

## 创建预测用数据框
pengframe <- with(penguin.data, data.frame(ckage = xvals.peng))

## 尝试预测但失败
pengframe$weight <- predict(fm2,newdata=pengframe,level=0)

运行后出现报错：

Error in eval(predvars, data, env) : object 'year' not found

报错原因是预测数据框pengframe缺少模型中的固定效应变量year。我通过rbind()拼接不同year水平的数据集解决了问题，成功得到各year水平的预测值及对应的95%置信区间：

### 分别创建含不同year水平的数据框后合并（year只有两个水平）
pengframe1 <- with(penguin.data,data.frame(ckage=xvals.peng))
pengframe1$year <- as.factor('2000')
pengframe2 <- with(penguin.data,data.frame(ckage=xvals.peng))
pengframe2$year <- as.factor('2002')

pengframe <- rbind(pengframe1, pengframe2)

# 现在可以按year预测weight
pengframe$weight <- predict(fm2,newdata=pengframe,level=0)

head(pengframe)

sampfun <- function(fitted,data,idvar="bandid") {
  pp <- predict(fitted,levels=1)
  rr <- residuals(fitted)
  dd <- data.frame(data,pred=pp,res=rr)
  ## 有放回抽样分组
  iv <- levels(data[[idvar]])
  bsamp1 <- sample(iv,size=length(iv),replace=TRUE)
  bsamp2 <- lapply(bsamp1,
                   function(x) {
                     ## 组内有放回抽样残差
                     ddb <- dd[dd[[idvar]]==x,]
                     ## Bootstrap响应值 = 预测值 + 抽样残差
                     ddb$height <- ddb$pred +
                       sample(ddb$res,size=nrow(ddb),replace=TRUE)
                     return(ddb)
                   })
  res <- do.call(rbind,bsamp2)  ## 合并结果
  if (is(data,"groupedData"))
    res <- groupedData(res,formula=formula(data))
  return(res)
}

pfun <- function(fm) {
  predict(fm,newdata=pengframe,level=0)
}

yvals2 <- replicate(nresamp, 
                    pfun(update(fm2, 
                                data = sampfun(fm2, 
                                               penguin.data, 
                                               "bandid"))))

peng2 <- get_CI(yvals2,"boot_")
head(peng2)
pengframe <- data.frame(pengframe,peng2)
head(pengframe)

ggplot(pengframe, aes(ckage, weight, color = year)) + 
  geom_smooth() + # 仅作示例，实际数据集使用geom_func()
  geom_ribbon(pengframe, mapping = aes(x = ckage, ymin = boot_lwr, ymax =boot_upr, group=year, fill = year), alpha = 0.3)

现在有两个问题：

1. 这个解决方法是否具备统计合理性？
2. 是否存在更简洁的实现方式（比如在生成xvals.peng时直接添加year变量）？后续我会用sex替代year，且模型包含嵌套随机因子(1 | group/id)。

问题解答 #

1. 解决方法的统计合理性 #

你的方法在统计上是完全合理的：

• 预测时必须提供模型中所有固定效应变量（这里是year），你通过构建包含所有year水平的预测数据框，确保模型能正确计算各水平的群体预测值（level=0），这一步逻辑正确。
• 采用的Bootstrap方法属于残差Bootstrap结合分组抽样，针对混合效应模型，这种方法通过有放回抽样分组（bandid），再在组内抽样残差生成新响应值，能较好保留原数据的分组结构和随机效应分布，符合Ben Bolker提出的思路，适合估计非线性混合模型的预测置信区间。
• 后续用sex替代year、模型加入嵌套随机因子(1 | group/id)时，只需调整抽样逻辑：按最上层的group进行有放回抽样，再在每个group内的id层级处理残差，核心逻辑依然适用。

2. 更简洁的实现方式 #

可以用expand.grid()直接生成包含所有ckage和year组合的预测数据框，无需手动拆分再rbind，代码更简洁：

# 获取year的所有水平
year_levels <- levels(penguin.data$year)
# 生成所有ckage和year的组合
pengframe <- expand.grid(
  ckage = xvals.peng,
  year = year_levels
)
# 直接预测
pengframe$weight <- predict(fm2, newdata = pengframe, level = 0)

如果后续替换为sex，只需把year_levels换成levels(penguin.data$sex)即可。针对嵌套随机因子(1 | group/id)，调整抽样函数示例如下：

sampfun_nested <- function(fitted, data, group_var = "group", id_var = "id") {
  pp <- predict(fitted, levels = 2)  # level对应嵌套结构的最高层，需根据模型调整
  rr <- residuals(fitted)
  dd <- data.frame(data, pred = pp, res = rr)
  
  # 有放回抽样group
  group_levels <- levels(data[[group_var]])
  bsamp_groups <- sample(group_levels, size = length(group_levels), replace = TRUE)
  
  bsamp_list <- lapply(bsamp_groups, function(g) {
    # 抽取该group下的所有id数据
    group_data <- dd[dd[[group_var]] == g, ]
    id_levels <- unique(group_data[[id_var]])
    # 有放回抽样该group内的id
    bsamp_ids <- sample(id_levels, size = length(id_levels), replace = TRUE)
    
    # 处理每个抽样到的id
    id_list <- lapply(bsamp_ids, function(i) {
      id_data <- group_data[group_data[[id_var]] == i, ]
      # 生成Bootstrap响应值
      id_data$height <- id_data$pred + sample(id_data$res, size = nrow(id_data), replace = TRUE)
      return(id_data)
    })
    do.call(rbind, id_list)
  })
  res <- do.call(rbind, bsamp_list)
  if (is(data, "groupedData")) {
    res <- groupedData(res, formula = formula(data))
  }
  return(res)
}

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