[R] nls() and loop

Martin Maechler maechler at stat.math.ethz.ch
Sat Oct 21 17:58:16 CEST 2017


Dear Vangi,

>>>>> Evangelina Viotto <evangelinaviotto at gmail.com>
>>>>>     on Fri, 20 Oct 2017 11:37:12 -0300 writes:

    > Hello I´m need fitt growth curve with data length-age. I want to evaluate
    > which is the function that best predicts my data, to do so I compare the
    > Akaikes of different models. I'm now need to evaluate if changing the
    > initial values changes the parameters and which do not allow to estimate
    > the model.

    > To do this I use the function nls(); 

good!

    > and I randomize the initial values (real positive number).

Not a very good idea, I'm sorry to say:
You have enough observations to fit such a simple 3-parameter model:

I'm using your data showing you two models, both provided with R
as "self starting models" already.
Self starting means the use the data (and math and linear least
squares) to get good initial (aka "starting") values for the parameters.

The first model, SSasymp() is equivalent yours but more smartly
parametrized: it uses exp(lrc) (!) -- see   help(SSasymp)   in R, 	      
the 2nd model assumes the true curve goes through the origin ( = (0,0)
and hence uses one parameter less.
As we will see, both models fit ok, but the more simple models
may be preferable.

Here is the (self contained) R code, including your data at the beginning:



ANO <- c( 1.65, 1.69, 1.72, 1.72, 1.72, 1.72, 1.73, 2.66 ,2.66, 2.66, 2.66, 2.76, 2.76, 2.76 ,2.76, 2.78, 2.8, 3.65, 3.65 ,3.65, 3.78, 3.78, 5.07, 7.02, 7.1, 7.81, 8.72, 8.74, 8.78, 8.8, 8.8, 8.83, 8.98, 9.1, 9.11, 9.75, 9.82, 9.84, 9.87, 9.87, 10.99, 11.67, 11.8, 11.81, 13.93, 14.83, 15.82, 15.99, 16.87, 16.88, 16.9, 17.68, 17.79, 17.8, 17.8)

SVL <- c(26.11,29.02,41.13,30.96,37.74,29.02,33.38,24.18,34.35,35.8,29.99,42.59,27.57,47.43,46.95,30.47,29.75,35.8,40.65,36.29,34.83,29.02,43.5,75,68,70,67.5,80,77.5,68,68,73.84,72.14,68,64.5,58.5,72.63,78.44,71.17,70.69,77,79,78,68.5,69.72,71.66,77,77,79,76.5,78.5,79,73,80,69.72)

d.SA <- data.frame(SVL, ANO) # creo data frame {but do _not_ call it 'data'}
str(d.SA) ## 55 x 2
summary(d.SA) # to get an idea;  the plot below is more useful

## MM: just do this: it's equivalent to your model (but better parametrized!)
fm1 <- nls(SVL ~ SSasymp(ANO, Asym, lrc, c0), data = d.SA)
summary(fm1)

## Compute nicely spaced predicted values for plotting:
ANO. <- seq(-1/2, 30, by = 1/8)
SVL. <- predict(fm1, newdata = list(ANO = ANO.))

plot(SVL ~ ANO, d.SA, xlim = range(ANO, ANO.), ylim = range(SVL, SVL.))
lines(SVL. ~ ANO., col="red", lwd=2)
abline(h = coef(fm1)[["Asym"]], col = "tomato", lty=2, lwd = 1.5)
abline(h=0, v=0, lty=3, lwd = 1/2)

## Both from summary  (where 'lrc' has large variance) and because of the fit,
## Trying the Asymptotic through the origin ==> 1 parameter less instead :
fmO <- nls(SVL ~ SSasympOrig(ANO, Asym, lrc), data = d.SA)
summary(fmO)## (much better determined)
SVL.O <- predict(fmO, newdata = list(ANO = ANO.))
lines(SVL.O ~ ANO., col="blue", lwd=2)# does go through origin (0,0)
abline(h = coef(fmO)[["Asym"]], col = "skyblue", lty=2, lwd = 1.5)

## look close, I'd probably take the simpler one:
## and classical statistical inference also does not see a significant
## difference between the two fitted models :
anova(fm1, fmO)
## Analysis of Variance Table

## Model 1: SVL ~ SSasymp(ANO, Asym, lrc, c0)
## Model 2: SVL ~ SSasympOrig(ANO, Asym, lrc)
##   Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
## 1     52     2635.1
## 2     53     2702.6 -1 -67.55   1.333 0.2535

---------------

I see that the 10 nice self-starting models that came with nls
already in the  1990's   are not known and probably not
understood by many.

I'm working at making their help pages nicer, notably by
slightly enhancing  the nice model-visualizing plot, you already now
get in R when you run

    example(SSasymp)
or
    example(SSasympOrig)

(but unfortunately, they currently use 'lwd = 0' to draw the asymptote
 which shows fine on a PDF but not on a typical my screen graphics device.)


Martin Maechler
ETH Zurich and R Core team



More information about the R-help mailing list