**fishR » R**, and kindly contributed to R-bloggers)

In my last post, I suggested that the Francis parameterization of the von Bertalanffy growth model may be used in cases where the typical parameterization did not converge (likely due to issues related to highly correlated parameters and data with lengths that are highly variable within ages and the full curvature of the model is not readily apparent because the data are truncated for some reason (e.g., high mortality rates, size-selective gear)). A follow-up question to that post is how to compare parameter estimates between sexes when using the Francis parameterization. This is largely the same as the demonstration for the typical parameterization in the Von Bertalanffy Growth – Intro Vignette on the fishR page except, of course, that the user must write out the more bulky Francis parameterization. This post is a demonstration of the required code (with few comments as most of this is generally described in the Von Bertalanffy Growth – Intro Vignette).

## Thank You

I would like to fish publicly thank Yihui Xie, the author of the knitr package for his quick response to my request for an additional feature to his `knit2wp()`

function which allows one to easily create blog posts such as this on wordpress.com in RStudio using the markdown language. The knitr package is amazing, but Yihui’s attention to his users is even more amazing. Thanks Yihui!

## Preliminaries and Data Manipulation

Begin by loading the **FSA** package, …

library(FSA)

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"`

and remove the four unknown sex individuals …

df1 <- Subset(df,Sex!="") dim(df1)

`## [1] 51 12`

## Model Preparation

Declare all possible models, …

frGen <- Length~L1[Sex]+(L3[Sex]-L1[Sex])*(1-((L3[Sex]-L2[Sex])/(L2[Sex]-L1[Sex]))^(2*(Age-t1)/(t3-t1)))/(1-((L3[Sex]-L2[Sex])/(L2[Sex]-L1[Sex]))^2) fr12 <- Length~L1[Sex]+(L3-L1[Sex])*(1-((L3-L2[Sex])/(L2[Sex]-L1[Sex]))^(2*(Age-t1)/(t3-t1)))/(1-((L3-L2[Sex])/(L2[Sex]-L1[Sex]))^2) fr13 <- Length~L1[Sex]+(L3[Sex]-L1[Sex])*(1-((L3[Sex]-L2)/(L2-L1[Sex]))^(2*(Age-t1)/(t3-t1)))/(1-((L3[Sex]-L2)/(L2-L1[Sex]))^2) fr23 <- Length~L1+(L3[Sex]-L1)*(1-((L3[Sex]-L2[Sex])/(L2[Sex]-L1))^(2*(Age-t1)/(t3-t1)))/(1-((L3[Sex]-L2[Sex])/(L2[Sex]-L1))^2) fr1 <- Length~L1[Sex]+(L3-L1[Sex])*(1-((L3-L2)/(L2-L1[Sex]))^(2*(Age-t1)/(t3-t1)))/(1-((L3-L2)/(L2-L1[Sex]))^2) fr2 <- Length~L1+(L3-L1)*(1-((L3-L2[Sex])/(L2[Sex]-L1))^(2*(Age-t1)/(t3-t1)))/(1-((L3-L2[Sex])/(L2[Sex]-L1))^2) fr3 <- Length~L1+(L3[Sex]-L1)*(1-((L3[Sex]-L2)/(L2-L1))^(2*(Age-t1)/(t3-t1)))/(1-((L3[Sex]-L2)/(L2-L1))^2) frCom <- Length~L1+(L3-L1)*(1-((L3-L2)/(L2-L1))^(2*(Age-t1)/(t3-t1)))/(1-((L3-L2)/(L2-L1))^2)

choose the youngest (`t1`

) and oldest (`t3`

) ages to use in the models, …

t1 <- 5 t3 <- 12

find starting values for the “all parameters in common” model, …

( svCom <- vbStarts(Length~Age,data=df1,type="Francis",tFrancis=c(t1,t3),methEV="poly") )

`## $L1`

## [1] 443.8

##

## $L2

## [1] 542.2

##

## $L3

## [1] 602.7

and expand those starting values for use with each possible model …

svGen <- lapply(svCom,rep,2) sv12 <- mapply(rep,svCom,c(2,2,1)) sv13 <- mapply(rep,svCom,c(2,1,2)) sv23 <- mapply(rep,svCom,c(1,2,2)) sv1 <- mapply(rep,svCom,c(2,1,1)) sv2 <- mapply(rep,svCom,c(1,2,1)) sv3 <- mapply(rep,svCom,c(1,1,2))

## Fit All Models

Fit the general, one parameter in common, two parameters in common, and all parameters in common models declared above to the data …

