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

[Gabor Grothendieck]

With a small number of parameters just use brute force on grid
to calculate starting values.  See nls2 package.

On Mon, Sep 21, 2009 at 12:17 PM, Keith Jewell <k.jewell at campden.co.uk> wrote:
> 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
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>