[R] Comparing Von Bertalanffy Growth Curves

[Ben Bolker]

April Lindeman <aprillindeman <at> yahoo.com> writes:


> I am trying to compare Vbert growth curves from several years of
> fish data. I am following the code provided by:
> http://www.ncfaculty.net/dogle/fishR/gnrlex/VonBertalanffy/VonBertalanffy.pdf.
> Specifically the section on VBGM Comparisons between groups.


> This code is pretty cut and dry. I am able to run it perfectly with
> the "fake" data that is provided. But when I run it with my own data
> I get stuck with this line: fitGen <-
> nls(vbGen,data=LMB,start=svGen)

> I get this error code: Error in numericDeriv(form[[3L]], names(ind), env) :
Missing value or an infinity
> produced when evaluating the model.

> Does anyone know how to fix it? I have no missing values and do not know how
to fix the infinity produced. 

> Here is my data set: structure(list(Age = c(1L, 2L, 3L, 5L, 1L, 2L,
> 3L, 4L, 1L, 2L, 3L, 4L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L,
> 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L), MM =
> c(155.79, 296.37, 325.77, 374.65, 181.46, 258.31, 321.88, 355.6,
> 139.75, 230.72, 319.61, 344.84, 130.92, 236.34, 290.53, 360.33,
> 400.61, 155.33, 240.87, 315.46, 345.05, 378.2, 134.71, 256.66,
> 333.71, 362.99, 381.46, 153.91, 217.21, 287.8, 357.28, 385.62,
> 222.25, 288.93, 294.05, 332.79, 367.39), Year = c(2005L, 2005L,
> 2005L, 2005L, 2006L, 2006L, 2006L, 2006L, 2007L, 2007L, 2007L,
> 2007L, 2008L, 2008L, 2008L, 2008L, 2008L, 2009L, 2009L, 2009L,
> 2009L, 2009L, 2010L, 2010L, 2010L, 2010L, 2010L, 2011L, 2011L,
> 2011L, 2011L, 2011L, 2012L, 2012L, 2012L, 2012L, 2012L)), .Names =
> c("Age", "MM", "Year"), class = "data.frame", row.names = c(NA,
> -37L))

> In case it's helpful here is all my code before this point: 
> str(LMB) 
> MMi=as.integer(MM) 
> Yearf=as.factor(Year) 
> Agei=as.integer(Age) 
> ( svCom <- vbStarts(MMi~Agei,data=LMB)) 
> ( svGen <- lapply(svCom,rep,2) ) 
> vbGen <- MMi~Linf[Yearf]*(1-exp(-K[Yearf]*(Age-t0[Yearf]))+error) 
> fitGen <- nls(vbGen,data=LMB,start=svGen) 


  I believe that the reason you're running into trouble is that
you have a very small number of data points per year.  You're trying
to fit a three-parameter model to at *most* 5 data points per year,
and at least 4 (2005-2007 have only 4 ages, 2008-2012 have 5 ages).
In principle this should be possible, but nls is a little bit
fragile and so it gets in trouble.  I wasn't able to run all your
code because some of the packages you use seem to be unavailable
for the development version of R I'm using ...

  The following solution works for me.  You may be better off
fitting a nonlinear mixed model, but I don't have time to work
out that example at the moment ...

LMB <- structure(list(Age = c(1L, 2L, 3L, 5L, 1L, 2L, 3L, 4L, 1L, 2L,
3L, 4L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L,
1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L), MM = c(155.79, 296.37,
325.77, 374.65, 181.46, 258.31, 321.88, 355.6, 139.75, 230.72, 319.61,
344.84, 130.92, 236.34, 290.53, 360.33, 400.61, 155.33, 240.87,
315.46, 345.05, 378.2, 134.71, 256.66, 333.71, 362.99, 381.46, 153.91,
217.21, 287.8, 357.28, 385.62, 222.25, 288.93, 294.05, 332.79,
367.39), Year = c(2005L, 2005L, 2005L, 2005L, 2006L, 2006L, 2006L,
2006L, 2007L, 2007L, 2007L, 2007L, 2008L, 2008L, 2008L, 2008L, 2008L,
2009L, 2009L, 2009L, 2009L, 2009L, 2010L, 2010L, 2010L, 2010L, 2010L,
2011L, 2011L, 2011L, 2011L, 2011L, 2012L, 2012L, 2012L, 2012L,
2012L)), .Names = c("Age", "MM", "Year"), class = "data.frame",
row.names = c(NA, -37L))

library(ggplot2)
LMB$Yearf <- factor(LMB$Year)
(g0 <- ggplot(LMB,aes(Age,MM))+geom_point(aes(colour=factor(Year))))
(g2 <- g0 + geom_smooth(aes(colour=factor(Year)),se=FALSE,method="auto",span=2))
n1 <- nls(MM~SSasympOff(Age,Asym,lrc,c0),data=LMB)
pframe <- data.frame(Age=seq(1,5,length.out=50))
pframe$MM <- predict(n1,newdata=pframe)
(g1 <- g0+geom_line(data=pframe))
library(nlme)
nyear <- length(unique(LMB$Year))
svec <-  c(rbind(coef(n1),
           matrix(0,nrow=nyear-1,ncol=length(coef(n1)))))
fit2 <- gnls(MM~SSasympOff(Age,Asym,lrc,c0),data=LMB,
     params=Asym+lrc+c0~Yearf,
     start=svec)
printCoefmat(summary(fit2)$tTable)
pframe2 <- expand.grid(Age=seq(1,5,length.out=50),
                       Year=unique(LMB$Year))
pframe2$Yearf <- factor(pframe2$Year)
pframe2$MM <- predict(fit2,newdata=pframe2)
g0 + geom_line(data=pframe2,
               aes(colour=factor(Year)))