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.
______________________________________________
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.
