>>>>> 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
