Hi, I am trying to bootstrap the difference between each parameters among two non linear regression (distributed loss model) as following: # data.frame> Raies[1:10,]Tps SolA Solb 1 0 32.97 35.92 2 0 32.01 31.35 3 1 21.73 22.03 4 1 23.73 18.53 5 2 19.68 18.28 6 2 18.56 16.79 7 3 18.79 15.61 8 3 17.60 13.43 9 4 14.83 12.76 10 4 17.33 14.91 etc... # non linear model (work well) RaiesLossA.nls<-nls(SolA~a/(1+b*Tps)^c,start=c(a=32,b=0.5,c=0.6)) RaiesLossB.nls<-nls(Solb~a/(1+b*Tps)^c,start=c(a=33,b=1.5,c=0.5)) # bootstrap library(boot) dif.param<-function(data,i){ RaiesLossA.nls<-nls(SolA[i]~a/(1+b*Tps[i])^c,start=c(a=31,b=0.5,c=0.6)) RaiesLossB.nls<-nls(Solb[i]~a/(1+b*Tps[i])^c,start=c(a=33,b=1.4,c=0.5)) RaiesLossA.nls$m$getPars()-RaiesLossB.nls$m$getPars() }> myboot<-boot(Raies,dif.param,R=1000)Error in numericDeriv(form[[3]], names(ind), env) : Missing value or an Infinity produced when evaluating the model It seems that the init values (start=) may come not to be suitable while bootstraping. Data can be sent offline to whom wanted to try on the dataset. Any hint welcome! Best regards, Patrick Giraudoux University of Franche-Comt?? Department of Environmental Biology EA3184 af. INRA F-25030 Besan??on Cedex tel.: +33 381 665 745 fax.: +33 381 665 797 http://lbe.univ-fcomte.fr
1. Have you studied the documentation, including working the
examples, for "nls" and "boot"? I don't see
"data" as an argument to
"nls", as I'm used to, and your use of "i" also seems
disconcerting to
me. Something like the following seems to me to be closer to the
standard syntax and therefore safer:
nls(SolA~a/(1+b*Tps)^c,start=c(a=31,b=0.5,c=0.6), data=data[i,])
This may not help with this problem, but it might help you with others
-- possibly even in understanding how R is structured more generally.
2. If this error occurs in only a very few cases and you can
afford to ignore such cases, then you could use "try". I didn't
remember exactly the way to do this, so I performed the following
experiment:
> Er <- try(if(NA)2)
Error in if (NA) 2 : missing value where TRUE/FALSE needed
> class(Er)
[1] "try-error"
This leads me to suggest the following:
dif.param<-function(data,i){
RaiesLossA.nls<-try(nls(SolA[i]~a/(1+b*Tps[i])^c,start=c(a=31,b=0.5,c=0.6)))
RaiesLossB.nls<-try(nls(Solb[i]~a/(1+b*Tps[i])^c,start=c(a=33,b=1.4,c=0.5)))
if((class(RaiesLossA.nls)=="try-error") |
(class(RaiesLossB.nls)=="try-error"))
return("whatever you want in this case")
else RaiesLossA.nls$m$getPars()-RaiesLossB.nls$m$getPars()
}
# UNTRIED
3. What is the nature of your data and the physical context? Can
the deviations from this model reasonably be considered to be
independent normal with constant variance? If yes, then you have a
reasonable model. Alternatively, if all numbers are positive, have you
considered the following:
log(SolA) = ln.a + c*log(1+b*Tps)
If the deviations from this model seem close to being normally
distributed, then you might try using this form with method="plinear".
If you do this, you only need starting values for "b". Whether you
use
"plinear" or not, have you tried getting starting values as follows:
Recall that log(1+x) = x+x^2/2+O(x^2). If b*Tps is usually small, then
the following might work to produce reasonable starting values:
fit0 <- lm(SolA~Tps + I(Tps^2), data=data[i,])
Then coef(fit0) estimates log(a), c*b and c*b^2/2. Thus, b <-
2*coef(fit0)[3]/coef(fit0)[2], etc.
hope this helps. spencer graves
p.s. PLEASE do read the posting guide!
"http://www.R-project.org/posting-guide.html". It may help you answer
questions for yourself, and when it doesn't, it can make it easier for
someone else to understand your problem and respond helpfully.
Patrick Giraudoux wrote:
>Hi,
>
>I am trying to bootstrap the difference between each parameters among two
non linear regression (distributed loss model) as
>following:
>
># data.frame
>
>
>>Raies[1:10,]
>>
>>
> Tps SolA Solb
>1 0 32.97 35.92
>2 0 32.01 31.35
>3 1 21.73 22.03
>4 1 23.73 18.53
>5 2 19.68 18.28
>6 2 18.56 16.79
>7 3 18.79 15.61
>8 3 17.60 13.43
>9 4 14.83 12.76
>10 4 17.33 14.91
>etc...
>
># non linear model (work well)
>
>RaiesLossA.nls<-nls(SolA~a/(1+b*Tps)^c,start=c(a=32,b=0.5,c=0.6))
>RaiesLossB.nls<-nls(Solb~a/(1+b*Tps)^c,start=c(a=33,b=1.5,c=0.5))
>
>
># bootstrap
>library(boot)
>
>dif.param<-function(data,i){
>RaiesLossA.nls<-nls(SolA[i]~a/(1+b*Tps[i])^c,start=c(a=31,b=0.5,c=0.6))
>RaiesLossB.nls<-nls(Solb[i]~a/(1+b*Tps[i])^c,start=c(a=33,b=1.4,c=0.5))
>RaiesLossA.nls$m$getPars()-RaiesLossB.nls$m$getPars()
>}
>
>
>
>> myboot<-boot(Raies,dif.param,R=1000)
>>
>>
>
>Error in numericDeriv(form[[3]], names(ind), env) :
> Missing value or an Infinity produced when evaluating the model
>
>It seems that the init values (start=) may come not to be suitable while
bootstraping. Data can be sent offline to whom wanted to
>try on the dataset.
>
>Any hint welcome!
>
>Best regards,
>
>Patrick Giraudoux
>
>
>University of Franche-Comt??
>Department of Environmental Biology
>EA3184 af. INRA
>F-25030 Besan??on Cedex
>
>tel.: +33 381 665 745
>fax.: +33 381 665 797
>http://lbe.univ-fcomte.fr
>
>______________________________________________
>R-help at stat.math.ethz.ch mailing list
>https://www.stat.math.ethz.ch/mailman/listinfo/r-help
>PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
>
>
Hi,
Just a recap on the trials done based on Spencer Grave, Bervin A Turlach and
Christian Ritz's advise.
On their suggestion, the trouble mentioned has well been turned using the
function try() (using the algorithm "plinear" unstead of
"default" was unsuccessful) in the following way:
library(boot)
dif.param<-function(data,i){
RaiesLossA.nls<-try(nls(SolA~a/(1+b*Tps)^c,start=c(a=31,b=0.5,c=0.6),data[i,]),silent=TRUE)
RaiesLossB.nls<-try(nls(Solb~a/(1+b*Tps)^c,start=c(a=33,b=1.4,c=0.5),data[i,]),silent=TRUE)
if ( (inherits(RaiesLossA.nls,"try-error")) |
(inherits(RaiesLossB.nls,"try-error"))) {return(NA)} else {
return(RaiesLossA.nls$m$getPars()-RaiesLossB.nls$m$getPars())}
}
myboot<-boot(Raies,dif.param,R=1000)
The boot objet "myboot" that one get cannot however be handled, eg
with boot.ci(), because "myboot$t" has NA values. This can be
corrected in this way:
myboot$t<-na.omit(myboot$t)
myboot$R<-length(myboot$t[,1])
boot.ci(myboot,index=1,
type=c("norm","basic","perc", "bca"))
Studentized CI appear to be quite unstable here and were not used.
That's it!
The effect of omitting some bootstrap replicates (the less than 10% that could
not lead to a propper fit) is however questionable:
may this
biase the result in a way? I am not too much worried about it for my current
purpose (crude comparison of parameters), but I
guess that it may be an issue for true statisticians...
Many thanks for the most helpful hints provided,
Patrick Giraudoux
>Hi,
>
>I am trying to bootstrap the difference between each parameters among two
non linear regression (distributed loss model) as
>following:
>
># data.frame
>
>
>>Raies[1:10,]
>>
>>
> Tps SolA Solb
>1 0 32.97 35.92
>2 0 32.01 31.35
>3 1 21.73 22.03
>4 1 23.73 18.53
>5 2 19.68 18.28
>6 2 18.56 16.79
>7 3 18.79 15.61
>8 3 17.60 13.43
>9 4 14.83 12.76
>10 4 17.33 14.91
>etc...
>
># non linear model (work well)
>
>RaiesLossA.nls<-nls(SolA~a/(1+b*Tps)^c,start=c(a=32,b=0.5,c=0.6))
>RaiesLossB.nls<-nls(Solb~a/(1+b*Tps)^c,start=c(a=33,b=1.5,c=0.5))
>
>
># bootstrap
>library(boot)
>
>dif.param<-function(data,i){
>RaiesLossA.nls<-nls(SolA[i]~a/(1+b*Tps[i])^c,start=c(a=31,b=0.5,c=0.6))
>RaiesLossB.nls<-nls(Solb[i]~a/(1+b*Tps[i])^c,start=c(a=33,b=1.4,c=0.5))
>RaiesLossA.nls$m$getPars()-RaiesLossB.nls$m$getPars()
>}
>
>
>
>> myboot<-boot(Raies,dif.param,R=1000)
>>
>>
>
>Error in numericDeriv(form[[3]], names(ind), env) :
> Missing value or an Infinity produced when evaluating the model
>
>It seems that the init values (start=) may come not to be suitable while
bootstraping. Data can be sent offline to whom wanted to
>try on the dataset.
>
>Any hint welcome!
>
>Best regards,
>
>Patrick Giraudoux