# R/growth_length_age.R In TropFishR: Tropical Fisheries Analysis

#### Documented in growth_length_age

#' @title Estimation of growth parameter using length-at-age data
#'
#' @description  This function estimates growth parameters from
#'    length-at-age data. It
#'    allows to perform different methods: Gulland and Holt, Ford Walford plot,
#'    Chapman's method, Bertalanffy plot, or non linear least squares method (LSM).
#'
#' @param param a list consisting of following parameters:
#' \itemize{
#'   \item \strong{age}: age measurements,
#'   \item \strong{length}: corresponding lengths in cm.
#' }
#' @param method indicating which of following methods should be applied:
#'    \code{"GullandHolt"},
#'    \code{"FordWalford"}, \code{"Chapman"}, \code{"BertalanffyPlot"},
#'    or \code{"LSM"}
#' @param Linf_est BertalanffyPlot requires an estimate for Linf to derive K and t0
#' @param Linf_init initital parameter of Linf for non-linear sqaures fitting (default 10)
#' @param K_init initital parameter of K for non-linear sqaures fitting (default 0.1)
#' @param t0_init initital parameter of t0 for non-linear sqaures fitting (default 0)
#' @param CI logical; Should confidence intervals be calculated? This option only
#'    works for the LSM method. Default is FALSE.
#' @param ci.level required confidence level (for LSM method only)
#' @param age_plot sequence with ages used for plotting (LSM method only). By default
#'      age_plot = seq(min(param$age),max(param$age),0.1)
#' @param do.sim logical. Should Monte Carlo simulation be applied? Default = FALSE
#' @param nsim the number of Monte Carlo simulations to be performed,
#'    minimum is 10000 (default).
#'
#' @examples
#' # synthetical length at age data
#' dat <- list(age = rep(1:7,each = 5),
#'    length = c(rnorm(5,25.7,0.9),rnorm(5,36,1.2),rnorm(5,42.9,1.5),rnorm(5,47.5,2),
#'    rnorm(5,50.7,0.4),rnorm(5,52.8,0.5),rnorm(5,54.2,0.7)))
#' growth_length_age(dat, method = "GullandHolt")
#'
#' # Bertalaffy plot
#' growth_length_age(dat, method = "BertalanffyPlot", Linf_est = 50)
#'
#' # non linear least squares method
#' output <- growth_length_age(param = dat, method = "LSM",
#'      Linf_init = 30, CI = TRUE, age_plot=NULL)
#' summary(output$mod) #' #' @details #' Gulland and Holt plot assumes #' infinitestimal delta t (only reasonable approximation of growth parameters if delta t #' is small). Ford Walford plot and Chapman assume constant time intervals between ages #' (delta t). The Bertalanffy plot is a robust method, however it requires an estimate of Linf. As #' long as this estimate is reasonable the resulting estimate of K is reasonable. For #' a first estimate of Linf the Powell Wetherall method \link{powell_wetherall} can #' be used. Otherwise, the largest fish or the average of the ten largest fish can be #' used for a small or large sample, respectively. All lengths have to be smaller than #' Linf as otherwise the logarithm is not defined. Oldest fish (if larger than Linf) have #' to be omitted. Non-linear least squares fitting is the preferred method to estimate #' growth parameters according to Sparre and Venema (1998). If \code{CI = TRUE} the #' confidence interval of parameters is calculated and plotted. For plotting the #' confidence interval the \code{\link{predictNLS}} from the \link{propagate} package #' is applied. #' #' @return A list with the input parameters and following parameters: #' \itemize{ #' \item \strong{x}: independent variable used for regression analysis, #' \item \strong{y}: dependent variable used for regression analysis, #' \item \strong{mod}: (non) linear model, #' \item \strong{Linf}: infinite length for investigated species in cm [cm], #' \item \strong{K}: growth coefficent for investigated species per year [1/year], #' \item \strong{t0}: theoretical time zero, at which individuals of this species hatch #' (only for Bertalanffy plot and LSM method). #' \item \strong{estimates}: dataframe with growth parameters and confidence intervals #' (only if LSM method was applied). #' } #' #' @importFrom propagate predictNLS #' @importFrom graphics abline lines plot segments #' @importFrom stats lm nls predict aggregate confint #' @importFrom grDevices adjustcolor #' @import MASS #' #' @references #' Sparre, P., Venema, S.C., 1998. Introduction to tropical fish stock assessment. #' Part 1. Manual. \emph{FAO Fisheries Technical Paper}, (306.1, Rev. 2). 407 p. #' #' @export growth_length_age <- function(param, method, Linf_est = NA, Linf_init = 10, K_init = 0.1, t0_init = 0, CI = FALSE, ci.level = 0.95, age_plot = NULL, do.sim = FALSE, nsim = 10000){ res <- param t <- res$age
Lt <- res$length switch(method, "GullandHolt"={ # Can not deal with multiple observations (lengths) per age if(length(unique(t)) == length(t)){ delta_t <- diff(t) }else{ delta_t <- diff(unique(t)) Lt <- aggregate(Lt,list(t = t),mean,na.rm= TRUE)$x
}

# growth rate
y <- rep(NA,length(delta_t))
for(i in 1:length(delta_t)){
y[i] <- (Lt[i+1] - Lt[i])/delta_t[i]
}

#mean Lt
x <- rep(NA,(length(Lt)-1))
for(i in 1:(length(Lt)-1)){
x[i] <- (Lt[i] + Lt[i+1])/ 2
}

mod <- lm(y ~ x)
sum_mod <- summary(lm(y ~ x))
a <- sum_mod$coefficients[1] b <- sum_mod$coefficients[2]
R2 <- sum_mod$r.squared K <- -b Linf <- -a/b plot(y ~ x, pch = 16, ylim = c(0,max(y,na.rm=TRUE)*1.1), xlab = "mean L(t)", ylab=expression(paste(delta,"L / ",delta,"t")), main = "Gulland and Holt") abline(a,b) tmp <- list(delta_t=delta_t,x=x,y=y,R2=R2) }, "FordWalford"={ # Can not deal with multiple observations (lengths) per age if(length(unique(t)) == length(t)){ delta_t <- diff(t) }else{ delta_t <- diff(unique(t)) Lt <- aggregate(Lt,list(t = t),mean,na.rm= TRUE)$x
}

