[R] Confidence intervals with nls()

[Ben Bolker]

  My best guess is that the code you were copying
was from an objective (fitted) function that computed
the gradient itself, so the "predict" result had
a "gradient" attribute.  Here's an alternative that
(I think) works OK, but it comes with the usual
warranty (none).

## simulate "data"
O.age = rep(1:12,each=5)
k = 0.2
Linf = 700
t0 = 0
O.length = round(rnorm(length(O.age),
      Linf*(1-exp(-k*O.age-t0)),sd=5))
OS = data.frame(O.age,O.length)

plot(O.length~O.age, data=OS)

## fit model
Oto = nls(O.length~Linf*(1-exp(-k*(O.age-t0))), data=OS,
      start=list(Linf=1000, k=0.1, t0=0.1), trace=TRUE)
      mod <- seq(0, 12)
mod=seq(1,12, by=0.01)
lines(mod, predict(Oto, list(O.age = mod)))

## use emdbook::deltavar to use delta method
##   to compute standard errors
library(emdbook)
se.fit = sqrt(deltavar(
   Linf*(1-exp(-k*mod-t0)),
   meanval=coef(Oto),
   Sigma=vcov(Oto)))

plot(O.length~O.age, data=OS)
matlines(mod, predict(Oto,list(O.age=mod))+
               outer(se.fit,qnorm(c(.5, .025,.975))),type="l",
               lty=1,col=c("black","gray","gray"))

   I notice that I have messed up the VB formula
here -- it should be -k*(mod-t0) -- but you can
correct that