R - Confidence bands for exponential model (nls) in basic graphics

Question

I'm trying to plot a exponential curve (nls object), and its confidence bands. I could easily did in ggplot following the Ben Bolker reply in this post. But I'd like to plot it in the basic graphics style, (also with the shaped polygon)

df <- 
structure(list(x = c(0.53, 0.2, 0.25, 0.36, 0.46, 0.5, 0.14, 
0.42, 0.53, 0.59, 0.58, 0.54, 0.2, 0.25, 0.37, 0.47, 0.5, 0.14, 
0.42, 0.53, 0.59, 0.58, 0.5, 0.16, 0.21, 0.33, 0.43, 0.46, 0.1, 
0.38, 0.49, 0.55, 0.54), 
y = c(63, 10, 15, 26, 34, 32, 16, 31,26, 37, 50, 37, 7, 22, 13, 
21, 43, 22, 41, 43, 26, 53, 45, 7, 12, 25, 23, 31, 19, 
37, 24, 50, 40)), 
.Names = c("x", "y"), row.names = c(NA, -33L), class = "data.frame")

m0 <- nls(y~a*exp(b*x), df, start=list(a= 5, b=0.04))
summary(m0)

coef(m0)
#   a        b 
#9.399141 2.675083 

df$pred <- predict(m0)
library("ggplot2"); theme_set(theme_bw())
g0 <- ggplot(df,aes(x,y))+geom_point()+
        geom_smooth(method="glm",family=gaussian(link="log"))+
        scale_colour_discrete(guide="none")

Thanks in advance!

Answer 1

This seems more of a question about statistics than R. It's very important that you understand where the "confidence interval" comes from. There are many ways of constructing one.

For the purposes of drawing a shaded area plot in R, I'm going to assume that we can add/subtract 2 "standard errors" from the nls fitted values to produce the plot. This procedure should be checked.

df <- 
  structure(list(x = c(0.53, 0.2, 0.25, 0.36, 0.46, 0.5, 0.14, 
                       0.42, 0.53, 0.59, 0.58, 0.54, 0.2, 0.25, 0.37, 0.47, 0.5, 0.14, 
                       0.42, 0.53, 0.59, 0.58, 0.5, 0.16, 0.21, 0.33, 0.43, 0.46, 0.1, 
                       0.38, 0.49, 0.55, 0.54), 
                 y = c(63, 10, 15, 26, 34, 32, 16, 31,26, 37, 50, 37, 7, 22, 13, 
                       21, 43, 22, 41, 43, 26, 53, 45, 7, 12, 25, 23, 31, 19, 
                       37, 24, 50, 40)), 
            .Names = c("x", "y"), row.names = c(NA, -33L), class = "data.frame")

m0 <- nls(y~a*exp(b*x), df, start=list(a= 5, b=0.04))
df$pred <- predict(m0)
se = summary(m0)$sigma
ci = outer(df$pred, c(outer(se, c(-1,1), '*'))*1.96, '+')
ii = order(df$x)
# typical plot with confidence interval
with(df[ii,], plot(x, pred, ylim=range(ci), type='l'))
matlines(df[ii,'x'], ci[ii,], lty=2, col=1)
# shaded area plot
low = ci[ii,1]; high = ci[ii,2]; base = df[ii,'x']
polygon(c(base,rev(base)), c(low,rev(high)), col='grey')
with(df[ii,], lines(x, pred, col='blue'))
with(df, points(x, y))

But I think the following plot is much nicer: