[R] nls problem

William Dunlap wdun|@p @end|ng |rom t|bco@com
Fri Apr 3 16:46:27 CEST 2020


If you will be copying the printed coefficients into your function (instead
of just using fitted() or predict()), then use dput(coef(m)) to get them
printed to full precision.

Also, if you regress on pH-7 instead of pH you don't have to worry so much
about the roundoff or cancellation error.  This is akin to what
poly(pH,degree) does.

Bill Dunlap
TIBCO Software
wdunlap tibco.com


On Thu, Apr 2, 2020 at 9:44 PM Troels Ring <tring using gvdnet.dk> wrote:

> Thank you – yes indeed, on suggestion from Rui I noticed that the problem
> was mainly that I took the parameters from “m” instead of coef(m) i.e
>
> a        b        d
>
> -1630.71   457.68   -32.11
>
> Instead of
>
>           a           b           d
>
> -1630.70993   457.67555   -32.11469
>
>
>
> And that alone saves me.
>
> BW
> Troels
>
>
>
> *Fra:* William Dunlap <wdunlap using tibco.com>
> *Sendt:* 2. april 2020 18:50
> *Til:* Troels Ring <tring using gvdnet.dk>
> *Cc:* r-help mailing list <r-help using r-project.org>
> *Emne:* Re: [R] nls problem
>
>
>
> 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=/>
>
>         [[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.
>
>
> 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=smtp&url=/>
>

	[[alternative HTML version deleted]]



More information about the R-help mailing list