Skip to content

nls does not accept start values

7 messages · Brian Ripley, Uwe Ligges, PIKAL Petr

PIKAL Petr
# 
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
Uwe Ligges
#
Petr PIKAL wrote:
If you really want to orgnize it that way, why not simpler as in:

form <- list(y ~ 1/(a - x), a = mean(y))
fit <- nls(form[[1]], data.frame(x=x, y=y), start = form[2])


Uwe Ligges
Brian Ripley
# 
It is not very clear what you are trying to do here, and
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 
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.
On Mon, 2 Mar 2009, Petr PIKAL wrote:

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      PIKAL Petr
    

    

      
      #
    

  


  

      
Thank you

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

            

      
      
      
    

        
Error in numericDeriv(form[[3]], names(ind), env) : 
  Missing value or an infinity produced when evaluating the model

however

 nls(fol[[1]], data.frame(x=x, y=y), start=list(a=mean(y)))

gives correct result. Therefore I started think about how to add a 
"better" starting value for some fits as a second part of my formula list 
to define structure like>

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

I wonder If you can push me to better direction.

Thanks again
Best regards
Petr




Uwe Ligges <ligges at statistik.tu-dortmund.de> napsal dne 02.03.2009 
09:41:45:

            

      
      
      
    

        
"list(a=mean(y))")),

            

      
      
      
    

        
start=list(a=mean(y)))

            

      
      
      
    

        
[[1]]

            

      
      
      
    

        
eval,

            

      
      
      
    

        
it

            

      
      
      
    

        

      
      
    

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

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Uwe Ligges
    

    

      
      #
    

  


  

      
Petr PIKAL wrote:

            

      
      
      
    

        
You can make up a list of lists (each containing one formula and its 
starting values) or specify formulas in one list and starting values in 
a corresponding second list.
You need just the corresponding subsetting in your call to nls such as 
in the simple case I suggested already.

Best,
Uwe

            

      
      
      
    

        

      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      PIKAL Petr
    

    

      
      #
    

  


  

      
Hi

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

            

      
      
      
    

        
"list(a=mean(y))")),

            

      
      
      
    

        
You are correct as usually.

            

      
      
      
    

        
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

            

      
      
      
    

        
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

            

      
      
      
    

        
"list(a=mean(y))")),

            

      
      
      
    

        
start=list(a=mean(y)))

            

      
      
      
    

        
[[1]]

            

      
      
      
    

        
eval,

            

      
      
      
    

        
it

            

      
      
      
    

  
  


  


      

    

    
1 day later

    

      
      

        


  

        

      PIKAL Petr
    

    

      
      #
    

  


  

      
Hallo

Here is my "final" solution, but i still wonder if it could not be 
improved somehow.

fit <- modely(fol, data.frame(x=x, y=y1)
fit
or 
fit <- modely(fol, data.frame(x=x, y=y2)
fit

shows that it arrives to some results. However selection starting values 
as means of response variable seems to me suboptimal.
Is it possible to have some list of functions which will for each item in 
formula list  compute starting values for nls? Or is it this mean guess as 
good as any other?

Petr
petr.pikal at precheza.cz


# formulas
fol <- list(y~1/(a-x), y~a*x^2+b*log(x))
set.seed(111)
x<-1:10
y1 <- 1/(0.5-x) +rnorm(length(x))/10
y2 <- 2*x^2 + 3*log(x) + rnorm(length(x))

# here is a function

modely <- function(formula, data, ...){
ll <- length(formula) # no of items in formula list

# initiation of vectors
result2 <- vector("list", ll)
result1 <- rep(NA, ll)

for(i in 1:ll) {

# some start values ***this shall be improved somehow***
n <- length(all.vars(fol[[i]])) - 2
start <- as.list(rep(mean(data[,2]),n))
names(start) <- letters[1:n]

# fit
fit<-try(nls(formula[[i]], data, start))
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)
}

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


Uwe Ligges <ligges at statistik.tu-dortmund.de> napsal dne 02.03.2009 
10:54:16:

            

      
      
      
    

        
function

            

      
      
      
    

        
which

            

      
      
      
    

        
tool

            

      
      
      
    

        
it

            

      
      
      
    

        
list

            

      
      
      
    

        
try

            

      
      
      
    

        
other

            

      
      
      
    

        
get

            

      
      
      
    

        
them

            

      
      
      
    

        
list,