[R] Confidence intervals for predicted values in nls

joerg van den hoff j.van_den_hoff at fz-rossendorf.de
Mon Jun 7 13:36:01 CEST 2004

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?
>Cristina Silva
>IPIMAR - Departamento de Recursos Marinhos
>Av. de Brasília, 1449-006 Lisboa
>Tel.: 351 21 3027096
>Fax: 351 21 3015948
>csilva at ipimar.pt
>R-help at stat.math.ethz.ch mailing list
>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)

fifu <- deriv(frml,c("a","b"),function(a,b,x){})

#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 <-

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

ff <- outer(c(1,1),f)
kk <- outer(c(1,1),k1)*c(1-amp,1+amp)

xylim <- par('usr')
xpos <- .1*(xylim[2]-xylim[1])+xylim[1]
ypos <- xylim[4] - .1*(xylim[4]-xylim[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!


More information about the R-help mailing list