如何在R中为gnls拟合的广义非线性模型计算并绘制95%置信区间？

背景与问题 #

我用nlme包的gnls()构建了指数衰减形式的广义非线性最小二乘模型，选择它而非基础nls()是因为支持异方差建模，不用对数据做转换处理。模型示例如下：

model <- gnls(Response ~ C * exp(k * Explanatory1) + A, 
              start = list(C = c(C1,C1), k = c(k1,k1), A = c(A1,A1)), 
              params = list(C ~ Explanatory2, k ~ Explanatory2, A ~ Explanatory2), 
              weights = varPower(), data = Data)

之前在gls()上能用的置信区间计算方法（基于模型矩阵和方差协方差矩阵），对gnls()完全无效。查资料后发现investr::predFit()或drc::drm()+predict()也不支持gnls对象。我试过用confint()获取单个参数的置信区间，再代入参数极值计算预测上下限，但不确定这是不是最优方案。

可行解决方案：非参数Bootstrap #

参考Ben Bolker的思路，我用非参数Bootstrap解决了这个问题——既能处理带分类解释变量的复杂gnls模型，又能得到可靠的预测置信区间，结合ggplot2的可视化效果也很好。

步骤1：模拟带异方差的数据集（含分类变量） #

# 生成数据
set.seed(2)
x <- rep(sample(1:100, 9), 12)
set.seed(15)
r <- rnorm(x, mean = 0, sd = 200*x^-0.8)
y <- c(200, 300) * exp(c(-0.08, -0.05)*x) + c(120, 100) + r
df <- data.frame(x = x, y = y, group = rep(letters[1:2], length.out = length(x)))

步骤2：拟合gnls模型 #

先通过nls()的SSasymp()获取初始值，再构建带异方差权重和分类参数的gnls模型：

# 寻找初始值
m <- nls(y ~ SSasymp(x, A, C, logk))
summary(m) # A = 108.9860, C = 356.6851, k = -2.9356

# 拟合广义非线性最小二乘模型
require(nlme)
mgnls <- gnls(y ~ C * exp(k * x) + A, 
              start = list(C = c(356.6851,356.6851), 
                          k = c(-exp(-2.9356),-exp(-2.9356)), 
                          A = c(108.9860,108.9860)), 
              params = list(C ~ group, k ~ group, A ~ group), 
              weights = varExp(), data = df)
plot(mgnls) # 方差更均匀，说明异方差建模有效

步骤3：用Bootstrap计算置信区间 #

编写Bootstrap函数，重复采样并拟合模型，生成预测值的分布，最后取分位数得到置信上下限：

# 准备新数据用于预测
new <- data.frame(x = c(1:100, 1:100), group = rep(letters[1:2], each = 100))
new$fit <- predict(mgnls, newdata = new)

# Bootstrap函数
bootfun <- function(newdata) {
  start <- coef(mgnls)
  # 有放回采样数据集
  dfboot <- df[sample(nrow(df), size = nrow(df), replace = TRUE),]
  # 更新模型拟合，捕获拟合失败的情况
  bootfit <- try(update(mgnls, start = start, data = dfboot), silent = TRUE)
  if(inherits(bootfit, "try-error")) return(rep(NA, nrow(newdata)))
  predict(bootfit, newdata)
}

# 运行500次Bootstrap
set.seed(10)
bmat <- replicate(500, bootfun(new))
# 计算95%置信区间的上下限
new$lwr <- apply(bmat, 1, quantile, 0.025, na.rm = TRUE)
new$upr <- apply(bmat, 1, quantile, 0.975, na.rm = TRUE)

步骤4：ggplot2可视化 #

把原始数据、预测曲线和置信区间一起绘制：

require(ggplot2)
ggplot() + 
  geom_point(data = df, aes(x, y, colour = group)) + 
  geom_ribbon(data = new, aes(x = x, ymin = lwr, ymax = upr, fill = group), alpha = 0.3) + 
  geom_line(data = new, aes(x, fit, colour = group)) + 
  theme_minimal()

为什么这个方法更好？ #

• 相比用单个参数置信区间极值计算的方法，Bootstrap考虑了参数之间的相关性，得到的置信区间更准确。
• 能灵活处理带分类变量、异方差的复杂gnls模型，适用性广。
• 可视化效果清晰，能直观展示不同组的拟合趋势和不确定性。

内容的提问来源于stack exchange，提问作者Luka Seamus Wright