if(!all(delta_t == 1)) stop("The Ford Walford method assumes constant time intervals!")

delta_t <- delta_t[1]

y <- Lt[2:length(Lt)]
x <- Lt[-length(Lt)]

mod <- lm(y ~ x)
sum_mod <- summary(lm(y ~ x))
a <- sum_mod$coefficients[1] b <- sum_mod$coefficients[2]
R2 <- sum_mod$r.squared K <- -1/delta_t * log(b) Linf <- a/(1-b) plot(y ~ x, pch = 16, ylim = c(0,max(y,na.rm=TRUE)*1.1), xlim = c(0,max(x,na.rm=TRUE)*1.1), xlab = "L(t)", ylab=expression(paste("L(t+",delta,"t)")), main = "Ford Walford plot") abline(a,b) abline(1,1, lty=2) segments(x0 = 0, y0 = Linf, x1 = Linf, y1 = Linf, lty=2) segments(x0 = Linf, y0 = 0, x1 = Linf, y1 = Linf, lty=2) tmp <- list(delta_t=delta_t,x=x,y=y,R2=R2) }, "Chapman"={ # Can not deal with multiple observations (lengths) per age if(length(unique(t)) == length(t)){ delta_t <- diff(t) }else{ delta_t <- diff(unique(t)) Lt <- aggregate(Lt,list(t = t),mean,na.rm= TRUE)$x
}

if(!all(delta_t == 1)) stop("The Chapman method assumes constant time intervals!")

delta_t <- delta_t[1]

Lt_short <- Lt[-length(Lt)]
y <- Lt[2:length(Lt)] - Lt_short
x <- Lt_short

