[R] standard errors for predict.nls?

[R Heberto Ghezzo, Dr]

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.