非线性混合效应模型的Bootstrap置信区间获取与绘图问题

问题概述 #

使用lme4包的nlmer()构建带固定渐近线的Gompertz非线性混合效应模型，模型本身收敛正常，但无法获取Bootstrap置信区间用于geom_ribbon()绘制曲线置信区间：

• 调用confint()时报错：Error: step factor reduced below 0.001 without reducing pwrss
• 调用bootMer()时报错：Error in model.matrix.default(eval(substitute(~foo, list(foo = x[[2]]))), : model frame and formula mismatch in model.matrix()

示例数据 #

# 创建结构相似的示例数据
x_data <- seq(1, 14, length.out = 200)
y_data <- 53 * exp(-exp(-0.14 *(x_data - 13))) + rnorm(length(x_data), mean = 0, sd = 0.5)
individ_ID <- c(1:200) # 个体ID
group_ID <- c(1:10) # 组ID 
group_ID2 <- rep(group_ID, each = 4)
mydata1 <- data.frame(x_data, y_data, individ_ID, group_ID2)

x_data <- seq(15, 28, length.out = 200)
y_data <- 53 * exp(-exp(-0.14 *(x_data - 13))) + rnorm(length(x_data), mean = 0, sd = 0.5)
individ_ID <- c(1:200)
group_ID <- c(1:10)
group_ID2 <- rep(group_ID, each = 4)
mydata2 <- data.frame(x_data, y_data, individ_ID, group_ID2)

mydata3 <- rbind(mydata1, mydata2)

模型构建代码 #

# 定义Gompertz模型函数（固定渐近线A=53）
gompertz_pp <- ~53*exp(-exp(-k*(d - Ti)))

# 生成带导数的模型函数，用于nlmer拟合
fn.gompertz.pp <- deriv(gompertz_pp,
                        namevec=c("k","Ti"), 
                        function.arg=c("d","k","Ti"))

# 设置初始值
startvec <- c(k = 0.15, Ti = 12)

# 拟合非线性混合效应模型
mod <- nlmer(y_data ~ fn.gompertz.pp(x_data, k, Ti) ~ (Ti|group_ID2/individ_ID) + (k|group_ID2/individ_ID), 
             data = mydata3, start = startvec)

解决方案 #

1. 修正bootMer()调用方式 #

非线性混合模型使用bootMer()时，直接传入fixef作为FUN可能触发公式匹配错误，需自定义FUN函数，并选择参数化Bootstrap（比非参数化更稳定）：

# 自定义Bootstrap函数，返回固定效应参数
boot_fun <- function(model) {
  return(fixef(model))
}

# 运行参数化Bootstrap，设置并行加速（可选）
library(parallel)
boot_result <- bootMer(mod, FUN = boot_fun, nsim = 100, type = "parametric",
                       parallel = "multicore", ncpus = detectCores()-1)

# 提取Bootstrap置信区间
boot_ci <- confint(boot_result, level = 0.95)

2. 优化confint()拟合稳定性 #

confint()报错是因为Bootstrap迭代中部分模型拟合失败，可通过更换优化器、增加迭代次数解决：

# 重新拟合模型，指定更稳定的优化器
mod_optim <- nlmer(y_data ~ fn.gompertz.pp(x_data, k, Ti) ~ (Ti|group_ID2/individ_ID) + (k|group_ID2/individ_ID), 
                   data = mydata3, start = startvec,
                   control = nlmerControl(optimizer = "optimx", 
                                          optCtrl = list(method = "L-BFGS-B", maxit = 1000)))

# 调用confint进行Bootstrap
confint_ci <- confint(mod_optim, method = "boot", nsim = 100, 
                      parallel = "multicore", ncpus = detectCores()-1)

3. 基于参数抽样的替代置信区间方法 #

如果Bootstrap始终失败，可通过从固定效应协方差矩阵中抽样参数，结合模型生成预测值的分布来计算置信区间：

library(MASS)
library(dplyr)
library(ggplot2)

# 提取固定效应估计值和协方差矩阵
fixef_vals <- fixef(mod)
vcov_mat <- vcov(mod)

# 从多元正态分布中抽样参数（1000次）
set.seed(123)
param_samples <- MASS::mvrnorm(n = 1000, mu = fixef_vals, Sigma = vcov_mat)

# 生成用于预测的x序列
new_x <- seq(min(mydata3$x_data), max(mydata3$x_data), length.out = 100)

# 生成所有抽样参数对应的预测值
pred_df <- expand.grid(x = new_x, sample_id = 1:1000)
pred_df$y_pred <- apply(pred_df, 1, function(row) {
  k_val <- param_samples[row["sample_id"], "k"]
  ti_val <- param_samples[row["sample_id"], "Ti"]
  # 代入Gompertz模型计算预测值
  53 * exp(-exp(-k_val*(row["x"] - ti_val)))
})

# 计算每个x点的95%置信区间
ci_df <- pred_df %>%
  group_by(x) %>%
  summarise(
    fit = mean(y_pred),
    lower_ci = quantile(y_pred, 0.025),
    upper_ci = quantile(y_pred, 0.975)
  )

# 绘制原始数据、拟合曲线和置信区间
ggplot() +
  geom_point(data = mydata3, aes(x = x_data, y = y_data), alpha = 0.3, size = 0.8) +
  geom_line(data = ci_df, aes(x = x, y = fit), color = "#E63946", linewidth = 1) +
  geom_ribbon(data = ci_df, aes(x = x, ymin = lower_ci, ymax = upper_ci), 
              fill = "#E63946", alpha = 0.2) +
  labs(x = "x", y = "y") +
  theme_minimal()

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