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

[Gabor Grothendieck]

With good starting values it often won't need to try anything outside
feasible region.

If after trying better starting values you find that that it is still
attempting to move outsie the feasible region then another approach is
to use the boundary value plus some monotonic function of the distance
from the boundary for points outside the feasible region.   If you
were using fixed large values then that won't work as well as smoother
function.

On Tue, Sep 22, 2009 at 4:50 AM, Keith Jewell <k.jewell at campden.co.uk> wrote:
> Thanks Gabor, but my problem isn't finding reasonable starting parameter
> values, it's preventing nls "giving up" when it tries parameter values
> resulting in NA or Inf.
>
> I know queries are often over-specific and the appropriate response is
> "don't start there", so I'm trying to balance between simplifying to
> essentials, and providing enough background. Here's some more background.
>
> I've had this whole system working for several years in S-plus where it's
> routinely and successfully applied to a variety of models and data sets. I'm
> gradually porting it all into R
>
> I'm using several sigmoidal models; the 4-optimisable-parameter Baranyi
> model I mentioned is just one of them. I have algorithms to find starting
> values for each of the sigmoidal models. However, even with reasonable
> starting values, nls sometimes tries "illegal" parameter values resulting in
> a "crash" of the nls process.
>
> Later I fit a more complicated model in which each of the parameters in the
> sigmoidal is a function of several other variables; these models typically
> have dozens of parameters. I use the "simple" sigmoidal fits to generate
> starting values for the complicated model. Even when the simple sigmoidal
> fits have all converged, nls sometimes tries "illegal" parameter values when
> fitting the complicated model, so (even) better starting values for the
> original sigmoidals wouldn't solve my ultimate problem.
>
> I think I've correctly identified my problem as preventing nls "giving up"
> when it tries parameter values resulting in NA or Inf - but I have been
> wrong before :-}
>
> I've considered modifying the sigmoidals to return extreme numeric values
> instead of NA or Inf; but
> a) it's non-trivial to choose an appropriate extreme
> b) it might "break" numericDeriv
> c) it offends me (!) to return a  (wrong!) value when the right answer is NA
>
> Any suggestions/comments gratefully received.
>
> Keith Jewell
> ===========================
> "Gabor Grothendieck" <ggrothendieck at gmail.com> wrote in message
> news:971536df0909210923r3fd13fb0we72850bf73232d41 at mail.gmail.com...
>
> 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.
>>
>
> ______________________________________________
> 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.
>