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

Question

我使用nlmeR包和包含的函数构建了几个广义非线性最小二乘模型（指数衰减）。我不简单地使用基函数构建非线性最小二乘模型的原因是因为我希望能够对异方差进行建模以避免变换。我的模型看起来像这样：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()是广义最小二乘法，使用nlmegls()函数运行运行。

现在我想计算置信区间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()对象，因此我查阅了教科书以及以前提出的问题，但似乎都不符合我的需要。我发现的唯一类似的问题是 如何计算 r 中非线性最小二乘的置信区间？。在最上面的答案中，建议使用 或investr::predFit()来构建模型drc::drm()，然后使用常规predict()函数。这些解决方案都没有帮助我gnls()。

我当前的最佳解决方案是使用该函数计算所有三个参数（C、k、A）的 95% 置信区间confint()，然后为置信上限和下限编写两个单独的函数，即一个使用 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，看起来很宽 \xe2\x80\x94 我认为这是因为这种形式的引导无法维持 x 分布的平衡（例如通过重新采样您可能最终无法观察到 x) 的小值。（也有可能我的代码中有一个错误。）

\n

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

\n

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

\n

df$res <- df$y-df$fit\nbootfun <- function(newdata=df, perturb=0, boot_res=FALSE) {\n    start <- coef(mgnls)\n    ## if we start exactly from the previously fitted coefficients we end\n    ## up getting all-identical answers? Not sure what\'s going on here, but\n    ## we can fix it by perturbing the starting conditions slightly\n    if (perturb>0) {\n        start <- start * runif(length(start), 1-perturb, 1+perturb)\n    }\n    if (!boot_res) {\n        ## bootstrap raw data\n        dfboot <- df[sample(nrow(df),size=nrow(df), replace=TRUE),]\n    } else {\n        ## bootstrap residuals\n        dfboot <- transform(df,\n                            y=fit+sample(res, size=nrow(df), replace=TRUE))\n    }\n    bootfit <- try(update(mgnls,\n                      start = start,\n                      data=dfboot),\n                   silent=TRUE)\n    if (inherits(bootfit, "try-error")) return(rep(NA,nrow(newdata)))\n    predict(bootfit,newdata=newdata)\n}\n

\n

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

\n

现在画图：

\n

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

\n

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

\n

\n

尝试一下 delta 方法也不错，但是（1）它比自举法做出了更强的假设/近似值，（2）我没时间了。

\n