-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 That does seem reasonable.
But: the boot package has several fancy (or at least let's just say non-naive) methods for computing confidence intervals that are different from the simple quantile approach. I haven't read the Davison and Hinkley book the "boot" package is based on, don't know how much trouble one can get into with the naive approach (which is what I have always done in the past) -- ??? Hell, I wouldn't even have thought to do the residual-based bootstrapping rather than naive bootstraping. cheers Ben Bolker R Heberto Ghezzo, Dr wrote: > Sorry if I am being too naive but if: > after > puro.boot=boot(rs,puro.bf2,R=1000) > # > # we do > # > qq <- apply(puro.boot$t,2,quantile) > matplot(t(qq),type="l") > points(st$conc,st$rate) > points(st$conc,predict(puro1),pch=2,col="red") > # > isn't that what we want graphically and t(qq) gives us the values > obviously quantile can be used with any set of p's > [EMAIL PROTECTED] > > > -----Original Message----- > From: [EMAIL PROTECTED] on behalf of Ben Bolker > Sent: Tue 11/11/2008 3:09 PM > To: [EMAIL PROTECTED] > Cc: r-help@r-project.org > Subject: Re: [R] standard errors for predict.nls? > > > Here's how far I've gotten. It's a start, but > I can't so far find a way to extract the actual > thing we want -- the bootstrap confidence intervals > on the predictions. > > I apologize for not taking the time to go > read up on the boot library etc. etc. (I *have* > RTFM, but I haven't backtracked to the book > TFM is based on ...) and understand exactly > how it works -- hopefully someone can supply > the missing pieces at the end? > > > p1 <- Puromycin[1:12,] > puro1<-nls(rate~a*conc/(b+conc), data=p1, > start=list(a=200, b=1)) #set up nls model > > > ## data frame for residual bootstrapping > st=cbind(Puromycin[1:12,],fit=predict(puro1)) > > ## prediction concentration vector > cvec <- seq(0,1.2,length=100) > > rs=scale(resid(puro1),scale=FALSE) > library(boot) > puro.bf2=function(rs,i){ > st$rate=st$fit+rs[i] > ## intercept errors > n1 <- try(nls(rate~a*conc/(b+conc), st, start=coef(puro1))) > if (class(n1)=="try-error") r <- rep(NA,length(cvec)) > else r <- predict(n1,newdata=list(conc=cvec)) > r > } > puro.boot=boot(rs,puro.bf2,R=1000) > > ## from here on things get a bit bad -- this > ## doesn't work ... how do I extract bootstrap > ## CIs for each concentration value? > tmpbootfun <- function(j,type="perc") { > b <- boot.ci(puro.boot2,type=type,index=j) > b[(length(b)-1):length(b)] > } > bootvals <- t(sapply(1:100,tmpbootfun)) > > > with(p1,plot(rate~conc,ylim=c(0,200))) > invisible(replicate(200,lines(cvec,puro.bf2(rs,sample(12)), > col=rgb(1,0.7,0.7,alpha=0.1)))) > matlines(cvec,bootvals,lty=2) > > > 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? >> >> >> Many thanks and best wishes, >> Christoph >> >> >> >> >> >> Prof Brian Ripley schrieb: >>> On Mon, 3 Nov 2008, Ben Bolker wrote: >>> >>>> Prof Brian Ripley wrote: >>>>>> Christoph Scherber <Christoph.Scherber <at> agr.uni-goettingen.de> >>>>>> writes: >>>>>> >>>>>>> Dear all, >>>>>>> >>>>>>> Is there a way to retrieve standard errors from nls models? >>>>>>> The help page tells me that arguments >>>>>>> such as se.fit are ignored... >>>>>>> >>>>>>> Many thanks and best wishes >>>>>>> Christoph >>>>> In general using the delta method (which is I guess what you mean, >>>>> local >>>>> linearization via derivatives) is nowhere near accurate enough to be >>>>> useful. That's why it has not been done on several occasions in the >>>>> past. >>>>> If you think it might be, see ?delta.method in package alr3. >>>>> >>>>> I would suggest using simulation/bootsrapping to explore the >>>>> uncertainty. >>>>> There is an example in MASS of doing so (and that illustrates some of >>>>> the pitfalls). >>>> Hmmm. By an example, do you mean an example of using bootstrapping to >>>> explore uncertainty in general, or of using bootstrapping to get >>>> standard errors of predictions from nonlinear regressions? I looked >>>> through my copy of MASS (4th ed.) and found only section 5.7 >>>> (bootstrapping in general) and chapter 8 (nonlinear and smooth >>>> regression, esp. p. 225 "bootstrapping" for getting bootstrap c.i.'s on >>>> parameter estimates). I didn't find anything *specifically* covering >>>> s.e./c.i. for nls predictions, but maybe that's not what you meant. >>> I meant the example on p.225 on bootstrapping a nls fit (and that you >>> needed to bootstrap residuals in some cases). You can use almost >>> identical code to set s.e./c.i. for nls predictions. >>> >>>> And yes, I meant "delta method" rather than "delta function" in my >>>> original post. Oops. >>>> >>>> I might add something quick/dirty/naive to the wiki giving >>>> some examples of delta method/bootstrap approaches ... >>>> >>>> If there is no intention to add confidence interval calculation >>>> to predict.se in the foreseeable future might I suggest that the details >>>> under "Value" as to what "se.fit" will do when it is implemented be >>>> removed? (And perhaps even a statement to the effect [as you say >>>> above] that delta method is considered unreliable?) As written it's a >>>> bit of a tease ... >>> I didn't write that ... and its author might have other opinions. >>> >>>> cheers >>>> Ben Bolker >>>> > > ______________________________________________ > 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. > -----BEGIN PGP SIGNATURE----- Version: GnuPG v1.4.9 (GNU/Linux) Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org iEYEARECAAYFAkka5bsACgkQc5UpGjwzenNXKQCeMApgMjsJwk6KkRPDCKufoEuC XhEAn0D16Ww2tIlCWJIpppa73PHA8YsS =V3mt -----END PGP SIGNATURE----- ______________________________________________ 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.