#simulated corvina data. Derek Bolser

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

#load packages
library(FSA)
library(nlstools)

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

#import data

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

#sub sample 1 yr data by age
age1 <- data1[sample(1:117, 117,
                      replace=FALSE),]
age2 <- data1[sample(118:272, 155,
                      replace=FALSE),]
age3 <- data1[sample(273:354, 82,
                      replace=FALSE),]
age4 <- data1[sample(355:501, 147,
                      replace=FALSE),]
age5 <- data1[sample(502:662, 161,
                      replace=FALSE),]
age6 <- data1[sample(663:721, 59,
                      replace=FALSE),]
age7 <- data1[sample(722:744, 23,
                      replace=FALSE),]
age8 <- data1[sample(745:749, 5,
                      replace=FALSE),]

#obtain stdevs
sd1=sd(age1$TL)
sd2=sd(age2$TL)
sd3=sd(age3$TL)
sd4=sd(age4$TL)
sd5=sd(age5$TL)
sd6=sd(age6$TL)
sd7=sd(age7$TL)
sd8=sd(age8$TL)

#means
m1=mean(age1$TL)
m2=mean(age2$TL)
m3=mean(age3$TL)
m4=mean(age4$TL)
m5=mean(age5$TL)
m6=mean(age6$TL)
m7=mean(age7$TL)
m8=mean(age8$TL)

#function that returns a matrix, and takes column vectors as arguments for mean and sd
normv <- function( n , mean , sd ){
  out <- rnorm( n*length(mean) , mean = mean , sd = sd )
  return( matrix( out , ncol = n , byrow = FALSE ) )
}

#generate data. normal distribution.
set.seed(1)
sim1=normv( 83 , m1 , sd1 )
set.seed(1)
sim2=normv( 45 , m2 , sd2 )
set.seed(1)
sim3=normv( 118 , m3 , sd3 )
set.seed(1)
sim4=normv( 53 , m4 , sd4 )
set.seed(1)
sim5=normv( 39 , m5 , sd5 )
set.seed(1)
sim6=normv( 141 , m6 , sd6 )
set.seed(1)
sim7=normv( 177 , m7 , sd7 )
set.seed(1)
sim8=normv( 195 , m8 , sd8 )

#turn in to character string to paste in to excel
s1=toString(sim1)
s2=toString(sim2)
s3=toString(sim3)
s4=toString(sim4)
s5=toString(sim5)
s6=toString(sim6)
s7=toString(sim7)
s8=toString(sim8)

s1
s2
s3
s4
s5
s6
s7
s8

#import new data

#merge datasets;done in excel

simdata=simnorm_oneyr
simage=simdata$Age
simtl=simdata$TL

#plot and model data
plot(simtl~simage)

#gompertz
svG1 <- list(Linf=1000,gi=0.5,ti=1.5)
fitG1 <- nlsLM(simtl~g1(simage,Linf,gi,ti),data=simdata,start=svG1)
bootG1 <- nlsBoot(fitG1)
warnings()
cbind(Ests=coef(fitG1),confint(bootG1))
gvals=coef(fitG1)


#logistic
svL1 <- list(Linf=1000,gninf=0.5,ti=2)
fitL1 <- nlsLM(simtl~l1(simage,Linf,gninf,ti),data=simdata,start=svL1)
bootL1 <- nlsBoot(fitL1)
warnings()
cbind(Ests=coef(fitL1),confint(bootL1))
lvals=Ests=coef(fitL1)

#von bertalanffy
svvb <- list(Linf=1000,K=0.2,t0=-0.1)
fitvb <- nlsLM(simtl~vb1(simage,Linf,K,t0),data=simdata,start=svvb)
bootvb <- nlsBoot(fitvb)
warnings()
cbind(Ests=coef(fitvb),confint(bootvb))
vbvals=Ests=coef(fitvb)

#schnute
svs1 <- list(a=0.3, b=1)
fitS1 <- nlsLM(simtl~Schnute(simage,case=1,t1=1,t3=8,L1=141,L3=1013,a,b),data=simdata,start=svs1)
bootS1 <- nlsBoot(fitS1)
warnings()
cbind(Ests=coef(fitS1),confint(bootS1))
svals=Ests=coef(fitS1)

#schnute-richards (Schnute 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)
fitSRsim1 <- nlsLM(simtl~(a*(1+b*exp(-c*simage^d))^(1/f)),data=simdata,start=svsr1,control=nls.lm.control(maxiter=500),trace=TRUE)
bootSRsim1 <- nlsBoot(fitSRsim1,nls.control(maxiter=500))
warnings()
cbind(coef(fitSRsim1),confint(bootSRsim1))
srvals=coef(fitSRsim1)

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

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

#view parameter estimates
coef(fitG1)
coef(fitL1)
coef(fitvb)
coef(fitS1)
coef(fitSRsim1)

#plot
plot(simtl~simage,data=simdata,pch=19,col=rgb(0,0,0,1/4),
     xlab="Age (years)",ylab="Total Length (mm)",xlim=c(0.25,8),ylim=c(0,1100))
curve(vb1(x,unlist(vbvals)),from=-1,to=10,add=TRUE,lwd=2,col="green")
curve(g1(x,unlist(gvals)),from=-1,to=10,add=TRUE,lwd=2,col="blue")
curve(l1(x,unlist(lvals)),from=-1,to=10,add=TRUE,lwd=2,col="yellow")
curve(Schnute(x,case=1,t1=1,t3=8,L1=116.7836,L3=1091.636094,a=-0.2019429,b=2.5196610),from=-1,to=10,add=TRUE,lwd=2,col="pink")
sr=function(simage){938.797*(1+-0.00458594*exp(-0.666398*simage^0.72038))^(1/0.00188908)}
curve(sapply(x,sr),from=-2,to=10,add=TRUE,lwd=2,col="orange")
legend(5.75,450, c("Schnute", "von Bertalanffy","Schnute-Richards","Gompertz","Logistic"),lwd=c(2.5,2.5),col=c("pink","green","orange","blue","yellow"))