[R] nls problem

J C Nash pro|jcn@@h @end|ng |rom gm@||@com
Thu Apr 2 20:06:01 CEST 2020


Yes, I was waiting to see how long before it would be noticed that this is
not the sort of problem for which nls() is appropriate.

And I'll beat the drum again that nls() uses a simple (and generally
deprecated) forward difference derivative approximation that gets into
trouble a VERY high proportion of applications. Duncan Murdoch and I
published nlsr to overcome several of the weaknesses in nls(). nls()
was very advanced when it first appeared, but I've always felt it was
intended for use by very advanced craftsmen rather than general users.
Indeed, it can do some things nlsr won't do, but usually only if handled
skillfully.

JN


On 2020-04-02 12:50 p.m., William Dunlap via R-help wrote:
> Roundoff/cancelation error: compare the following.  The first is equivalent
> to your function, the last to fitted().
> 
>> with(aedf, t(cbind(1, pH, pH^2) %*% round(coef(m), digits=2)))
>             [,1]      [,2]       [,3]         [,4]       [,5]       [,6]
>     [,7]       [,8]       [,9]     [,10]
> [1,] -0.01635965 0.1076024 0.07477337 -0.007166685 -0.1337865 -0.2851877
> -0.4804638 -0.6865135 -0.9870137 -1.398027
>> with(aedf, t(cbind(1, pH, pH^2) %*% round(coef(m), digits=4)))
>            [,1]        [,2]       [,3]       [,4]       [,5]       [,6]
>   [,7]      [,8]      [,9]     [,10]
> [1,] -0.2196286 -0.09864071 -0.1345213 -0.2180792 -0.3463237 -0.4991861
> -0.6959856 -0.903394 -1.205598 -1.618608
>> with(aedf, t(cbind(1, pH, pH^2) %*% round(coef(m), digits=16)))
>            [,1]        [,2]       [,3]       [,4]       [,5]      [,6]
>  [,7]       [,8]      [,9]     [,10]
> [1,] -0.2208705 -0.09989926 -0.1357969 -0.2193638 -0.3476174 -0.500488
> -0.697296 -0.9047119 -1.206926 -1.619947
> 
> Note that your model is linear and could be fitted with
>    lm(data=aedf, Flux ~ pH + I(pH^2))
>    lm(data=aedf, Flux ~ poly(pH, 2))
> The latter uses a more stable parameterization.
> 
> Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
> 
> 
> On Thu, Apr 2, 2020 at 4:15 AM Troels Ring <tring using gvdnet.dk> wrote:
> 
>> Dear friends - I'm on Win10 with R 6.3.1 and have a very simple problem
>> with
>> nls which apparently gives a decent fit to the parable below, even without
>> starting values. But when I then think I know the meaning of the three
>> parameters a, b, and d it goes very wrong. I guess I am again overlooking
>> something easy but cannot spot it.
>>
>> BW
>> Troels Ring,
>>
>> Aalborg, Denmark
>>
>>
>>
>>
>>
>> aedf <- structure(list(Flux = c(-0.141256, -0.154709, -0.215247,
>> -0.302691,
>>
>>             -0.32287, -0.511211, -0.605381, -0.813901, -1.11659, -1.76906
>>
>> ), pH = c(7.06273, 7.11182, 7.16182, 7.18818, 7.21455, 7.23818,
>>
>>           7.26273, 7.28455, 7.31182, 7.34364)), class = "data.frame",
>> row.names = c(NA,
>>
>>
>> -10L))
>>
>> m <- with(aedf,nls(Flux~a + b*pH + d*pH^2))
>>
>>
>>
>> with(aedf,plot(pH,Flux))
>>
>> with(aedf,lines(pH,fitted(m),lty=1,col="red",lwd=3))
>>
>>
>>
>> m
>>
>> # a        b        d
>>
>> # -1630.70   457.67   -32.11
>>
>>
>>
>> fitted(m)
>>
>> # 1] -0.22087053 -0.09989926 -0.13579690 -0.21936385 -0.34761742
>> -0.50048799
>>
>> # [7] -0.69729602 -0.90471194 -1.20692552 -1.61994657
>>
>>
>>
>> FPG <- function(pH) -1630.70 + 457.67*pH -32.11*pH^2
>>
>>
>>
>> FPG(aedf$pH)
>>
>> # [1] -0.016359649  0.107602395  0.074773375 -0.007166685 -0.133786467
>>
>> # [6] -0.285187665 -0.480463769 -0.686513537 -0.987013685 -1.398026917
>>
>>
>>
>> # So why aren't fitted(m) and FPG(aedf$pH) not closer ("equal")?
>>
>>
>>
>>
>> This email has been scanned by BullGuard antivirus protection.
>> For more info visit www.bullguard.com
>> <
>> http://www.bullguard.com/tracking.aspx?affiliate=bullguard&buyaffiliate=smt
>> p&url=/
>> <http://www.bullguard.com/tracking.aspx?affiliate=bullguard&buyaffiliate=smtp&url=/>>
>>
>>
>>         [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> 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.
>>
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> 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.
>



More information about the R-help mailing list