Hi

--- Susan Imholt <

[hidden email]> wrote:

>

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

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)

pdd.lik.coef<-function(data,i){

d<-data[i,]

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

list(a = 30, b = 36), trace=FALSE,

control=list(maxiter=100))

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

}

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

R=1000)

pdd.boot$t

hope this helps

Marco

______________________________________________

[hidden email] mailing list

https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide!

http://www.R-project.org/posting-guide.html