Re: [Rd] Seeking opinions on possible change to nls() code
> J C Nash
> on Fri, 20 Aug 2021 11:41:26 -0400 writes:
> Thanks Martin. I'd missed the intention of that option,
> but re-reading it now it is obvious.
> FWIW, this problem is quite nasty, and so far I've found
> no method that reveals the underlying dangers well. And
> one of the issues with nonlinear models is that they
> reveal how slippery the concept of inference can be when
> applied to parameters in such models.
> JN
Indeed.
Just for the public (and those reading the archives in the future).
When Doug Bates and his phd student José Pinheiro wrote
"the NLME book" (<==> Recommended R package {nlme}
https://cran.R-project.org/package=nlme )
José C. Pinheiro and Douglas M. Bates
Mixed-Effects Models in S and S-PLUS
Springer-Verlag (January 2000)
DOI: 10.1007/b98882 --> https://link.springer.com/book/10.1007%2Fb98882
They teach quite a bit about non-linear regression much of which
seems not much known nor taught nowadays.
NOTABLY they teach self-starting models, something phantastic,
available in R together with nls() but unfortunately *also* not
much known nor taught!
I have improved the help pages, notably the examples for these,
in the distant past I vaguely remember.
Your present 9-point example can indeed also be solved beautiful
by R's builtin SSbiexp() [Self-starting bi-exponential model]:
NLSdata <- data.frame(
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))
## Once you realize that the above is the "simple" bi-exponential model,
## you should remember SSbiexp(), and then
"everything is easy "
try4 <- nls(conc ~ SSbiexp(time, A1, lrc1, A2, lrc2), data = NLSdata,
trace=TRUE, control=nls.control(warnOnly=TRUE))
## --> converges nicely and starts much better anyway:
## 0.1369091 (2.52e+00): par = (-0.7623289 -2.116174 -2.339856 2.602446)
## 0.01160784 (4.97e-01): par = (-0.1988961 -1.974059 -3.523139 2.565856)
## 0.01016776 (1.35e-01): par = (-0.3653394 -1.897649 -3.547569 2.862685)
## 0.01005199 (3.22e-02): par = (-0.3253514 -1.909544 -3.55429 2.798951)
## 0.01004574 (8.13e-03): par = (-0.336659 -1.904219 -3.559615 2.821439)
## 0.01004534 (2.08e-03): par = (-0.3338447 -1.905399 -3.558815 2.816159)
## 0.01004532 (5.30e-04): par = (-0.3345701 -1.905083 -3.559067 2.817548)
## 0.01004531 (1.36e-04): par = (-0.3343852 -1.905162 -3.559006 2.817195)
## 0.01004531 (3.46e-05): par = (-0.3344325 -1.905142 -3.559022 2.817286)
## 0.01004531 (8.82e-06): par = (-0.3344204 -1.905147 -3.559018 2.817263)
## 0.01004531 (7.90e-06): par = (-3.559018 -0.3344204 2.817263 -1.905147)
## even adding central differences and 'scaleOffset' .. but that's not making
big diff.:
try5 <- nls(conc ~ SSbiexp(time, A1, lrc1, A2, lrc2), data = NLSdata,
trace=TRUE, control=nls.control(warnOnly=TRUE, nDcentral=TRUE,
scaleOffset = 1))
## 0.1369091 (1.43e-01): par = (-0.7623289 -2.116174 -2.339856 2.602446)
##
## 0.01004531(5.43e-06): par = (-3.559006 -0.3343852 2.817195 -1.905162)
fitted(try5)
## [1] 0.6880142 1.2416734 1.3871354 1.3503718 1.1051246 0.8451185 0.6334280
0.4717800
## [9] 0.2604932
all.equal( coef(try4), coef(try5)) # "Mean relative difference: 1.502088e-05"
all.equal(fitted(try4), fitted(try5)) # "Mean relative difference: 2.983784e-06"
## and a nice plot:
plot(NLSdata, ylim = c(0, 1.5), pch=21, bg="red")
abline(h=0, lty=3, col="gray")
lines(NLSdata$time, fitted(try5), lty=2, lwd=1/2, col="orange")
tt <- seq(0, 17, by=1/8)
str(pp <- predict(try5, newdata = list(time = tt)))
## num [1:137] -0.7418 -0.4891 -0.2615 -0.0569 0.1269 ...
## - attr(*, "gradient")= num [1:137, 1:4] 1 0.914 0.836 0.765 0.699 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:4] "A1" "lrc1" "A2" "lrc2"
lines(tt, pp, col=4)
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Seeking opinions on possible change to nls() code
Thanks Martin. I'd missed the intention of that option, but re-reading
it now it is obvious.
FWIW, this problem is quite nasty, and so far I've found no method
that reveals the underlying dangers well. And one of the issues with
nonlinear models is that they reveal how slippery the concept of
inference can be when applied to parameters in such models.
JN
On 2021-08-20 11:35 a.m., Martin Maechler wrote:
>> 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 > 6 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
>lrc1lrc2 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
>
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] Seeking opinions on possible change to nls() code
> 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 > 6 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
lrc1lrc2 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
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] Seeking opinions on possible change to nls() code
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 > 6 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. Best, John Nash __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
