On 7/11/2008, at 11:33 PM, Christoph Scherber wrote:
Dear all,
I would like to get standard errors (or confidence intervals) for
*predicted* values from an nls fit.
I have tried to implement code from p.225 in MASS (bootstrapping a
nls fit), but this gives only the
confidence intervals of the parameter estimates, but not an overall
confidence interval for the
predicted value(s):
puro1<-nls(rate~a*conc/(b+conc), data=Puromycin[1:12,], start=list
(a=200, b=1)) #set up nls model
# assume only one predicted value is obtained using predict
(puro1,list(conc=0.02)):
st=cbind(Puromycin[1:12,],fit=predict(puro1,list(conc=0.02)))
puro.bf=function(rs,i){
st$rate=st$fit+rs[i]
coef(nls(rate~a*conc/(b+conc), st, start=coef(puro1)))
}
rs=scale(resid(puro1),scale=F)
(puro.boot=boot(rs,puro.bf,R=100))
boot.ci(puro.boot,index=1,type=c("norm"))
boot.ci$t0
boot.ci$norm
###
How do I have to change the code to get the c.i. for the predicted
value?
First: What you want is ***NOT*** a confidence interval, it is a
prediction
interval.
There are in effect two components --- a ``confidence part'' and a
``predict a new
observation'' part.
The ``confidence part'' is actually the stumbling block.
You are fitting a model
Y = f(x,theta) + E
where the function f() is non-linear in theta. You are probably
willing to assume that
the E_i are i.i.d. Gaussian, mean 0, standard deviation sigma.
The nls fit gives you an estmate of theta and of the standard error
in this estimate.
It also gives you an estimate of sigma.
A 95% confidence level *prediction interval* for a new value of Y at
a given value of x
is given (if the sample size is fairly large) by
f(x,theta-hat) +/- 1.96 * sqrt( (SE-fit)^2 + s^2 )
where s^2 is your (readily available) estimate of sigma^2 and SE-fit
is the estimate
of the standard error of f(x,theta-hat). The problem is that SE-fit
is *not* readily
available (even though the standard errors of theta-hat *are*) since f
() is non-linear
in theta.
As was previously discussed in this thread, the ``usual'' method of
getting at SE-fit
is the so-called ``delta method'' and as Professor Ripley has pointed
out this is in
general insufficiently accurate to be really useful.
Bootstrapping has been suggested as an alternative approach to
getting at SE-fit.
To use boot-strapping for this you must calculate and retain f
(x,theta-hat*) for each bootstrap
sample and then calculate the standard error of the values of f
(x,theta-hat*).
(This is probably somewhat naive --- there may be subtleties here
that I've overlooked
and more sophisticated steps that one can take to address these
subtleties. I am
definitely not an expert on bootstrapping, so look elsewhere for
better advice here.)
If your sample size is ``reasonably large'' then it is ``reasonably
safe'' to assume
that SE-fit is small (negligible) in comparison with s and can be
ignored, i.e. set
to 0. This is the approach that is (almost?) universally adopted in
constructing prediction
intervals in (ARMA) time series analysis.
cheers,
Rolf Turner
######################################################################
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
