[R] Confidence intervals for predicted values in nls

[joerg van den hoff]

Cristina Silva wrote:

>Dear all
>
>I have tried to estimate the confidence intervals for predicted values of a
>nonlinear model fitted with nls. The function predict gives the predicted
>values and the lower and upper limits of the prediction, when the class of
>the object is lm or glm. When the object is derived from nls, the function
>predict (or predict.nls) gives only the predicted values. The se.fit and
>interval aguments are just ignored.
>
>Could anybody tell me how to estimate the confidence intervals for the
>predicted values (not the model parameters), using an object of class nls?
>
>Regards,
>
>Cristina
>
>------------------------------------------
>Cristina Silva
>IPIMAR - Departamento de Recursos Marinhos
>Av. de Brasília, 1449-006 Lisboa
>Portugal
>Tel.: 351 21 3027096
>Fax: 351 21 3015948
>csilva at ipimar.pt
>
>______________________________________________
>R-help at stat.math.ethz.ch mailing list
>https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
>  
>
maybe this example helps:

==============cut here===============
#define a model formula (a and b are the parameters, "f" is "x"):
frml <- k1 ~ f*(1-a*exp(-b/f))

#simulate some data:
a0 <- .6
b0 <- 1.2
f  <- seq(0.01,4,length=20)
k1true<- f*(1-a0*exp(-b0/f))
#include some noise
amp <- .1
k1 <- rnorm(k1true,k1true,amp*k1true)

#fit:
fifu <- deriv(frml,c("a","b"),function(a,b,x){})
rr<-nls(k1~fifu(a,b,f),start=list(a=.5,b=1))

#the derivatives and variance/covariance matrix:
#(derivs could be extracted from fifu, too)
dk1.da <- D(frml[[3]],'a')
dk1.db <- D(frml[[3]],'b')
covar <- vcov(rr)

#gaussian error propagation:
a <- coef(rr)['a']
b <- coef(rr)['b']
vark1 <-
       eval(dk1.da)^2*covar[1,1]+
       eval(dk1.db)^2*covar[2,2]+
       2*eval(dk1.da)*eval(dk1.db)*covar[1,2]

errk1 <- sqrt(vark1)
lower.bound <- fitted(rr)-2*errk1      #two sigma !
upper.bound <- fitted(rr)+2*errk1      #dito

plot(f,k1,pch=1)
ff <- outer(c(1,1),f)
kk <- outer(c(1,1),k1)*c(1-amp,1+amp)
matlines(ff,kk,lty=3,col=1)

matlines(f,cbind(k1true,fitted(rr),lower.bound,upper.bound),col=c(1,2,3,3),lty=c(1,1,2,2))
xylim <- par('usr')
xpos <- .1*(xylim[2]-xylim[1])+xylim[1]
ypos <- xylim[4] - .1*(xylim[4]-xylim[3])
legend(xpos,ypos,
      c(     
         'data',
         'true',
         'fit', 
         'confidence'
       ),     
       pch=c(1,-1,-1,-1),
       lty=c(0,1,1,2),
       col=c(1,1,2,3)
)
==============cut here===============
if you put this in a file and source it a few times from within R you'll 
get an impression how often the fit deviates from the 'true' curve more 
than expected from
the shown confidence limits.

I believe this approach is 'nearly' valid as long as gaussian error 
probagation is valid (i.e. only to first order in covar and therefore 
for not too large std. errors, am I right?).
to my simple physicist's mind this should suffice to get 'reasonable' 
(probably, in strict sense, not completely correct?) confidence 
intervals for the fit/the prediction.
If somebody objects, please let me know!


joerg