Reading the data

Data <- read.csv("Alfonsino.csv")
head(Data)
##    Age Length Sex
## 1 10.5   40.0   M
## 2 12.5   40.2   M
## 3 12.5   40.5   F
## 4  8.5   34.3   M
## 5  6.5   31.4   F
## 6  8.5   39.4   M
names(Data)
## [1] "Age"    "Length" "Sex"
Age <- Data$Age
Length <- Data$Length
Sex <- as.numeric(Data$Sex)

Visual presenation of data

mark <- c(19, 2)
color <- c("red","blue")
plot(Length~Age, pch=mark[Sex],
     col=color[Sex], xlim=c(0,25), 
     ylim=c(0,50))
legend(18,20,pch=mark,col=color,
       legend=c("Female","Male"), bty="n")

Model

We will employ a famous growth curve, which the von Bertalanffy growth curve defined as \[ L_t = L_{\infty} (1-exp(-K (t-t0))) \]

growth.VB <- function(t, Linf, K, t0){
  Linf*(1-exp(-K*(t-t0)))
}

Estimation of parameters under an assumption of common parameters between female and male

start <- list(Linf=max(Length), K=0.1, t0=0)

res.000 <- nls(Length~growth.VB(t=Age, Linf, K, t0), start=start)

summary(res.000)
## 
## Formula: Length ~ growth.VB(t = Age, Linf, K, t0)
## 
## Parameters:
##      Estimate Std. Error t value Pr(>|t|)    
## Linf 47.95842    2.91479  16.453  < 2e-16 ***
## K     0.08936    0.02413   3.704  0.00027 ***
## t0   -6.38637    2.14489  -2.977  0.00324 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.564 on 215 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 5.63e-06
confint(res.000)
## Waiting for profiling to be done...
##              2.5%      97.5%
## Linf  44.10249477 58.3766050
## K      0.04465578  0.1360594
## t0   -12.34296242 -3.3029532
# Visual presentation
taxis <- seq(0,25,0.5); tlen <- length(taxis);
pred.000 <- predict(res.000, list(Age=taxis))
plot(Length~Age,xlim=c(0,25),ylim=c(0,50), main="Results under res.000 with prediction CI")
points(taxis, pred.000, type="l")
points(taxis, pred.000-1.96*sigma(res.000), type="l", lty=2)
points(taxis, pred.000+1.96*sigma(res.000), type="l", lty=2)

# Diagnostics 
e <- residuals(res.000)/sigma(res.000)
plot(fitted(res.000), e, xlab="Fitted values", ylab="Standardized residuals")
abline(0,0); abline(-1.96,0,lty=2); abline(1.96,0,lty=2);

Estimation of parameters under an assumption of different parameters between female and male

# Setting a set of initial values for female and male simultaneously
start <- list(Linf=rep(max(Length),2),K=c(0.1,0.1),t0=c(0,0))

# Estimation
res.111 <- nls(Length~growth.VB(t=Age,Linf[Sex], K[Sex], t0[Sex]), start=start)

# Showing the estimation summary
summary(res.111)
## 
## Formula: Length ~ growth.VB(t = Age, Linf[Sex], K[Sex], t0[Sex])
## 
## Parameters:
##       Estimate Std. Error t value Pr(>|t|)    
## Linf1 48.15050    4.52391  10.644  < 2e-16 ***
## Linf2 47.70053    3.75965  12.687  < 2e-16 ***
## K1     0.09797    0.04315   2.270  0.02418 *  
## K2     0.08646    0.02995   2.887  0.00429 ** 
## t01   -5.12851    3.38334  -1.516  0.13106    
## t02   -6.94091    2.80413  -2.475  0.01410 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.567 on 212 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 3.178e-06
# Confidence interval for the growth paramters
confint(res.111)
## Waiting for profiling to be done...
##               2.5%      97.5%
## Linf1  43.09967561 81.7105984
## Linf2  43.18112305 66.9658334
## K1      0.02164258  0.1833475
## K2      0.03031102  0.1453427
## t01   -18.02360053 -0.9526554
## t02   -16.13978541 -3.1310768
# Visual presentation
taxis <- seq(0,25,0.5); tlen <- length(taxis);
pred.female.111 <- predict(res.111, list(Age=taxis, Sex=rep(1,tlen)))
pred.male.111 <- predict(res.111, list(Age=taxis, Sex=rep(2,tlen)) )

plot(Length~Age,pch=mark[Sex],col=color[Sex],xlim=c(0,25),ylim=c(0,50), 
     main="Results under res.111 with prediction CI")
points(taxis, pred.female.111, type="l", col="red", lty=1)
points(taxis, pred.female.111-1.96*sigma(res.111), type="l", col="red", lty=2)
points(taxis, pred.female.111+1.96*sigma(res.111), type="l", col="red", lty=2)
points(taxis, pred.male.111, type="l", col="blue", lty=1)
points(taxis, pred.male.111-1.96*sigma(res.111), type="l", col="blue", lty=2)
points(taxis, pred.male.111+1.96*sigma(res.111), type="l", col="blue", lty=2)
legend(18,20,pch=mark,col=color,lty=c(1,2),legend=c("Female","Male"),bty="n")

# Diagnostics for female
e <- residuals(res.111)/sigma(res.111)
plot(fitted(res.111)[Sex==1], e[Sex==1], col="red", 
     xlab="Fitted values", ylab="Standardized residuals")
abline(0,0); abline(-1.96,0,lty=2); abline(1.96,0,lty=2);

