#' @title Production models with time series fitting
#'
#' @description This function applies the production models under non-equilibrium
#' conditions by applying time series fitting using non-linear least squares
#' minimisation.
#'
#' @param data a dataframe of parameters
#' \itemize{
#' \item \code{year} years,
#' \item \code{yield} catch in weight of fishery per year,
#' \item \code{effort} fishing effort per year,
#' \item \code{CPUE} catch per unit of effort per year (optional).
#' }
#' @param method indicating if Schaefer or Fox model should be applied. First assumes a
#' logistic relationship between growth rate and biomass, whereas second assumes it to
#' foolow the Gompertz distribution (Richards 1959). Default is the dynamic Schaefer model.
#' @param B0_init numeric; if realistic initial estimate for virgin biomass is available.
#' If NA initial estimate for virgin biomass is set to two times average yield of all
#' or part of yield values (see \code{B0_est}).
#' @param B0_est intital value of virgin biomass estimating using all yield values (NA) or
#' first years of time series, then provide numerical representing number of years
#' @param effort_unit multiplication factor for the unit of effort. Default is 1.
#' @param plot logical; if TRUE (default) a graph is displayed
#'
#' @keywords function biomass MSY production
#'
#' @examples
#' data(emperor)
#' prod_mod_ts(emperor, method = "Schaefer")
#' prod_mod_ts(emperor, method = "Fox")
#'
#' @return A list with the input parameters and following list objects:
#' \itemize{
#' \item \strong{Bvec}: biomass vector,
#' \item \strong{CPUE_hat}: predicted CPUE,
#' \item \strong{K}: carrying capacity,
#' \item \strong{r}: population growth rate,
#' \item \strong{q}: catchability coefficient,
#' \item \strong{MSY}: maximum sustainabale yield (MSY),
#' \item \strong{Bmsy}: biomass at MSY,
#' \item \strong{Emsy}: fishing effort at MSY
#' \item \strong{Fmsy}: fishing mortality at MSY,
#' }
#'
#' @details Either catch per unit of effort (CPUE) is inserted
#' into the model directly (by a column \code{CPUE}) or CPUE is calculated from
#' the catch and effort, then these two vectors should have required units.
#' Whenever a good estimate for the virigin biomass is available, this estimate
#' should be inserted for \code{B_init}. The default approach for the initial
#' estimate of the virgin biomass is to multiply the average yield by 2 (Dharmendra
#' and Solmundsson, 2005). Alternatively, just a part of the time series of
#' yield values can be choosen to represent the virgin biomass.
#' The minimisation procedure is based on least error sum of squares (SSE). For
#' the logistic (Schaefer) method the standard calculation of SSE is applied
#' (\code{sum((CPUE - predicted CPUE)^2)}), for
#' the method with Gompertz distribution (Fox) SSE is calculated according to
#' the Thiel's U statistic \code{sqrt(sum(CPUE - predicted CPUE)/sum(CPUE(t) - CPUE(t-1)))}
#' (Wittink, 1988).
#'
#' @importFrom graphics plot
#' @importFrom stats nlm
#'
#' @references
#' Dharmendra, D., Solmundsson, J., 2005. Stock assessment of the offshore Mauritian banks
#' using dynamic biomass models and analysis of length frequency of the Sky Emperor
#' (\emph{Lethrinus mahsena}). Fisheries Training Program The United Nations University, 61
#'
#' Hilborn, R. and Walters, C., 1992. Quantitative Fisheries Stock Assessment: Choice,
#' Dynamics and Uncertainty. Chapman and Hall, New York
#'
#' Prager, M. H., 1994. A suite of extensions to a non-equilibrium surplus production
#' model. \emph{Fishery Bulletin} 92: 374-389
#'
#' Richards, F. J., 1959. A flexible growth function for empirical use.
#' \emph{Journal of experimental Botany}, 10(2), 290-301.
#'
#' Wittink, D. R., 1988. The application of regression analysis. Allyn and Bacon. Inc.
#' Boston. MA. 324p.
#'
#'
#' @export
prod_mod_ts <- function(data, method = "Schaefer",
B0_init = NA, B0_est = NA,
effort_unit = 1, plot = TRUE){
res <- data
year <- res$year
Y <- res$Y
if("f" %in% names(res)){
f <- res$f
}else f <- Y/res$CPUE
if("CPUE" %in% names(res)){
CPUE <- res$CPUE
}else CPUE <- Y/f
# Initialisation of parameters (K, B0, r with fixed q)
if(is.na(B0_init)){
if(is.na(B0_est)){
B0 <- 2 * mean(Y)
}else B0 <- 2 * mean(Y[1:B0_est])
}else B0 <- B0_init
K <- B0 * 1.3
r <- 1
q <- mean(CPUE)/B0
input <- c(K, B0, r)
B <- B0
# Evaluating model fit
ssefn<-function(input){
K <- input[1]
B0 <- input[2]
r <- input[3]
B <- B0
Yvec <- NULL
Bvec <- NULL
CPUE_hat <- NULL # predicted CPUE
yrs <- 1:length(Y)
for (y in yrs){
if(method == "Schaefer") SY <- r * B * (1 - B/K)
if(method == "Fox") SY <- r * B * log(K/B)
Bvec <- c(Bvec, B)
CPUE_hat <- c(CPUE_hat, q * B)
B <- B + SY - Y[y]
B <- ifelse(B < 0, 0, B)
}
switch(method,
"Schaefer"={
SSE <- sum((CPUE - CPUE_hat)^2, na.rm=TRUE)
},
"Fox"={
CPUE_diff <- rep(NA,length(yrs))
for(i in 2:length(yrs)){
CPUE_diff[i] <- CPUE[i] - CPUE[i-1]
}
SSE <- sqrt(sum((CPUE - CPUE_hat)^2, na.rm=TRUE) /
(sum((CPUE_diff)^2, na.rm=TRUE)))
}
)
return(SSE)
}
estA <- nlm(ssefn, input, typsize = input, iterlim = 1000)
estA <- nlm(ssefn, estA$est, typsize = input, iterlim = 1000)
# Estimation of r and q
K <- estA$est[1]
B0 <- estA$est[2]
r <- estA$est[3]
q <- mean(CPUE) / B0
input2 <- c(r, q)
B <- B0
ssefn2 <- function(input){
r <- input[1]
q <- input[2]
B <- B0
Yvec <- NULL
Bvec <- NULL
CPUE_hat <- NULL
yrs <- 1:length(Y)
for (y in yrs){
if(method == "Schaefer") SY <- r * B * (1 - B/K)
if(method == "Fox") SY <- r * B * log(K/B)
Bvec <- c(Bvec, B)
CPUE_hat <- c(CPUE_hat, q * B)
B <- B + SY - Y[y]
B <- ifelse(B<0,0,B)
}
switch(method,
"Schaefer"={
SSE <- sum((CPUE - CPUE_hat)^2, na.rm=TRUE)
},
"Fox"={
CPUE_diff <- rep(NA,length(yrs))
for(i in 2:length(yrs)){
CPUE_diff[i] <- CPUE[i] - CPUE[i-1]
}
SSE <- sqrt(sum((CPUE - CPUE_hat)^2, na.rm=TRUE) /
(sum((CPUE_diff)^2, na.rm=TRUE)))
}
)
return(SSE)
}
estB <- nlm(ssefn2, input2, typsize=input2, iterlim=1000)
estB <- nlm(ssefn2, estB$est, typsize=estB$est, iterlim=1000)
# calculate predicted biomass for all age classes
K <- estA$est[1]
B0 <- estA$est[2]
r <- estB$est[1]
q <- estB$est[2]
B <- B0
Yvec <- NULL
Bvec <- NULL
CPUE_hat <- NULL
yrs <- 1:length(Y)
for (y in yrs){
if(method == "Schaefer") SY <- r * B * (1 - B/K)
if(method == "Fox") SY <- r * B * log(K/B)
Bvec <- c(Bvec, B)
CPUE_hat <- c(CPUE_hat, q * B)
B <- B + SY - Y[y]
}
# calculate MSY, Bmsy, Fmsy and Emsy
if(method == "Schaefer"){
MSY <- r * K / 4
Bmsy <- K / 2
Fmsy <- r / 2
Emsy <- r / (2 * q) * effort_unit
}
if(method == "Fox"){
MSY <- (K * r) / exp(1)
Bmsy <- K / exp(1)
Fmsy <- r / exp(1)
Emsy <- r / (exp(1) * q) * effort_unit
}
# create ouput list
ret <- list(
year = year,
Y = Y,
f = f,
method = method,
CPUE = CPUE,
Bvec = Bvec,
CPUE_hat = CPUE_hat,
K = K,
r = r,
q = q,
MSY = MSY,
Bmsy = Bmsy,
Emsy = Emsy,
Fmsy = Fmsy)
class(ret) <- "prod_mod_ts"
# create plot
if(plot) plot(ret)
return(ret)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.