fitGen <- nls(frGen,data=df1,start=svGen) fit12 <- nls(fr12,data=df1,start=sv12) fit13 <- nls(fr13,data=df1,start=sv13) fit23 <- nls(fr23,data=df1,start=sv23) fit1 <- nls(fr1,data=df1,start=sv1) fit2 <- nls(fr2,data=df1,start=sv2) fit3 <- nls(fr3,data=df1,start=sv3) fitCom <- nls(frCom,data=df1,start=svCom)

## Model Comparisons

Compare each pair of “one parameter in common” models …

anova(fit12,fitGen)

`## Analysis of Variance Table`

##

## Model 1: Length ~ L1[Sex] + (L3 - L1[Sex]) * (1 - ((L3 - L2[Sex])/(L2[Sex] - L1[Sex]))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3 - L2[Sex])/(L2[Sex] - L1[Sex]))^2)

## Model 2: Length ~ L1[Sex] + (L3[Sex] - L1[Sex]) * (1 - ((L3[Sex] - L2[Sex])/(L2[Sex] - L1[Sex]))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3[Sex] - L2[Sex])/(L2[Sex] - L1[Sex]))^2)

## Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)

## 1 46 243467

## 2 45 239038 1 4429 0.83 0.37

anova(fit13,fitGen)

`## Analysis of Variance Table`

##

## Model 1: Length ~ L1[Sex] + (L3[Sex] - L1[Sex]) * (1 - ((L3[Sex] - L2)/(L2 - L1[Sex]))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3[Sex] - L2)/(L2 - L1[Sex]))^2)

## Model 2: Length ~ L1[Sex] + (L3[Sex] - L1[Sex]) * (1 - ((L3[Sex] - L2[Sex])/(L2[Sex] - L1[Sex]))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3[Sex] - L2[Sex])/(L2[Sex] - L1[Sex]))^2)

## Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)

## 1 46 268444

## 2 45 239038 1 29406 5.54 0.023 *

## ---

## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(fit23,fitGen)

`## Analysis of Variance Table`

##

## Model 1: Length ~ L1 + (L3[Sex] - L1) * (1 - ((L3[Sex] - L2[Sex])/(L2[Sex] - L1))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3[Sex] - L2[Sex])/(L2[Sex] - L1))^2)

## Model 2: Length ~ L1[Sex] + (L3[Sex] - L1[Sex]) * (1 - ((L3[Sex] - L2[Sex])/(L2[Sex] - L1[Sex]))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3[Sex] - L2[Sex])/(L2[Sex] - L1[Sex]))^2)

## Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)

## 1 46 248698

## 2 45 239038 1 9660 1.82 0.18

From this it is seen that there is no difference in the or parameters but there may be in the parameter. The model with in common (i.e., `fit12`

) fits slightly better (lower RSS) then the model with in common, so the following tests will compare the two two parameter in common models that also have in common to the model with only in common …

anova(fit1,fit12)

`## Analysis of Variance Table`

##

## Model 1: Length ~ L1[Sex] + (L3 - L1[Sex]) * (1 - ((L3 - L2)/(L2 - L1[Sex]))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3 - L2)/(L2 - L1[Sex]))^2)

## Model 2: Length ~ L1[Sex] + (L3 - L1[Sex]) * (1 - ((L3 - L2[Sex])/(L2[Sex] - L1[Sex]))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3 - L2[Sex])/(L2[Sex] - L1[Sex]))^2)

## Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)

## 1 47 269584

## 2 46 243467 1 26117 4.93 0.031 *

## ---

## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(fit2,fit12)

`## Analysis of Variance Table`

##

## Model 1: Length ~ L1 + (L3 - L1) * (1 - ((L3 - L2[Sex])/(L2[Sex] - L1))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3 - L2[Sex])/(L2[Sex] - L1))^2)

## Model 2: Length ~ L1[Sex] + (L3 - L1[Sex]) * (1 - ((L3 - L2[Sex])/(L2[Sex] - L1[Sex]))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3 - L2[Sex])/(L2[Sex] - L1[Sex]))^2)

## Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)

## 1 47 260197

## 2 46 243467 1 16730 3.16 0.082 .

## ---

## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

These results suggest that differs between sexes and, perhaps, that also differs between sexes depending on the level of that one is using. It is probably reasonable to use an of 0.1 with this type of data because of the high degree of variability in lengths and the likely interest in whether there is even a slight difference between sexes. Of course, should have been set way before this stage (i.e., way before looking at the results). I will continue assuming that was 0.05 and that these results show that only differs between sexes.

Now compare the model with and , but not , in common to the model with no differences between sexes to confirm that differs between the sexes …

anova(fitCom,fit2)

`## Analysis of Variance Table`

##

## Model 1: Length ~ L1 + (L3 - L1) * (1 - ((L3 - L2)/(L2 - L1))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3 - L2)/(L2 - L1))^2)

