[R] How to use nls when [selfStart] function returns NA or Inf??

[Keith Jewell]

Hi Everyone,

I posted this a couple of weeks ago with no responses. My interface (via 
gmane) seemed a bit flakey at the time, so I'm venturing to repost with some 
additional information.

I'm trying to write selfStart non-linear models for use with nls. In these 
models some combinations of parameter values are illegal; the function value 
is undefined.
That's OK when calling the function directly [e.g.  SSmodel(x, pars...)]; it 
is appropriate to return an appropriate non-value such as NA or Inf.
However, when called from nls [e.g. nls(y~SSmodel(x, pars...), ...)] those 
non-values lead to errors such as (but not limited to):
    Error in numericDeriv(form[[3L]], names(ind), env) :
      Missing value or an infinity produced when evaluating the model
or (if I provide a gradient attribute)
    Error in qr.default(.swts * attr(rhs, "gradient")) :
      NA/NaN/Inf in foreign function call (arg 1)

A toy example demonstrating my problem (legal values of param are >1):
#-----------
SSexample<-selfStart(
 model=function(x, param) x^log(param-1),
 initial = function(mCall, data, LHS){
   val<- 1.001
   names(val) <- mCall[c("param")]
   val
   },
 parameters=c("param")
)
#----------------
nls(y~SSexample(x, par), data=data.frame(x=1:10,y=rnorm(10)))
#---------

(repeat the last line a few times and you'll get the error).

I can't see a way of making nls either
a) stick to legal parameter values (which I'd have trouble specifying 
anyway), or
b) (ideally) accept NA/NaN/Inf as indicating "bad" parameter values, 
equivalent to very large errors in 'y' values

I really do want to use nls rather than a bounded optimisation tool (such as 
optim) because this fits into a much bigger picture predicated on nls.

I'd appreciate any suggestions.

Keith Jewell
==================================
sessionInfo() is given below.

I think the toy example above is enough to demonstrate my problem, but in 
case it is relevant (I don't think it is) here is some more info about the 
models I'm fitting:

I'm fitting sigmoidal models to microbial growth over time. The specific 
model giving me problems at the moment is only one of a whole class of such 
models which I need to work with. I specify it here to illustrate that it is 
not always obvious what the bounds on the parameters are.

SSbaranyiBR94<-selfStart(
model=function(Time, y0, ymax, mu, lambda, m = 1, v = mu)
{
 #+
 # From the paper Baranyi J. & Roberts T. A. (1994). A dynamic approach to 
predicting bacterial growth in food. Int J. of Fd. Micro. 23. 277-294
 # Papers equations (6), (5b), (4b)
 # eq 4b: y(Time) = y0 + mu' * A(Time) - ln(1+(exp(m' * mu' * A(Time)) - 
1)/exp(m' * (ymax-y0)))/m'
 # eq 5b: A(Time) = Time + ln((exp(-v * Time) + q0)/(1+ q0))/v
 # from eq 6: q0 = 1/(exp(mu'*lambda)-1)
 #
 #   m' = m*ln(10)
 #  mu' = mu/ln(10)
 #
 # Cited paper gives defaults m = 1 and v = mu
 #-
  m <- m * log(10)
  mu <- mu/log(10)
 .value <- ymax - log(1 + (exp(m * (ymax - y0)) - 1)/exp(m * mu * (Time +
    log((exp( - v * Time) + (1/(exp(v * lambda) - 1)))/(1 + (1/(exp(v *
    lambda) - 1))))/v)))/m
   .value
}
, initial = function(mCall, data, LHS)
{
##### <snip>
},
parameters=c("y0", "ymax", "mu", "lambda")
)
====================
> sessionInfo()
R version 2.9.1 (2009-06-26)
i386-pc-mingw32

locale:
LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United
Kingdom.1252;LC_MONETARY=English_United
Kingdom.1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252

attached base packages:
[1] stats     graphics  grDevices datasets  tcltk     utils     methods
base

other attached packages:
[1] xlsReadWrite_1.3.3 svSocket_0.9-43    svMisc_0.9-48      TinnR_1.0.3
R2HTML_1.59-1      Hmisc_3.6-1

loaded via a namespace (and not attached):
[1] cluster_1.12.0  grid_2.9.1      lattice_0.17-25 stats4_2.9.1
VGAM_0.7-9