[Rd] Seeking opinions on possible change to nls() code

[Martin Maechler]

>>>>> J C Nash 
>>>>>     on Fri, 20 Aug 2021 11:06:25 -0400 writes:

    > In our work on a Google Summer of Code project
    > "Improvements to nls()", the code has proved sufficiently
    > entangled that we have found (so far!)  few
    > straightforward changes that would not break legacy
    > behaviour. One issue that might be fixable is that nls()
    > returns no result if it encounters some computational
    > blockage AFTER it has already found a much better "fit"
    > i.e. set of parameters with smaller sum of squares.  Here
    > is a version of the Tetra example:

    time=c( 1,  2,  3,  4,  6 , 8, 10, 12, 16)
    conc = c( 0.7, 1.2, 1.4, 1.4, 1.1, 0.8, 0.6, 0.5, 0.3)
    NLSdata <- data.frame(time,conc)
    NLSstart <-c(lrc1=-2,lrc2=0.25,A1=150,A2=50) # a starting vector (named!)
    NLSformula <-conc ~ A1*exp(-exp(lrc1)*time)+A2*exp(-exp(lrc2)*time)
    tryit <- try(nls(NLSformula, data=NLSdata, start=NLSstart, trace=TRUE))
    print(tryit)

    > If you run this, tryit does not give information that the
    > sum of squares has been reduced from > 60000 to < 2, as
    > the trace shows.

    > Should we propose that this be changed so the returned
    > object gives the best fit so far, albeit with some form of
    > message or return code to indicate that this is not
    > necessarily a conventional solution? Our concern is that
    > some examples might need to be adjusted slightly, or we
    > might simply add the "try-error" class to the output
    > information in such cases.

    > Comments are welcome, as this is as much an infrastructure
    > matter as a computational one.

Hmm...  many years ago, we had introduced the  'warnOnly=TRUE'
option to nls()  i.e., nls.control()  exactly for such cases,
where people would still like to see the solution:

So,

------------------------------------------------------------------------------
> try2 <- nls(NLSformula, data=NLSdata, start=NLSstart, trace=TRUE, 
              control = nls.control(warnOnly=TRUE))
61215.76    (3.56e+03): par = (-2 0.25 150 50)
2.175672    (2.23e+01): par = (-1.9991 0.3171134 2.618224 -1.366768)
1.621050    (7.14e+00): par = (-1.960475 -2.620293 2.575261 -0.5559918)
Warning message:
In nls(NLSformula, data = NLSdata, start = NLSstart, trace = TRUE,  :
  singular gradient

> try2
Nonlinear regression model
  model: conc ~ A1 * exp(-exp(lrc1) * time) + A2 * exp(-exp(lrc2) * time)
   data: NLSdata
   lrc1    lrc2      A1      A2 
 -22.89   96.43  156.70 -156.68 
 residual sum-of-squares: 218483

Number of iterations till stop: 2 
Achieved convergence tolerance: 7.138
Reason stopped: singular gradient

> coef(try2)
      lrc1       lrc2         A1         A2 
 -22.88540   96.42686  156.69547 -156.68461 


> summary(try2)
Error in chol2inv(object$m$Rmat()) : 
  element (3, 3) is zero, so the inverse cannot be computed
>
------------------------------------------------------------------------------

and similar for  vcov(), of course, where the above error
originates.

{ I think  GSoC (andr other)  students should start by studying and
  exploring relevant help pages before drawing conclusions
  ......
  but yes, I've been born in the last millennium ...
}

;-)

Have a nice weekend!
Martin