R help - Aug 2008 - Confidence intervals with nls()

I have data that looks like

O.lengthO.age
176	1
179	1
182	1
...
493	5
494	5
514	5
606	5
462	6
491	6
537	6
553	6
432	7
522	7
625	8
661	8
687	10
704	10
615	12
(truncated)

with a simple VonB growth model from within nls():

plot(O.length~O.age, data=OS)
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.001)
lines(mod, predict(Oto, list(O.age = mod)))

I'm trying to put 95% confidence intervals on the nls() regression with code
that I found in the R-help archive:

se.fit <- sqrt(apply(attr(predict(Oto,list(O.age=mod)),"gradient",
col="blue"),1,
                  function(x) sum(vcov(Oto)*outer(x,x))))
matplot(mod, predict(Oto,list(O.age=mod))+
               outer(se.fit,qnorm(c(.5, .025,.975))),type="l")

Unfortunately, I get the error:

Error in apply(attr(predict(Oto, list(O.age = mod)), "gradient", col =
"blue"),  :
  dim(X) must have a positive length

after the se.fit statement is submitted.  My lack of R programming knowledge
prohibits me from debugging this on my own.  My biggest problem is not knowing
what, exactly, is going on in the se.fit or matplot() statements.

Any ideas as to how I can get this to work?

Thanks,

SR  



Steven H. Ranney
Graduate Research Assistant (Ph.D)
USGS Montana Cooperative Fishery Research Unit
Montana State University
PO Box 173460
Bozeman, MT 59717-3460

phone: (406) 994-6643
fax:   (406) 994-7479


	[[alternative HTML version deleted]]

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