计算并绘制广义非线性模型的 95% 置信区间

Question

我已经使用R 包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)

简单的nls() 模型的主要区别在于weights 参数，它可以通过解释变量对异方差建模。与gnls() 等价的线性是广义最小二乘法，它使用nlme 的gls() 函数运行。

现在我想计算R 中的置信区间，并将它们与我的模型一起绘制在ggplot()（ggplot2 包）中。我为gls() 对象执行此操作的方式是：

NewData <- data.frame(Explanatory1 = c(...), Explanatory2 = c(...)) 
NewData$fit <- predict(model, newdata = NewData)

到目前为止，一切正常，我的模型适合我。

modmat <-  model.matrix(formula(model)[-2], NewData)
int <- diag(modmat %*% vcov(model) %*% t(modmat))
NewData$lo <- with(NewData, fit - 1.96*sqrt(int))
NewData$hi <- with(NewData, fit + 1.96*sqrt(int))

这部分不适用于gnls()，所以我无法获得我的上下模型预测。

由于这似乎不适用于gnls() 对象，因此我查阅了教科书以及以前提出的问题，但似乎没有一个适合我的需要。我发现的唯一类似问题是How to calculate confidence intervals for Nonlinear Least Squares in r?。在最佳答案中，建议使用investr::predFit() 或使用drc::drm() 构建模型，然后使用常规predict() 函数。这些解决方案都不能帮助我处理gnls()。

我目前的最佳解决方案是使用confint() 函数计算所有三个参数（C、k、A）的 95% 置信区间，然后为置信上限和下限编写两个单独的函数，即一个使用 Cmin， kmin 和 Amin 和一个使用 Cmax、kmax 和 Amax。然后我使用这些函数来预测值，然后用ggplot() 绘制这些值。但是，我对结果并不完全满意，并且不确定这种方法是否最佳。

这是一个最小的可重现示例，为简单起见，忽略了第二个分类解释变量：

# generate data
set.seed(10)
x <-  rep(1:100,2)
r <- rnorm(x, mean = 10, sd = sqrt(x^-1.3))
y <- exp(-0.05*x) + r
df <-  data.frame(x = x, y = y)

# find starting values
m <- nls(y ~ SSasymp(x, A, C, logk))
summary(m) # A = 9.98071, C = 10.85413, logk = -3.14108
plot(m) # clear heteroskedasticity

# fit generalised nonlinear least squares
require(nlme)
mgnls <- gnls(y ~ C * exp(k * x) + A, 
              start = list(C = 10.85413, k = -exp(-3.14108), A = 9.98071),
              weights = varExp(),
              data = df)
plot(mgnls) # more homogenous

# plot predicted values 
df$fit <- predict(mgnls)
require(ggplot2)
ggplot(df) +
  geom_point(aes(x, y)) +
  geom_line(aes(x, fit)) +
  theme_minimal()

按照 Ben Bolker 的回答进行编辑

应用于第二个模拟数据集的标准非参数自举解决方案，它更接近我的原始数据并包含第二个分类解释变量：

# generate data
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)))

# find starting values
m <- nls(y ~ SSasymp(x, A, C, logk))
summary(m) # A = 108.9860, C = 356.6851, k = -2.9356
plot(m) # clear heteroskedasticity

# fit generalised nonlinear least squares
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) # more homogenous

# calculate predicted values 
new <- data.frame(x = c(1:100, 1:100),
                  group = rep(letters[1:2], each = 100))
new$fit <- predict(mgnls, newdata = new)

# calculate bootstrap confidence intervals
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)
}

set.seed(10)
bmat <- replicate(500, bootfun(new))
new$lwr <- apply(bmat, 1, quantile, 0.025, na.rm = TRUE)
new$upr <- apply(bmat, 1, quantile, 0.975, na.rm = TRUE)

# plot data and predictions
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()

这是生成的情节，看起来很整洁！

Answer 1

我实施了一个引导解决方案。我最初进行了标准的非参数引导，对观察结果进行了重新采样，但这会产生 95% 的 CI，看起来很宽——我认为这是因为这种引导形式无法维持x 分布中的平衡（例如，通过重新采样，您最终可能无法观察到 x 的小值）。 （也有可能是我的代码中只有一个错误。）

作为第二个镜头，我切换到从初始拟合重新采样 残差，并将它们添加到预测值；这是一种相当标准的方法，例如在引导时间序列中（尽管我忽略了残差中自相关的可能性，这需要块引导）。

这是基本的引导重采样器。

df$res <- df$y-df$fit
bootfun <- function(newdata=df, perturb=0, boot_res=FALSE) {
    start <- coef(mgnls)
    ## if we start exactly from the previously fitted coefficients we end
    ## up getting all-identical answers? Not sure what's going on here, but
    ## we can fix it by perturbing the starting conditions slightly
    if (perturb>0) {
        start <- start * runif(length(start), 1-perturb, 1+perturb)
    }
    if (!boot_res) {
        ## bootstrap raw data
        dfboot <- df[sample(nrow(df),size=nrow(df), replace=TRUE),]
    } else {
        ## bootstrap residuals
        dfboot <- transform(df,
                            y=fit+sample(res, 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=newdata)
}

set.seed(101)
bmat <- replicate(500,bootfun(perturb=0.1,boot_res=TRUE))   ## resample residuals
bmat2 <- replicate(500,bootfun(perturb=0.1,boot_res=FALSE)) ## resample observations
## construct envelopes (pointwise percentile bootstrap CIs)
df$lwr <- apply(bmat, 1, quantile, 0.025, na.rm=TRUE)
df$upr <- apply(bmat, 1, quantile, 0.975, na.rm=TRUE)
df$lwr2 <- apply(bmat2, 1, quantile, 0.025, na.rm=TRUE)
df$upr2 <- apply(bmat2, 1, quantile, 0.975, na.rm=TRUE)

现在画图：

ggplot(df, aes(x,y)) +
    geom_point() +
    geom_ribbon(aes(ymin=lwr, ymax=upr), colour=NA, alpha=0.3) +
    geom_ribbon(aes(ymin=lwr2, ymax=upr2), fill="red", colour=NA, alpha=0.3) +
    geom_line(aes(y=fit)) +
    theme_minimal()

粉红色/浅红色区域是观察级引导 CI（可疑）；灰色区域是剩余的引导 CI。

也可以尝试使用 delta 方法，但 (1) 无论如何它比自举方法做出更强大的假设/近似值，并且 (2) 我没时间了。