R help - Oct 2017 - nls() and loop

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(); and I randomize the initial values
(real positive number).  To that I put it inside a function that every time
q is executed it changes the initial parameters and affter and then do a
loop  y and  save a list of the results that interest me in the function.
this problem is does not converge by the initial values, the loop stops and
throws error.
I need to continue and  save initial values with the error that generates
those values


Cheers

Vangi



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)

data<-data.frame (SVL, ANO)# creo data frame
data
> Logiscorri<-function(){
+   a<-runif(1, min=0, max=150)#devuelve 1 al azar dentro de un max y un
min
+   b<-runif(1, min=0, max=100)
+   g<-runif (1, min=0, max=1)
+   d<-runif (1,min=0, max=100)
+
+   ## estimo la curva de distribucion de mis datos
+   caiman<-nls(SVL~DE+(alfa/(1+exp(-gamma*ANO))),
+               data=data,
+               start=list(alfa= a  ,gamma= g, DE= d),
+               control=nls.control(maxiter = 100, warnOnly=TRUE),
+               trace=FALSE)
+   caimansum<-summary(caiman)#ME DA LOS PARAMETROS ESTIMADO, EL NUM DE
ITERACIONES
+   ## analizamos akaike
+   akaike<-AIC(caiman)
+   Bayesiano<-BIC(caiman)
+   alfa<-coef(caiman)[1]
+   beta<-coef(caiman)[2]
+   gamma<- coef(caiman)[3]
+   DE<- coef(caiman)[4]
+   formu<-formula(caiman)
+
+   ValoresIniciales<-c(a, g, d)
+   resultados<-list(formu, caimansum, ValoresIniciales, akaike, Bayesiano)
+   return(resultados)
+ }
> Logiscorri()
[[1]]
SVL ~ DE + (alfa/(1 + exp(-gamma * ANO)))
<environment: 0x16a6b89c>

[[2]]

Formula: SVL ~ DE + (alfa/(1 + exp(-gamma * ANO)))

Parameters:
      Estimate Std. Error t value Pr(>|t|)
alfa  133.0765     6.9537  19.138  < 2e-16 ***
gamma   0.2746     0.0371   7.401 1.13e-09 ***
DE    -54.0467     7.1047  -7.607 5.34e-10 ***
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

Residual standard error: 6.821 on 52 degrees of freedom

Number of iterations to convergence: 30
Achieved convergence tolerance: 4.995e-06


[[3]]
[1] 112.2528283   0.4831461  38.5151401

[[4]]
[1] 372.2001

[[5]]
[1] 380.2294
> resultados<-list()
> resultados
list()
> for(i in 1:10){
+   resultados[i]<- list(Logiscorri())
+ }
Error in chol2inv(object$m$Rmat()) :
  element (2, 2) is zero, so the inverse cannot be computed
