[R] nls() and loop

Evangelina Viotto evangelinaviotto at gmail.com
Fri Oct 20 16:37:12 CEST 2017


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]]



More information about the R-help mailing list