[R] unstable results of nlxb fit

[Duncan Murdoch]

As John said, sums of exponentials are hard.  One thing that often helps 
a lot is to use the partially linear structure:  given a and b, you've 
got a linear model to compute A and B.  Now that you're down to two 
nonlinear parameters, you can draw a contour plot of nearby values to 
see how much of a mess you're dealing with.

Duncan Murdoch

On 07/05/2020 9:12 a.m., PIKAL Petr wrote:
> Dear all
> 
> I started to use nlxb instead of nls to get rid of singular gradient error.
> I try to fit double exponential function to my data, but results I obtain
> are strongly dependent on starting values.
> 
> tsmes ~ A*exp(a*plast) + B* exp(b*plast)
> 
> Changing b from 0.1 to 0.01 gives me completely different results. I usually
> check result by a plot but could the result be inspected if it achieved good
> result without plotting?
> 
> Or is there any way how to perform such task?
> 
> Cheers
> Petr
> 
> Below is working example.
> 
>> dput(temp)
> temp <- structure(list(tsmes = c(31, 32, 32, 32, 32, 32, 32, 32, 33,
> 34, 35, 35, 36, 36, 36, 37, 38, 39, 40, 40, 40, 40, 40, 41, 43,
> 44, 44, 44, 46, 47, 47, 47, 47, 48, 49, 51, 51, 51, 52, 53, 54,
> 54, 55, 57, 57, 57, 59, 59, 60, 62, 63, 64, 65, 66, 66, 67, 67,
> 68, 70, 72, 74, 76, 78, 81, 84, 85, 86, 88, 90, 91, 92, 94, 96,
> 97, 99, 100, 102, 104, 106, 109, 112, 115, 119, 123, 127, 133,
> 141, 153, 163, 171), plast = c(50, 51, 52, 52, 53, 53, 53, 54,
> 55, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 64, 64, 65, 65, 66,
> 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 75, 76, 76, 77, 77, 78,
> 78, 79, 80, 81, 82, 83, 84, 85, 85, 86, 86, 87, 88, 88, 89, 90,
> 91, 91, 93, 93, 94, 95, 96, 96, 97, 98, 98, 99, 100, 100, 101,
> 102, 103, 103, 104, 105, 106, 107, 107, 108, 109, 110, 111, 112,
> 112, 113, 113, 114, 115, 116)), row.names = 2411:2500, class = "data.frame")
> 
> library(nlsr)
> 
> fit <- nlxb(tsmes ~ A*exp(a*plast) + B* exp(b*plast), data=temp,
> start=list(A=1, B=15, a=0.025, b=0.01))
> coef(fit)
>             A            B            a            b
> 3.945167e-13 9.772749e+00 2.802274e-01 2.179781e-02
> 
> plot(temp$plast, temp$tsmes, ylim=c(0,200))
> lines(temp$plast, predict(fit, newdata=temp), col="pink", lwd=3)
> ccc <- coef(fit)
> lines(0:120,ccc[1]*exp(ccc[3]*(0:120)))
> lines(0:120,ccc[2]*exp(ccc[4]*(0:120)), lty=3, lwd=2)
> 
> # wrong fit with slightly different b
> fit <- nlxb(tsmes ~ A*exp(a*plast) + B* exp(b*plast), data=temp,
> start=list(A=1, B=15, a=0.025, b=0.1))
> coef(fit)
>             A            B            a            b
> 2911.6448377    6.8320597  -49.1373979    0.0261391
> lines(temp$plast, predict(fit, newdata=temp), col="red", lwd=3)
> ccc <- coef(fit)
> lines(0:120,ccc[1]*exp(ccc[3]*(0:120)), col="red")
> lines(0:120,ccc[2]*exp(ccc[4]*(0:120)), lty=3, lwd=2, col="red")
> 
> 
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>