Add Pauly et al. (1992) growth cessation model.

[droglenc]

Requested by Erica Newton.

[droglenc]

This code from Melanie Beguer-Pon may be useful as a start ...

#------------------------------------------------------------------------------------------
#Model developed by Pauly et al (1992)
#---------------------------------------------------------------------------------------------
#Set the parameters (with numbers that you think could fit but it doesn't changed
#  anything, you could put 0 for each)
Kprime=0.5  # year-1, growth parameter
Linf=63.26 # the asymptotic length, Linf=Lmax/0.95
Lo=9.3 # Size at age=0
Ts=3 #(year) date of starting higher growth
NGT=1 # (year), the no-growth time before the first period growth
p4 <- c(Kprime,Linf,Lo,Ts,NGT)

#creation of the growth function
VBGF4 <- function(x,p4) {
  if ( length(p4)==5) {
    res <- p4[2]*(1-exp(-(p4[1]*(x-(log(1-p4[3]/p4[2]))/p4[1])+(p4[1]/((2*pi)/(1-p4[5])))*(sin(((2*pi)/(1-p4[5]))*(x-p4[4]))-sin(((2*pi)/(1-p4[5]))*((log(1-p4[3]/p4[2])/p4[1])-p4[4])))))) }
  else {
    res <- "NA"
    print("dim of parameter p is incorrect") 
  }
  res
}

#unconstrained model
ResultsWithoutConstaints <- nls(formula =Size~Linf*(1-exp(-(Kprime*(Time-((log(1-Lo/Linf))/Kprime))+(Kprime/((2*pi)/(1-NGT)))*(sin(((2*pi)/(1-NGT))*(Time-Ts))-sin(((2*pi)/(1-NGT))*(((log(1-Lo/Linf))/Kprime)-Ts)))))),start = list(Lo=4,Kprime = 0.2,Ts=1,NGT=0.0,Linf=63.26),algorithm = "port",trace=TRUE) 
ResultsWithoutConstaints

coef(ResultsWithoutConstaints) 
AIC(ResultsWithoutConstaints)

#Plot the growth curve
coef(ResultsWithoutConstaints)
Lo=coef(ResultsWithoutConstaints)[1]
Kprime=coef(ResultsWithoutConstaints)[2]
Ts=coef(ResultsWithoutConstaints)[3]
NGT=coef(ResultsWithoutConstaints)[4]
Linf=coef(ResultsWithoutConstaints)[5]

p4 <- c(Kprime,Linf,Lo,Ts,NGT)

Time2=seq(2,8.2,0.2)
points(Time2,VBGF4(Time2,p4),col="red",type="l",lwd=2,lty=1) #well-asdjusted no?

jpeg(filename = "Growth_L1_PaulyEquation.jpg",width = 1000, height = 727, pointsize=70, bg = "white", res = 1200,quality=100, restoreConsole = TRUE)
dev.off()

GrowthPerformance=2*log10(Linf)+log10(Kprime)
GrowthPerformance

#---------------------------------------------------------------------
# You can try with some constraints if you want to set some parameters 
# you just have to fill the part "upper=..." and even set a lower=...
#To constrain some parameters, i.e. Linf (max size)
ResultsWithConstraint <- nls(formula =Size~Linf*(1-exp(-(Kprime*(Time-((log(1-Lo/Linf))/Kprime))+(Kprime/((2*pi)/(1-NGT)))*(sin(((2*pi)/(1-NGT))*(Time-Ts))-sin(((2*pi)/(1-NGT))*(((log(1-Lo/Linf))/Kprime)-Ts)))))),start = list(Lo=4,Kprime = 0.2,Ts=1,NGT=0.0,Linf=63.26), algorithm = "port",upper=list(Lo=4,Kprime = 0.2,Ts=1,NGT=0.0,Linf=63.26),trace=TRUE) 
ResultsWithConstraint

coef(ResultsWithConstraint)
Lo=coef(ResultsWithConstraint)[1]
Kprime=coef(ResultsWithConstraint)[2]
Ts=coef(ResultsWithConstraint)[3]
NGT=coef(ResultsWithConstraint)[4]
Linf=coef(ResultsWithConstraint)[5]

p4bis <- c(Kprime,Linf,Lo,Ts,NGT)
points(Time2,VBGF4(Time2,p4bis),col="blue",type="l",lwd=2,lty=1) 
#----------------------------------------------------------------------------------------------

[droglenc]

This function (code above) seems to be for a different model or parametrization than what is in Pauly et al. (1992). Note the use of L_0. Not clear what Kprime is, either.

[droglenc]

Not sure where to go with the equation in Pauly et al. (1992) because he defines tprime as "is obtained by subtracting from the real age (t) the total no-growth time occurring up to age t.” It seems circular to me … you need t-prime to fit the model and estimate NGT, but you need NGT to know what t-prime is.

[droglenc]

Added the model to vbFuns() and growthFunShow() in v0.8.8. Did not add to vbStarts() as I am not sure how to globally estimate starting values for ts and NGT.