## Model 2: Length ~ L1 + (L3 - L1) * (1 - ((L3 - L2[Sex])/(L2[Sex] - L1))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3 - L2[Sex])/(L2[Sex] - L1))^2)

## Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)

## 1 48 272966

## 2 47 260197 1 12769 2.31 0.14

This model suggests that all parameters in common is at least as good of a model as the model where differs. This, however, is inconsistent with what was seen above. Perhaps the inconsistency comes from the fact that both and should be allowed to differ. Thus, try comparing the model with separate and parameters to the model with all parameters in common.

anova(fitCom,fit12)

`## Analysis of Variance Table`

##

## Model 1: Length ~ L1 + (L3 - L1) * (1 - ((L3 - L2)/(L2 - L1))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3 - L2)/(L2 - L1))^2)

## Model 2: Length ~ L1[Sex] + (L3 - L1[Sex]) * (1 - ((L3 - L2[Sex])/(L2[Sex] - L1[Sex]))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3 - L2[Sex])/(L2[Sex] - L1[Sex]))^2)

## Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)

## 1 48 272966

## 2 46 243467 2 29499 2.79 0.072 .

```
```

There is some, but not significant evidence, for a difference in the and parameters between the sexes. Thus, this ultimately shows that a single model would best fit both sexes.

The AIC results suggest that the and parameters differ between the sexes …

AIC(fitGen,fit12,fit13,fit23,fit1,fit2,fit3,fitCom)

`## df AIC`

## fitGen 7 589.8

## fit12 6 588.7

## fit13 6 593.7

## fit23 6 589.8

## fit1 5 591.9

## fit2 5 590.1

## fit3 5 592.6

## fitCom 4 590.6

Thus, the AIC results suggest (though not strongly) that the mean length of males and females differs at age-8.5, potentially differ at age-5, but do not differ at age-12. The coefficient results from the best-fit model(s) give an indication of the difference in mean lengths at

summary(fit12)

`##`

## Formula: Length ~ L1[Sex] + (L3 - L1[Sex]) * (1 - ((L3 - L2[Sex])/(L2[Sex] -

## L1[Sex]))^(2 * (Age - t1)/(t3 - t1)))/(1 - ((L3 - L2[Sex])/(L2[Sex] -

## L1[Sex]))^2)

##

## Parameters:

## Estimate Std. Error t value Pr(>|t|)

## L11 389.8 46.0 8.48 5.9e-11 ***

## L12 466.8 19.1 24.47 < 2e-16 ***

## L21 502.6 17.4 28.85 < 2e-16 ***

## L22 548.8 14.3 38.45 < 2e-16 ***

## L3 576.0 19.4 29.73 < 2e-16 ***

## ---

## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##

## Residual standard error: 72.8 on 46 degrees of freedom

##

## Number of iterations to convergence: 12

## Achieved convergence tolerance: 8e-06

The two model fits can be visualized with …

xlbl <- "Age (yrs)" ylbl <- "Total Length (mm)" xlmt <- c(5,19) ylmt <- c(300,750) plot(Length~Age,data=Subset(df1,Sex=="Female"),pch=16,xlab=xlbl,ylab=ylbl, main="No differences",xlim=xlmt,ylim=ylmt) points(Length~Age,data=Subset(df1,Sex=="Male"),pch=16,col="red") legend("bottomright",c("Female","Male"),pch=16,col=c("black","red"),cex=0.75) vbF <- vbFuns("Francis") curve(vbF(x,L1=coef(fitCom),t1=t1,t3=t3),from=5,to=12,lwd=2,col="blue",add=TRUE) plot(Length~Age,data=Subset(df1,Sex=="Female"),pch=16,xlab=xlbl,ylab=ylbl, main="L1 and L2 differ",xlim=xlmt,ylim=ylmt) points(Length~Age,data=Subset(df1,Sex=="Male"),pch=16,col="red") legend("bottomright",c("Female","Male"),pch=16,col=c("black","red"),cex=0.75) vbF <- vbFuns("Francis") curve(vbF(x,L1=coef(fit12),t1=t1,t3=t3),from=5,to=12,lwd=2,add=TRUE) curve(vbF(x,L1=coef(fit12),t1=t1,t3=t3),from=5,to=12,lwd=2,col="red",add=TRUE)

## Notes

- I am not sure what to do with the two “older” females in the sample. It does not seem correct to the fit both models out to age-19 as there are no males older than age-11 and only the two females older than age-12. I decided to simply fit the model over the general ages where both sexes were observed.
- I use the two final plots above to demonstrate the results of the models. If I was preparing these plots for publication I would not constrain the two lines to have the exact same parameter. In other words, I would plot the two curves from the
`fitGen`

model and then note that the parameter was not statistically different.

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

**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...