#growth models with corvina. Derek Bolser

#set working directory
setwd('C:/Users/ROV/Documents/R/corvina_growth')

#load packages
library(FSA)
library(AICcmodavg)
library(minpack.lm)
library(nlstools)

#import data

#name data
data1=oneyrcorvinagrowthdata
tl=data1$TL
age=data1$Age

#name functions from FSA
l1 <- logisticFuns()
g1 <- GompertzFuns()
vb1 <- vbFuns()

#gompertz
svG1 <- list(Linf=1000,gi=0.5,ti=1.5)
fitG1 <- nlsLM(tl~g1(age,Linf,gi,ti),data=data1,start=svG1,trace=TRUE)
bootG1 <- nlsBoot(fitG1)
warnings()
cbind(coef(fitG1),confint(bootG1))
gvals=coef(fitG1)

#logistic
svL1 <- list(Linf=1000,gninf=0.5,ti=2)
fitL1 <- nlsLM(tl~l1(age,Linf,gninf,ti),data=data1,start=svL1,trace=TRUE)
bootL1 <- nlsBoot(fitL1)
warnings()
cbind(coef(fitL1),confint(bootL1))
lvals=coef(fitL1)
  
#von bertalanffy
svvb <- list(Linf=1000,K=0.2,t0=-0.1)
fitvb <- nlsLM(tl~vb1(age,Linf,K,t0),data=data1,start=svvb,trace=TRUE)
bootvb <- nlsBoot(fitvb)
warnings()
cbind(coef(fitvb),confint(bootvb))
vbvals=coef(fitvb)

#schnute
svs1 <- list(a=0.3, b=1)
fitS1 <- nlsLM(tl~Schnute(age,case=1,t1=1,t3=8,L1=141,L3=1013,a,b),data=data1,start=svs1,trace=TRUE)
bootS1 <- nlsBoot(fitS1)
warnings()
cbind(coef(fitS1),confint(bootS1))
svals=coef(fitS1)

#schnute-richards (1990); FSA does not contain this model, so inputted manually.
svsr1 <- list(a=800, b=0.3, c=2.5, d=0.002, f=0.9)
fitSR1 <- nlsLM(tl~(a*(1+b*exp(-c*age^d))^(1/f)),data=data1,start=svsr1,control=nls.lm.control(maxiter=500),trace=TRUE)
bootSR1 <- nlsBoot(fitSR1)
warnings()
cbind(coef(fitSR1),confint(bootSR1))
srvals=coef(fitSR1)

#create AIC table
aictab(list(fitvb,fitL1,fitG1,fitS1,fitSR1),c("von Bertalanffy","logistic","Gompertz","Schnute","Schnute-Richards"))

#obtain BIC values
BIC(fitSR1,fitL1,fitG1,fitvb,fitS1)

#view parameter ests
gvals
lvals
vbvals
svals
srvals

#plot
plot(tl~age,data=data1,pch=19,col=rgb(0,0,0,1/4),
     xlab="Age (years)",ylab="Total Length (mm)",xlim=c(0,8),ylim=c(0,1100))
curve(vb1(x,unlist(vbvals)),from=-1,to=9,add=TRUE,lwd=2,col="green")
curve(g1(x,unlist(gvals)),from=-1,to=9,add=TRUE,lwd=2,col="blue")
curve(l1(x,unlist(lvals)),from=-1,to=9,add=TRUE,lwd=2,col="yellow")
curve(Schnute(x,case=1,t1=1,t3=8,L1=141,L3=1013,a=-0.2019429,b=2.5196610),from=-1,to=9,add=TRUE,lwd=2,col="pink")
sr=sr=function(age){730.910*(1+-0.003*exp(-0.117*age^2.181))^(1/0.003)}
curve(sapply(x,sr),from=-2,to=9,add=TRUE,lwd=2,col="orange")
legend(5.75,450, c("Schnute","von Bertalanffy","Gompertz","Logistic","Schnute-Richards"),lwd=c(2.5,2.5),col=c("pink","green","blue","yellow","orange"))