mod <- lm(y ~ x)
sum_mod <- summary(lm(y ~ x))
a <- sum_mod$coefficients[1] b <- sum_mod$coefficients[2]
R2 <- sum_mod$r.squared c <- -b K <- -(1/delta_t) * log(1+b) Linf <- -a/b plot(y ~ x, pch = 16, ylim = c(0,max(y,na.rm=TRUE)*1.1), xlim = c(0,max(x,na.rm=TRUE)*1.1), xlab = "L(t)", ylab=expression(paste("L(t+",delta,"t)-L(t)")), main = "Chapman") abline(a,b) tmp <- list(delta_t=delta_t,x=x,y=y,R2=R2) }, "BertalanffyPlot"={ delta_t <- diff(t) #Linf should be given: otherwise using optim? if(is.na(Linf_est)) stop("For the Bertalanffy plot you have to provide an estimate for Linf_est!") y <- -log(1 - Lt / Linf_est) x <- t mod <- lm(y ~ x) sum_mod <- summary(lm(y ~ x)) a <- sum_mod$coefficients[1]
b <- sum_mod$coefficients[2] R2 <- sum_mod$r.squared
K <- b
t0 <- a/-K
Linf <- Linf_est

plot(y ~ x, pch = 16, ylim = c(0,max(y,na.rm=TRUE)*1.1),
xlim = c(0,max(x,na.rm=TRUE)*1.1),
xlab = "t",
ylab=expression(paste("-ln(1-L(t)/Linf)")),
main = "Bertalanffy plot")
abline(a,b)

tmp <- list(x=x,y=y,R2=R2,t0=t0)

if(max(Lt, na.rm=TRUE) > Linf_est) writeLines("Some lengths were larger as the Linf value which was \nprovided. These lengths are omitted.")
},
"LSM"={

nls_mod <- nls(Lt ~ (Linf * (1 - exp(-K * (t - t0)))),
start = list(Linf = Linf_init, K = K_init, t0 = t0_init))

if(!CI){
plot(Lt ~ t,
ylab = "L(t)",
xlab = "t(age)",
main = "Non-linear least squares method")
lines(t, predict(nls_mod))
}
sum_mod <- summary(nls_mod)
Linf <- sum_mod$coefficients[1] K <- sum_mod$coefficients[2]
t0 <- sum_mod$coefficients[3] mod <- nls_mod tmp <- list(t0 = t0) if(CI){ suppressMessages(cis <- confint(nls_mod, level = ci.level)) nls_res <- data.frame(Names=c("Linf","K","t0"), Value=c(round(Linf,2), round(K,2), round(t0,2)), Lower.CI = c(round(cis[1,1],2), round(cis[2,1],2), round(cis[3,1],2)), Upper.CI = c(round(cis[1,2],2), round(cis[2,2],2), round(cis[3,2],2))) tmp$nls_res = nls_res

# plot with confidence interval
if(is.null(age_plot)){
age_plot <- seq(min(floor(t)),max(ceiling(t)),0.1)
}else age_plot <- age_plot
# Taylor error propagation and Monte Carlo simulation for confidence interval
sink(tempfile())
pred_L <- suppressMessages(propagate::predictNLS(nls_mod,
do.sim = do.sim,
nsim = nsim,
newdata = data.frame(t = age_plot)))
sink()
# Taylor propagation
if(!do.sim){
predVals <- cbind(age_plot,as.data.frame(pred_L$summary[,c(1,5,6)])) } # Monte Carlo simulation if(do.sim){ predVals <- cbind(age_plot,as.data.frame(pred_L$summary[,c(7,11,12)]))
}
names(predVals) <- c("age_plot","fit","lower","upper")
plot(Lt ~ t, type = "n", ylab = "L(t)", xlab = "t(age)", xlim=c(min(age_plot),max(age_plot)),
main = "Non-linear least squares method")
polygon(c(age_plot,rev(age_plot)),
c(predVals$lower,rev(predVals$upper)),
border = FALSE, col = adjustcolor("dodgerblue", alpha.f = 0.4))
points(t, Lt)
lines(age_plot, predict(nls_mod,newdata = data.frame(t=age_plot)))
}
},

stop(paste("\n",method, "not recognised, possible options are \n",
"\"GullandHolt\", \"FordWalford\", \"Chapman\" \n",
"\"BertalanffyPlot\", and \"LSM\""))
)

res2 <- res
if(exists("delta_t", where = tmp)) res2$delta_t <- delta_t if(exists("x", where = tmp)) res2$x <- x
if(exists("y", where = tmp)) res2$y <- y if(exists("R2", where = tmp)) res2$r2 <- R2

ret <- c(res2,list(mod = mod,
Linf = Linf,
K = K))

if(exists("t0", where = tmp)) ret$t0 <- t0 if(exists("nls_res", where = tmp)) ret$estimates <- nls_res

return(ret)
}


## Try the TropFishR package in your browser

Any scripts or data that you put into this service are public.

TropFishR documentation built on June 29, 2018, 5 p.m.