A Problem Fitting the von Bertalanffy Growth Model with nls()

November 19, 2013
By

(This article was first published on fishR » R, and kindly contributed to R-bloggers)

Recently, a fishR user asked me the following (modified) question:

I have successfully fit the typical von Bertalanffy growth model to all of my fish. However, when I subset the data by sex and attempt to fit the model separately to males and females I receive an ‘infinity’ and ‘singular gradient’ error. Why does this happen and is there anything I can do about it?

The user sent me their data which I use below to confirm the problem stated above and to describe some observations that I have made from fitting many von Bertalanffy growth models (VBGMs).

Preliminaries and Data Manipulation

I begin by loading the FSA package …

library(FSA)

… and the data (note that the working direction would have been set before read.csv()) …

df <- read.csv("BKData.csv",header=TRUE)
str(df)
## 'data.frame': 55 obs. of 12 variables:
## $ Primary_Code : Factor w/ 3 levels "SH091613BK","SH091813BK",..: 3 3 3 3 3 1 1 1 1 1 ...
## $ Location : Factor w/ 1 level "SH": 1 1 1 1 1 1 1 1 1 1 ...
## $ Date : Factor w/ 3 levels "16-Oct-12","16-Sep-13",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ Transect : Factor w/ 6 levels "SHWPBTEF03","SHWPHONT02",..: 1 1 1 1 1 3 3 6 6 6 ...
## $ Replicate : int NA NA NA NA NA 1 1 1 1 1 ...
## $ Species : Factor w/ 1 level "CHC": 1 1 1 1 1 1 1 1 1 1 ...
## $ Length : int 526 227 214 226 508 501 486 732 630 588 ...
## $ Weight : int 1519 86 70 85 1223 1156 993 3655 2326 1528 ...
## $ ID_Code : Factor w/ 50 levels "A1","A10","A11",..: 5 6 7 8 1 6 9 19 21 22 ...
## $ Age : int 7 1 1 1 10 7 10 16 10 9 ...
## $ Sex : Factor w/ 3 levels "","Female","Male": 1 1 1 1 2 2 2 2 2 2 ...
## $ Collection_Technique: Factor w/ 2 levels "BTEF","HONT": 1 1 1 1 1 2 2 2 2 2 ...
levels(df$Sex)
## [1] "" "Female" "Male"

(Note: A blank is a poor way to code for “unknown” sex individuals.)

I separated the data into three data frames corresponding to male, female, and unknown sexed individuals …

df_M <- Subset(df,Sex=="Male")
dim(df_M)
## [1] 29 12
df_F <- Subset(df,Sex=="Female")
dim(df_F)
## [1] 22 12
df_U <- Subset(df,Sex=="")
dim(df_U)
## [1] 4 12

Fitting the Typical Von Bertalanffy Model to All Fish

I illustrate the fitting of the typical VBGM by using the vbFuns() convenience function in the FSA package. This function requires the “name” of the VBGM as its only argument and returns a function that can be used to compute fish length from age given parameter values for that model.

( vbT <- vbFuns("typical") )
## function(t,Linf,K=NULL,t0=NULL) {
## if (length(Linf)==3) {
## K <- Linf[2]
## t0 <- Linf[3]
## Linf <- Linf[1]
## } else if (length(Linf)!=1 | is.null(K) | is.null(t0)) {
## stop("One or more model parameters (Linf, K, t0) are missing or incorrect.",call.=FALSE)
## }
## Linf*(1-exp(-K*(t-t0)))
## }
## <environment: 0x1057e184>

The vbStarts() function is used to find reasonable starting values for a particular version of the VBGM. This function requires a formula containing the length and age data, a data frame in data=, and the “name” of the VBGM parameterization in type=.

( svTall <- vbStarts(Length~Age,data=df,type="typical") )
## $Linf
## [1] 555.7
##
## $K
## [1] 0.6162
##
## $t0
## [1] -0.9446

The VBGM is fit with nls() where a formula is the first argument that has the length variable on the right-hand-side and the model function saved from vbFuns() as the left-hand-side where that function has the age variable as its first argument followed by the names of the parameters as the remaining arguments. The nls() function also requires a data frame in data= and the list of starting values in start=.

fitAll <- nls(Length~vbT(Age,Linf,K,t0),data=df,start=svTall)
summary(fitAll)
## Formula: Length ~ vbT(Age, Linf, K, t0)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## Linf 625.4246 52.1546 11.99 <2e-16 ***
## K 0.1906 0.0628 3.03 0.0038 **
## t0 -1.4570 0.9892 -1.47 0.1468
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 72.8 on 52 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 1.55e-06

A quick fit of the model can be seen with

plot(Length~Age,data=df,pch=16,xlab="Age",ylab="Total Length (mm)")
curve(vbT(x,Linf=coef(fitAll)),from=0,to=19,col="red",lwd=2,add=TRUE)

plot of chunk VBGMfit

Fitting the Typical VBGM Separately to Males and Females

A similar process can be followed for fitting the typical VBGM to females …

svTF <- vbStarts(Length~Age,data=df_F,type="typical")
fitTF <- nls(Length~vbT(Age,Linf,K,t0),data=df_F,start=svTF)
## Error: Missing value or an infinity produced when evaluating the model

… but an error occurs. Similarly for males …

svTM <- vbStarts(Length~Age,data=df_M,type="typical")
fitTM <- nls(Length~vbT(Age,Linf,K,t0),data=df_M,start=svTM)
## Error: singular gradient

So, what happened? A plot separated by sex is instructive …

plot(Length~Age,data=df_F,pch=16,xlab="Age",ylab="Total Length (mm)",
     xlim=c(1,19),ylim=c(200,750))
