[R] nls does not accept start values

[Petr PIKAL]

Hi

Prof Brian Ripley <ripley at stats.ox.ac.uk> napsal dne 02.03.2009 10:24:52:

> It is not very clear what you are trying to do here, and
> 
> > form <- structure(list(a = list(quote(y ~ 1/(a - x)), 
"list(a=mean(y))")),
> > .Names = "a")
> 
> is using a historic anomaly (see the help page).
> 
> I am gussing you want to give nls an object containing a formula and 
> an expression for the starting value.  It seems you are re-inventing 

You are correct as usually. 

> self-starting nls models: see ?selfStart and MASS$ ca p. 216.
> One way to use them in your example is
> 
> mod <- selfStart(~ 1/(a - x), function(mCall, data, LHS) {
>      structure(mean(eval(LHS, data)), names="a")
> }, "a")
> 
> nls(y ~ mod(x, a))
> 
> But if you want to follow ypur route, youer starting values would be 
> better to be a list that you evaluate in an appropriate context 
> (which y is this supposed to be?).  nls() knows where it will find 
> variables, but it is not so easy for you to replicate its logic 
> without access to its evaluation frames.

It was simplified version of my problem. I want to elaborate a function 
which can take predefined list of formulas, some data and evaluate which 
formulas can fit the data. I was inspired by some article in Chemical 
engineering in which some guy used excel solver for such task. I was 
curious if I can do it in R too. I am not sure if nls is appropriate tool 
for such task but I had to start somewhere.

Here is a function which takes list of formulas and data and gives a 
result for each formula.

modely <- function(formula, data, ...){
ll <- length(formula)   #no of items in formula list
result2 <- vector("list", ll) #prepare results
result1 <- rep(NA, ll)
for(i in 1:ll) {
fit<-try(nls(formula[[i]], data))
if( class(fit)=="try-error") result1[i] <- NA  else result1[i] <- 
sum(resid(fit)^2)
if( class(fit)=="try-error") result2[[i]] <- NA  else result2[[i]] <- 
coef(fit)
}

ooo<-order(result1) #order results according to residual sum

#combine results into one list together with functions used

result <- mapply(c, "sq.resid" = result1, result2) 
names(result) <- as.character(formula)
# output
result[ooo]
}

# data
x <-1:10
y <-1/(.5-x)+rnorm(10)/100

# list of formulas
fol <- structure(list(a = y ~ 1/(a - x), b = y ~ a * x^2 + b * log(x), 
    c = y ~ x^a), .Names = c("a", "b", "c"))

modely(fol, data.frame(x=x, y=y)

does not use "correct" model because when using default start values it 
results in

> nls(fol[[1]], data.frame(x=x, y=y))
Error in numericDeriv(form[[3]], names(ind), env) : 
  Missing value or an infinity produced when evaluating the model

I tried to establish such structure to get more appropriate starting 
values

list(a= list(formula1, start.formula1), b=list(formula2, start.formula2), 
....)

But did not manage yet to get correct syntax for let say mean of response 
values. I try to look more closely what I can achieve with selfStart

Thank you again

Best regards
Petr


> 
> On Mon, 2 Mar 2009, Petr PIKAL wrote:
> 
> > Hi to all
> >
> > OK as I did not get any response and I really need some insight I try
> > again with different subject line
> >
> > I have troubles with correct evaluating/structure of nls input
> >
> > Here is an example
> >
> > # data
> > x <-1:10
> > y <-1/(.5-x)+rnorm(10)/100
> >
> > # formula list
> > form <- structure(list(a = list(quote(y ~ 1/(a - x)), 
"list(a=mean(y))")),
> > .Names = "a")
> >
> > # This gives me an error due to not suitable default starting value
> >
> > fit <- nls(form [[1]] [[1]], data.frame(x=x, y=y))
> >
> > # This works and gives me a result
> >
> > fit <- nls(form [[1]] [[1]], data.frame(x=x, y=y), 
start=list(a=mean(y)))
> >
> > *** How to organise list "form" and call to nls to enable to use other
> > then default starting values***.
> >
> > I thought about something like
> >
> > fit <- nls(form [[1]] [[1]], data.frame(x=x, y=y), start=get(form 
[[1]]
> > [[2]]))
> >            ^^^^^^^^^^^^^^^^^^^
> > but this gives me an error so it is not correct syntax. (BTW I tried 
eval,
> > assign, sustitute, evalq and maybe some other options but did not get 
it
> > right.
> >
> > I know I can put starting values interactively but what if I want them
> > computed by some easy way which is specified by second part of a list,
> > like in above example.
> >
> > If it matters
> > WXP,  R2.9.0 devel.
> >
> > Regards
> > Petr
> >
> > petr.pikal at precheza.cz
> 
> 
> -- 
> Brian D. Ripley,                  ripley at stats.ox.ac.uk
> Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> University of Oxford,             Tel:  +44 1865 272861 (self)
> 1 South Parks Road,                     +44 1865 272866 (PA)
> Oxford OX1 3TG, UK                Fax:  +44 1865 272595