[R] Bootstrapping gnls models

[Christoph Scherber]

Dear all,

Here comes a reproducible example, with the original data added.

The error message when running boot() is:
Error in gnls(response.variable ~ a * LD/(b + LD), params <- list(a +  :
   Step halving factor reduced below minimum in NLS step

boot() in this case only seems to work for very small values of R (say, within [1..5]), which is of 
course not desirable. Maybe this is due to problems with model convergence.

I would very much appreciate any help.

###
LD=c(4.087462841,2.321928095,4.087462841,1,1,4.087462841,2.321928095,1,1,2.321928095,4.087462841,2.321928095,5.930737338,2.321928095,5.930737338,1,1,2.321928095,2.321928095,4.087462841,1,1,2.321928095,4.087462841,4.087462841,1,2.321928095,1,4.087462841,2.321928095,1,2.321928095,5.930737338,4.087462841,1,4.087462841,2.321928095,4.087462841,5.930737338,4.087462841,1,2.321928095,2.321928095,1,2.321928095,1,1,4.087462841,4.087462841,2.321928095)
L=c(1,1,0,1,0,0,1,0,0,0,1,0,1,1,1,0,0,1,0,0,0,1,1,1,1,0,1,0,0,0,1,0,1,1,0,1,1,1,1,0,0,1,1,1,1,0,0,1,1,0)

response.variable=c(0.335737179487179,0.391025641025641,0.391826923076923,0.487980769230769,0.294070512820513,0.507211538461538,0.395833333333333,0.0825320512820513,0.443108974358974,0.290064102564103,0.59775641025641,0.514423076923077,0.65625,0.193910256410256,0.447916666666667,0,0.260416666666667,0.407852564102564,0.44150641025641,0.596153846153846,0.0600961538461538,0.21875,0.256410256410256,0.364583333333333,0.435096153846154,0.0793269230769231,0.249198717948718,0.0304487179487179,0.230769230769231,0.485576923076923,0.684294871794872,0.0737179487179487,0.490384615384615,0.599358974358974,0.215544871794872,0.219551282051282,0.602564102564103,0.907852564102564,0.391025641025641,0.43349358974359,0.0384615384615385,0.337339743589744,0.502403846153846,0.405448717948718,1,0.362980769230769,0.116185897435897,0.459134615384615,0.661057692307692,0.0769230769230769)

a=data.frame(LD,L,response.variable)

model=gnls(model = response.variable ~ a * LD/(b + LD), data = a,
     params = list(a + b ~ L), start = c(1, 1, 1, 1))

df<-cbind(a,fit<-as.numeric(predict(model,list(LD=1,L=0.5))))
rs<-scale(resid(model),scale=F)

model.bootfunc<-function(rs,i){
	df$response.variable<-df$fit+rs[i]
	as.numeric(predict(gnls(response.variable ~ a * LD/(b + LD),
       params <- list(a + b ~ L), start = coef(model), data=df)))
}

(model.boot<-boot(rs,model.bootfunc,R=100))
booted<-boot.ci(model.boot,index=1,type=c("norm"))
booted$t0
booted$normal

###

Best wishes,
Christoph.


Prof Brian Ripley schrieb:
> On Fri, 7 Nov 2008, Christoph Scherber wrote:
> 
>> Dear all,
>>
>> I am trying to bootstrap predictions from gnls models using the 
>> following code:
>>
>> # a is the dataframe with which I am working; it contains the variables
>> # response.variable,LD,L,G,P and F
> 
> And without it your code is not reproducible.
> 
>>
>> ###
>>
>> model=gnls(response.variable ~ a * LD/(b + LD),
>>        params = list(a + b ~ L), start = c(1,1,1,1), data=a)
>>
>> df=cbind(a,fit=predict(model,list(LD=1,L=0.5,G=0.5,P=0.46,F=2.2)))
>> model.bootfunc=function(rs,i){
>>     df$response.variable=df$fit+rs[i]
>>     as.numeric(predict(gnls(response.variable ~ a * LD/(b + LD),
>>        params = list(a + b ~ L), start = coef(model), data=df)))
>> }
>>
>> rs=scale(resid(model),scale=F)
>> (model.boot=boot(rs,model.bootfunc,R=1))
>> booted=boot.ci(model.boot,index=1,type=c("norm","basic","perc","bca"))
> 
> Do please try to make your code readable, using spaces and <- for 
> assignments.  I would have spotted the problem much sooner with legible, 
> reproducible code.
> 
>> ###
>>
>> The problem is that this code yields "NA" for the s.e. of the 
>> bootstrap statistics:
>>
>> Bootstrap Statistics :
>>      original       bias    std. error
>> t1*  0.1651658 -0.020663364          NA
>> t2*  0.1669592 -0.021759335          NA
>> t3*  0.1676765 -0.001858686          NA
>> t4*  0.1726982 -0.025321349          NA
>> t5*  0.1658092  0.024721214          NA
>>
>>
>> And hence the boot.ci function and others don?t work.
>>
>> Does anyone have an idea on that?
> 
> Yes: how can you estimate standard errors from a single sample (you set 
> R=1)?
> 
>> Many thanks and best wishes
>> Christoph
>>
>>
>>
>> -- 
>> Dr. rer.nat. Christoph Scherber
>> University of Goettingen
>> DNPW, Agroecology
>> Waldweg 26
>> D-37073 Goettingen
>> Germany
>>
>> phone +49 (0)551 39 8807
>> fax   +49 (0)551 39 8806
>>
>> Homepage http://www.gwdg.de/~cscherb1
>>
>> ______________________________________________
>> 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.
>>
>