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() > resultadoslist()> 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:10Error in names(resultados) <- 1:10 : 'names' attribute [10] must be the same length as the vector [4]> parametros<- t(sapply(LogisticoConCorri, "[[", "par?metros")) #estp loque 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 poritem 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