points(Length~Age,data=df_M,pch=16,col="red")
points(Length~Age,data=df_U,pch=15,col="blue")
legend("bottomright",c("Female","Male","Unknown"),pch=c(16,16,15),
       col=c("black","red","blue"),cex=0.75)

plot of chunk MFdata

From this it is seen that all of the fish less than age-3 were of “unknown” sex. Thus, when the data frame was reduced to just males or just females there were no fish less than age-3. Additionally, no female fish were less than age-6, only two female fish were older than age-12, and no males were older than age-11. In other words, when the data were subsetted, the age range of each subset was substantially narrowed.

Furthermore, over these constrained age ranges, the relationship between length and age for both females and males was largely linear. In other words, there is no distinctive curve towards the origin at “young” ages and no obvious asymptote at “older” ages. In addition, there is considerable variability in lengths within a given age (e.g., the lengths for age-6 range from approximately 350 to 625 mm).

These observations largely explain why the typical von Bertalanffy model does not “converge” properly for these subsets of data. In layman’s terms, it is very hard to fit a curved relationship to data with a largely linear relationship coupled with a high degree of natural (i.e., individual) variability.

An Alternative VBGM Helps

Two of the known problems with the typical parameterization of the VBGM is that the parameters are highly correlated (knowing L_{\infty} is very close to knowing K) and the parameters are often on very different scales L_{\infty} might be well into the 100s where as K is usually between 0 and 1). These problems are exacerbated with the data issues identified in the previous problem. As described in the Von Bertalanffy Growth – Intro Vignette on the fishR webpage, Francis (1988) reparameterized the traditional VBGM to have three parameters that are predicted lengths at three ages where the youngest and oldest age are chosen by the user and the third age is exactly between the two chosen ages. The Francis parameterization has the benefit of having parameters that are less correlated and have similar scales. Thus, the Francis parameterization is more likely to fit with the data issues presented here. On the down-side, the Francis parameterization will not provide estimates of the usual L_{\infty}, K, and t_{0} parameters. However, one can still predict mean lengths-at-age (thus, “growth” can be described) and compare growth between groups.

The Francis parameterization is defined with …

( vbF <- vbFuns("Francis") )
## function(t,L1,L2=NULL,L3=NULL,t1,t3=NULL) {
## if (length(L1)==3) {
## L2 <- L1[2]
## L3 <- L1[3]
## L1 <- L1[1]
## } else if (length(L1)!=1 | is.null(L2) | is.null(L3)) {
## stop("One or more model parameters (L1, L2, L3) are missing or incorrect.",call.=FALSE)
## }
## if (length(t1)==2) {
## t3 <- t1[2]
## t1 <- t1[1]
## } else if (length(t1)!=1 | is.null(t3)) {
## stop("One or more model definitions (t1, t3) are missing or incorrect.",call.=FALSE)
## }
## r <- (L3-L2)/(L2-L1)
## L1+(L3-L1)*((1-r^(2*((t-t1)/(t3-t1))))/(1-r^2))
## }
## <environment: 0x08c5ade8>

… and the “young” and “old” ages that I chose to model with are defined with …

ages <- c(3,12)

The model is fit to the female data (without an error) with …

( svFF <- vbStarts(Length~Age,data=df_F,type="Francis",tFrancis=ages,methEV="poly") )
## $L1
## [1] 295.5
##
## $L2
## [1] 473.7
##
## $L3
## [1] 583
fitFF <- nls(Length~vbF(Age,L1,L2,L3,t1=ages),data=df_F,start=svFF)
summary(fitFF)
## Formula: Length ~ vbF(Age, L1, L2, L3, t1 = ages)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## L1 315.6 102.4 3.08 0.0062 **
## L2 473.3 18.7 25.34 4.2e-16 ***
## L3 567.0 21.8 26.00 2.6e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 72 on 19 degrees of freedom
##
## Number of iterations to convergence: 9
## Achieved convergence tolerance: 5.11e-06

The model is fit to the male data (without an error) with …

svFM <- vbStarts(Length~Age,data=df_M,type="Francis",tFrancis=ages,methEV="poly")
fitFM <- nls(Length~vbF(Age,L1,L2,L3,t1=ages),data=df_M,start=svFM)
summary(fitFM)
## Formula: Length ~ vbF(Age, L1, L2, L3, t1 = ages)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## L1 365.0 30.5 11.9 4.6e-12 ***
## L2 537.1 18.0 29.8 < 2e-16 ***
## L3 616.3 53.7 11.5 1.1e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 73.5 on 26 degrees of freedom
##
## Number of iterations to convergence: 3
## Achieved convergence tolerance: 3.38e-06

The two fits can be visualized with …

plot(Length~Age,data=df_F,pch=16,xlab="Age",ylab="Total Length (mm)",
     xlim=c(1,19),ylim=c(200,750))
points(Length~Age,data=df_M,pch=16,col="red")
legend("bottomright",c("Female","Male"),pch=16,col=c("black","red"),cex=0.75)
curve(vbF(x,L1=coef(fitFF),t1=ages),from=3,to=12,lwd=2,add=TRUE)
curve(vbF(x,L1=coef(fitFM),t1=ages),from=3,to=12,lwd=2,col="red",add=TRUE)

plot of chunk VBGMfit2

Further details of fitting this VBGM and methods for constructing bootstrap confidence intervals and developing models to compare parameters between groups are described in the Von Bertalanffy Growth – Intro Vignette on the fishR webpage.


Filed under: Fisheries Science, R Tagged: Growth, R, von Bertalanffy

To leave a comment for the author, please follow the link and comment on his blog: fishR » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.