Re: [R] sampling and nls formula

From: Marco Geraci <>
Date: Tue 07 Feb 2006 - 12:05:20 EST


> Hello,
> I am trying to bootstrap a function that extracts
> the log-likelihood value and the nls coefficients
> from an nls object. I want to sample my dataset
> (pdd) with replacement and for each sampled dataset,
> I want to run nls and output the nls coefficients
> and the log-likelihood value.
> Code:
> x<-c(1,2,3,4,5,6,7,8,9,10)
> y<-c(10,11,12,15,19,23,26,28,28,30)
> pdd<-data.frame(x,y)
> i<-sample(10, replace=TRUE)

i don't think you need a starting sample of the 'index'. You can omit
> i<-sample(10, replace=TRUE)

> pdd.lik.coef<-function(data,i){
> d<-data[i,]
> pdd.nls<-nls(d$y~(a*(d$x)^2)/(b+(d$x)^2),
> data=pdd, start = list(a = 30, b = 36), trace=FALSE)

there's a programming error here
The argument 'data' of 'nls' must match the one that you pass through 'pdd.lik.coef' or the one that you redefine ('d'). Also, when a function like 'nls' has the argument 'data' you don't need to use '$' in the formula (providing that 'data' contains the variables named in the formula).

> pdd.logLik<-logLik(pdd.nls)
> coeff <- coef(pdd.nls)
> lik.coef <- c(pdd.logLik, coeff)

to gain some speed you might want to avoid assigning the result to 'lik.coef', 'coeff', and 'lik.coef'.

Summarizing, the following code should work. 'should' means that might not work. I tried with few boot samples, and it did work. When I tried with R=1000, the program stopped because 'nls' didn't like a particular sample of the data and it reached maxit=50 (default). I added a 'control=list(maxiter=100)' to 'nls' to your code. I suggest you to
work on that.




     pdd.nls<-nls(y~(a*x^2)/(b+x^2), data=d, start =
list(a = 30, b = 36), trace=FALSE,

     c(logLik(pdd.nls), coef(pdd.nls))

pdd.boot<-boot(data=pdd, statistic=pdd.lik.coef, R=1000)

hope this helps


> _______________________________________________
> Join Excite! -
> The most personalized portal on the Web!
> ______________________________________________
> mailing list
> PLEASE do read the posting guide!
> mailing list PLEASE do read the posting guide! Received on Tue Feb 07 12:13:07 2006

This archive was generated by hypermail 2.1.8 : Wed 08 Feb 2006 - 02:44:37 EST