# Why is my nonlinear least squares confidence band so wide?

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? ## 1 Answer 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.

• This is very helpful! Thank you. Can you just help me understand why you chose BCa interval? Jun 19 '17 at 15:53
• @CephBirk It usually provides better coverage: stats.stackexchange.com/q/43635/11849 Jun 20 '17 at 6:02