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

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

我使用 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,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),data = df)
plot(mgnls) # more homogenous

# calculate predicted values 
new <- data.frame(x = c(1:100,1:100),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,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,fit,colour = group)) +
  theme_minimal()

这是结果图，看起来很整洁。但是，当我对原始数据应用相同的方法时，尽管两组的数据分布相似，但其中一个置信区间出现了奇怪的凸起：

知道为什么会发生这种情况吗？是否有可能获得更平滑的置信区间？

解决方法

我实施了一个引导解决方案。我最初做了标准的非参数自举，它重新采样观察，但这产生了 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,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,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,0.975,na.rm=TRUE)
df$lwr2 <- apply(bmat2,na.rm=TRUE)
df$upr2 <- apply(bmat2,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",alpha=0.3) +
    geom_line(aes(y=fit)) +
    theme_minimal()

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

尝试 delta 方法也不错，但是 (1) 无论如何，它比引导法做出了更强的假设/近似，以及 (2) 我没时间了。