# Diagnostics for male
plot(fitted(res.111)[Sex==2], e[Sex==2], col="blue", 
     xlab="Fitted values", ylab="Standardized residuals")
abline(0,0); abline(-1.96,0,lty=2); abline(1.96,0,lty=2)

Model with res.100

start <- list(Linf=rep(max(Length),2), K=0.1, t0=0)

res.100 <- nls(Length~growth.VB(t=Age,Linf[Sex], K, t0), start=start)

summary(res.100)
## 
## Formula: Length ~ growth.VB(t = Age, Linf[Sex], K, t0)
## 
## Parameters:
##       Estimate Std. Error t value Pr(>|t|)    
## Linf1 48.54689    3.08114  15.756  < 2e-16 ***
## Linf2 47.90648    3.00673  15.933  < 2e-16 ***
## K      0.08754    0.02420   3.618 0.000371 ***
## t0    -6.58214    2.21800  -2.968 0.003343 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.559 on 214 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 7.413e-06

Model with res.010

start <- list(Linf=max(Length), K=c(0.1, 0.1), t0=0)

res.010 <- nls(Length~growth.VB(t=Age,Linf, K[Sex], t0), start=start)

summary(res.010)
## 
## Formula: Length ~ growth.VB(t = Age, Linf, K[Sex], t0)
## 
## Parameters:
##      Estimate Std. Error t value Pr(>|t|)    
## Linf 48.32799    3.11385  15.520  < 2e-16 ***
## K1    0.08744    0.02459   3.555 0.000465 ***
## K2    0.08513    0.02378   3.580 0.000426 ***
## t0   -6.73140    2.25591  -2.984 0.003177 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.561 on 214 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 5.365e-06

Model with res.001

start <- list(Linf=max(Length), K=0.1, t0=c(0,0))

res.001 <- nls(Length~growth.VB(t=Age,Linf, K, t0[Sex]), start=start)

summary(res.001)
## 
## Formula: Length ~ growth.VB(t = Age, Linf, K, t0[Sex])
## 
## Parameters:
##      Estimate Std. Error t value Pr(>|t|)    
## Linf 48.39942    3.15358  15.347  < 2e-16 ***
## K     0.08529    0.02402   3.551 0.000473 ***
## t01  -7.03254    2.32792  -3.021 0.002827 ** 
## t02  -6.67567    2.25953  -2.954 0.003483 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.564 on 214 degrees of freedom
## 
## Number of iterations to convergence: 5 
## Achieved convergence tolerance: 1.406e-06

Model with res.110

start <- list(Linf=rep(max(Length),2), K=rep(0.1,2), t0=0)

res.110 <- nls(Length~growth.VB(t=Age,Linf[Sex], K[Sex], t0), start=start)

summary(res.110)
## 
## Formula: Length ~ growth.VB(t = Age, Linf[Sex], K[Sex], t0)
## 
## Parameters:
##       Estimate Std. Error t value Pr(>|t|)    
## Linf1 49.66398    3.71760  13.359  < 2e-16 ***
## Linf2 46.99869    2.88861  16.270  < 2e-16 ***
## K1     0.08431    0.02384   3.536 0.000498 ***
## K2     0.09291    0.02562   3.627 0.000359 ***
## t0    -6.34966    2.16198  -2.937 0.003679 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.562 on 213 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 3.422e-06

Model with res.011

start <- list(Linf=max(Length), K=rep(0.1,2), t0=rep(0,2))

res.011 <- nls(Length~growth.VB(t=Age,Linf, K[Sex], t0[Sex]), start=start)

summary(res.011)
## 
## Formula: Length ~ growth.VB(t = Age, Linf, K[Sex], t0[Sex])
## 
## Parameters:
##      Estimate Std. Error t value Pr(>|t|)    
## Linf 47.90369    2.88967  16.578  < 2e-16 ***
## K1    0.10029    0.02990   3.354 0.000942 ***
## K2    0.08491    0.02261   3.755 0.000224 ***
## t01  -4.96476    2.45415  -2.023 0.044323 *  
## t02  -7.07769    2.26390  -3.126 0.002017 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.561 on 213 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 4.91e-06

Model with res.101

start <- list(Linf=rep(max(Length),2), K=0.1, t0=rep(0,2))

res.101 <- nls(Length~growth.VB(t=Age,Linf[Sex], K, t0[Sex]), start=start)

summary(res.101)
## 
## Formula: Length ~ growth.VB(t = Age, Linf[Sex], K, t0[Sex])
## 
## Parameters:
##       Estimate Std. Error t value Pr(>|t|)    
## Linf1 48.97747    3.13603  15.618  < 2e-16 ***
## Linf2 47.23314    2.86083  16.510  < 2e-16 ***
## K      0.09044    0.02452   3.688 0.000287 ***
## t01   -5.73684    2.22748  -2.575 0.010687 *  
## t02   -6.58757    2.19411  -3.002 0.002999 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.561 on 213 degrees of freedom
## 
## Number of iterations to convergence: 4 
## Achieved convergence tolerance: 7.627e-06

Comparison of models

AIC(res.000)
## [1] 1034.122
AIC(res.100)
## [1] 1034.308
AIC(res.010)
## [1] 1034.706
AIC(res.001)
## [1] 1035.137
AIC(res.110)
## [1] 1035.71
AIC(res.101)
## [1] 1035.587
AIC(res.011)
## [1] 1035.538
AIC(res.111)
## [1] 1037.531

Conclusion

Based on the comparison of AIC values, the best model is “res.000”, which means that there are no sex-difference in the growth of Alfonsino.