Why is my nonlinear least squares confidence band so wide?

Question

I have the following dataset:

d = structure(list(x = c(23.1073083778966, 22.9362216327734, 24.4504147133069, 
  24.7226685133887, 23.2710364752618, 22.5253827558421, 23.545305003427, 
  23.7449683042789, 24.8139647577093, 22.8804536162757, 24.3948588709677, 
  25.4304112554113, 25.5500819672131, 25.7243410214168, 26.6003943661972, 
  26.0698382492864), y = c(NA, 3.536, 3.867, NA, 4.482, 2.033, 
  NA, NA, 2.912, 3.958, 5.445, 6.973, NA, 5.115, 8.382, 4.438)), .Names = c("x", 
  "y"), class = "data.frame", row.names = 15:30)

and I am trying to plot the best exponential fit with a confidence band. I'm trying to fit the data to the following relationship:

 mod = nls(y ~ a * exp(b * x), data = d, start = list(a = 1, b = 0.05))

This gives me a nice fit:

preds = data.frame(x = seq(22, 27, by = 0.1))
preds$y = predict(mod, newdata = preds)

ggplot(d, aes(x, y)) + 
  geom_point() + 
  geom_line(mapping = aes(x, y), data = preds)

However, I'd like to include a confidence band. I found How do I define a confidence band for a custom (nonlinear) function? which gives a nice demonstration of how to do this. The problem is that my confidence band seems unrealistically large:

pa = propagate::predictNLS(mod, newdata = preds)

preds$lcl <- pa$summary[,5]
preds$ucl <- pa$summary[,6]

When I add on

geom_ribbon(aes(x = x, ymin = lcl, ymax = ucl), alpha = 0.3, data = preds, inherit.aes = FALSE)

I get this:

Why is this band so large?

Answer 1

I do not know what predictNLS does, but this result looks indeed not plausible. I would bootstrap the model residuals to derive confidence intervals.

d <- na.omit(d)
mod = nls(y ~ a * exp(b * x), data = d, start = list(a = 1, b = 0.05))

preds = data.frame(x = seq(22, 27, by = 0.1))
preds$y = predict(mod, newdata = preds)

library(boot)
set.seed(42)
myboot <- boot(d, function(d, ind, mod, preds) {
  d$y <- fitted(mod) + residuals(mod)[ind]
  tryCatch(predict(nls(y ~ a * exp(b * x), data = d, start = list(a = 1, b = 0.05)), newdata = preds),
           error = function(e) preds$x * NA)
}, mod = mod, preds = preds, R = 1e4)

CI <- t(sapply(seq_len(nrow(preds)), function(i) boot.ci(myboot, type = "bca", index = i)$bca[4:5]))
colnames(CI) <- c("lwr", "upr")
preds <- cbind(preds, CI)

library(ggplot2)
ggplot(d, aes(x, y)) + 
  geom_ribbon(aes(ymin = lwr, ymax = upr), data = preds, color = NA, fill = "grey80") +
  geom_line(mapping = aes(x, y), data = preds) +
  geom_point() 

Note that I calculate the CIs pointwise from the bootstrap results, which neglects dependencies.