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)
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")
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)))
}
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);
# 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)
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
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
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
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
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
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
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
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.