In addition: Warning message:
In nls(SVL ~ DE + (alfa/(1 + exp(-gamma * ANO))), data = data, start list(alfa =
a,  :
  singular gradient
> names(resultados)<- 1:10
Error in names(resultados) <- 1:10 :
  'names' attribute [10] must be the same length as the vector
[4]
> parametros<- t(sapply(LogisticoConCorri, "[[",
"par?metros")) #estp lo
que hace es ir item por item de la lista y sacar los par?metros
Error in FUN(X[[i]], ...) : subscript out of bounds
> colnames(parametros)<- c("alfa", "beta",
"gamma", "DE")
Error in dimnames(x) <- dn :
  length of 'dimnames' [2] not equal to array
extent
> akaikefinal<- sapply(LogisticoConCorri, "[[",
"akaike")#esto va item por
item de la lista y saca el akaike
Error in FUN(X[[i]], ...) : subscript out of bounds
> bayesfinal<- sapply(LogisticoConCorri, "[[",
"Bayesiano")
Error in FUN(X[[i]], ...) : subscript out of bounds
>
> --
Bi?l. Evangelina V. Viotto
Laboratorio Ecolog?a Animal
Centro de investigaciones Cient?ficas y de Transferencias de
Tecnolog?a Aplicada a la Producci?n
(CICyTTP-CONICET-UADER)
Diamante, Entre R?os
Argentina

	[[alternative HTML version deleted]]

?tryCatch
-- 
Sent from my phone. Please excuse my brevity.

On October 20, 2017 7:37:12 AM PDT, Evangelina Viotto <evangelinaviotto at
gmail.com> wrote:
>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(); and I randomize the initial values
>(real positive number).  To that I put it inside a function that every
>time
>q is executed it changes the initial parameters and affter and then do
>a
>loop  y and  save a list of the results that interest me in the
>function.
>this problem is does not converge by the initial values, the loop stops
>and
>throws error.
>I need to continue and  save initial values with the error that
>generates
>those values
>
>
>Cheers
>
>Vangi
>
>
>
>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)
>
>data<-data.frame (SVL, ANO)# creo data frame
>data
>> Logiscorri<-function(){
>+   a<-runif(1, min=0, max=150)#devuelve 1 al azar dentro de un max y
>un
>min
>+   b<-runif(1, min=0, max=100)
>+   g<-runif (1, min=0, max=1)
>+   d<-runif (1,min=0, max=100)
>+
>+   ## estimo la curva de distribucion de mis datos
>+   caiman<-nls(SVL~DE+(alfa/(1+exp(-gamma*ANO))),
>+               data=data,
>+               start=list(alfa= a  ,gamma= g, DE= d),
>+               control=nls.control(maxiter = 100, warnOnly=TRUE),
>+               trace=FALSE)
>+   caimansum<-summary(caiman)#ME DA LOS PARAMETROS ESTIMADO, EL NUM DE
>ITERACIONES
>+   ## analizamos akaike
>+   akaike<-AIC(caiman)
>+   Bayesiano<-BIC(caiman)
>+   alfa<-coef(caiman)[1]
>+   beta<-coef(caiman)[2]
>+   gamma<- coef(caiman)[3]
>+   DE<- coef(caiman)[4]
>+   formu<-formula(caiman)
>+
>+   ValoresIniciales<-c(a, g, d)
>+   resultados<-list(formu, caimansum, ValoresIniciales, akaike,
>Bayesiano)
>+   return(resultados)
>+ }
>> Logiscorri()
>[[1]]
>SVL ~ DE + (alfa/(1 + exp(-gamma * ANO)))
><environment: 0x16a6b89c>
>
>[[2]]
>
>Formula: SVL ~ DE + (alfa/(1 + exp(-gamma * ANO)))
>
>Parameters:
>      Estimate Std. Error t value Pr(>|t|)
>alfa  133.0765     6.9537  19.138  < 2e-16 ***
>gamma   0.2746     0.0371   7.401 1.13e-09 ***
>DE    -54.0467     7.1047  -7.607 5.34e-10 ***
>---
>Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
>
>Residual standard error: 6.821 on 52 degrees of freedom
>
>Number of iterations to convergence: 30
>Achieved convergence tolerance: 4.995e-06
>
>
>[[3]]
>[1] 112.2528283   0.4831461  38.5151401
>
>[[4]]
>[1] 372.2001
>
>[[5]]
>[1] 380.2294
>
>> resultados<-list()
>> resultados
>list()
>> for(i in 1:10){
>+   resultados[i]<- list(Logiscorri())
>+ }
>Error in chol2inv(object$m$Rmat()) :
>  element (2, 2) is zero, so the inverse cannot be computed
>In addition: Warning message:
>In nls(SVL ~ DE + (alfa/(1 + exp(-gamma * ANO))), data = data, start
>list(alfa = a,  :
>  singular gradient
>> names(resultados)<- 1:10
>Error in names(resultados) <- 1:10 :
>  'names' attribute [10] must be the same length as the vector [4]
>> parametros<- t(sapply(LogisticoConCorri, "[[",
"par?metros")) #estp
>lo
>que hace es ir item por item de la lista y sacar los par?metros
>Error in FUN(X[[i]], ...) : subscript out of bounds
>> colnames(parametros)<- c("alfa", "beta",
"gamma", "DE")
>Error in dimnames(x) <- dn :
>  length of 'dimnames' [2] not equal to array extent
>> akaikefinal<- sapply(LogisticoConCorri, "[[",
"akaike")#esto va item
>por
>item de la lista y saca el akaike
>Error in FUN(X[[i]], ...) : subscript out of bounds
>> bayesfinal<- sapply(LogisticoConCorri, "[[",
"Bayesiano")
>Error in FUN(X[[i]], ...) : subscript out of bounds
>>
>> --
>Bi?l. Evangelina V. Viotto
>Laboratorio Ecolog?a Animal
>Centro de investigaciones Cient?ficas y de Transferencias de
>Tecnolog?a Aplicada a la Producci?n
>(CICyTTP-CONICET-UADER)
>Diamante, Entre R?os
>Argentina
>
>	[[alternative HTML version deleted]]
>
>______________________________________________
>R-help at 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.

Yes, some form of try() is often needed with nls() to avoid scripts stopping.

You might also find nlxb() from package nlsr more reliable in finding solutions.
It uses analytic derivatives if
available if the model is given as an expression, and a Marquardt stabilized
solver. But do expect it to take more
iterations. The syntax is close, but not perfectly equivalent, to that of nls().

JN

On 2017-10-20 11:20 AM, Jeff Newmiller wrote:
> ?tryCatch
>

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