[R] standard errors for predict.nls?

[Ben Bolker]

-----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
> Heberto.Ghezzo at McGill.CA
> 
> 
> -----Original Message-----
> From: r-help-bounces at r-project.org on behalf of Ben Bolker
> Sent: Tue 11/11/2008 3:09 PM
> To: Christoph.Scherber at agr.uni-goettingen.de
> Cc: r-help at 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 at 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-----