Nothing
#' @title Virtual Population Analysis (VPA)
#'
#' @description This function applies the Virtual Population Analysis (VPA) or
#' Cohort analysis (CA). Methods used to estimate stock biomass and fishing
#' mortality per age/length group.
#'
#' @param param a list consisting of following parameters:
#' \itemize{
#' \item \code{midLengths} or \code{age}: midpoints of the length class
#' (length-frequency data) or ages (age composition data),
#' \item \code{Linf}: infinite length for investigated species in cm [cm],
#' \item \code{K}: growth coefficent for investigated species per year [1/year],
#' \item \code{t0}: theoretical time zero, at which individuals of this species hatch,
#' \item \code{M}: natural mortality [1/year] (numeric value or vector of identical
#' length than midLengths),
#' \item \code{a}: length-weight relationship coefficent (W = a * L^b; for kg/cm3),
#' \item \code{b}: length-weight relationship coefficent (W = a * L^b),
#' \item \code{catch}: catch as vector for pseudo cohort analysis,
#' or a matrix with catches of subsequent years to follow a real cohort.
#' For age-based VPA/CA catch has to be provided in numbers, e.g. '000 individuals
#' for length-based VPA/CA catch can also be provided in weight, e.g. kg (use
#' argument \code{catch_unit}).
#' }
#' @param catch_columns numerical; indicating the column of the catch matrix which should be
#' used for the analysis.
#' @param catch_unit optional; a character indicating if the catch is provided in weight
#' ("tons" or "kg") or in thousand individuals ("'000")
#' @param catch_corFac optional; correction factor for catch, in case provided
#' catch does spatially or temporarily not reflect catch for fishing ground of
#' a whole year.
#' @param terminalF the fishing mortality rate of the last age/length group.
#' @param terminalE the exploitation rate of the last age/length group.
#' @param analysis_type determines which type of assessment should be done,
#' options: "VPA" for age or length-based Virtual Population Analysis, "CA" for age- or length-based
#' Cohort Analysis. Default is "VPA".
#' @param algorithm an Algorithm to use to solve for fishing mortality. The default
#' setting \code{"new"} uses \code{\link[stats]{optimise}},
#' while \code{"old"} uses the algorithm described by Sparre and Venema (1998).
#' @param plus_group logical; indicating if the last length group is a plus group (default: TRUE).
#' @param plot logical; indicating whether a plot should be printed
#'
#' @details The main difference between virtual population analysis (VPA) and cohort
#' analysis (CA) is the step of calculating the fishing mortality per age class or
#' length group. While CA works with an approximation by assuming that all fish are
#' caught during a single day, which makes the calcualtion easier, VPA assumes that
#' the fish are caught continuously, which has to be solved by the trial and error
#' method (Sparre and Venema, 1998).
#' For the age-based VPA/CA the catch has to be provided in numbers (or '000 numbers),
#' while for the length-based VPA/CA the catch can also be provided in weight (tons or kg) by
#' using the argument \code{catch_unit}.
#' The catch has to be representative for fished species, that means there should not be
#' other fisheries fishing the same stock. If this is the case \code{catch_corFac} can
#' be used as a raising factor to account for the proportion of fish caught by other
#' fisheries.
#' When the model should follow a real cohort instead of a pseudo cohort, \code{catch}
#' has to be provided as matrix. The model then starts to follow the first age class
#' in the first column.
#' If \code{catch} matrix is shorter than the number of age classes, the age or length
#' classes without catch information are omitted. It is recommended to only
#' follow a real cohort if there is enough information for all age classes
#' (test with: \code{dim(catch)[1] <= dim(catch)[2]}).
#' If \code{plus_group} is TRUE a different calculation for the survivors of the last length group
#' is used (for more details please refer to Sparre & Venema (1998)).
#'
#' @keywords function VPA mortality F stock biomass cohort
#'
#' @examples
#' #_______________________________________________
#' # Virtual Popuation Analysis with age-composition data
#' data(whiting)
#' output <- VPA(param = whiting, catch_columns = 1, terminalE = 0.5, analysis_type = "VPA")
#' plot(output)
#'#_______________________________________________
#' # Pope's Cohort Analysis with age-composition data
#' data(whiting)
#' VPA(whiting, terminalE = 0.5, catch_columns = 3, analysis_type = "CA",
#' plot= TRUE, plus_group = TRUE)
#'
#'#_______________________________________________
#' # Virtual population analysis with length-composition data
#' data(hake)
#' VPA(hake, terminalE = 0.5, analysis_type = "VPA", plot = TRUE,
#' catch_unit = "'000", plus_group = TRUE)
#'#_______________________________________________
#' # Jones's Cohort Analysis with length-composition data
#' data(hake)
#' VPA(hake, terminalE = 0.5, analysis_type = "CA", plot = TRUE,
#' catch_unit = "'000", plus_group = TRUE)
#'
#' @return A list with the input parameters and following list objects:
#' \itemize{
#' \item \strong{classes.num}: numeric age classes or length groups (without plus sign),
#' \item \strong{catch.cohort}: a vector with the catch values which were used for
#' the analysis (exists only if catch was a matrix),
#' \item \strong{FM_calc}: a vector with the ifshing mortality (M),
#' \item \strong{Z}: a vector with the total mortality (Z),
#' \item \strong{survivors}: a vector with the number of fish surviving to the
#' next age class or length group (same unit than input catch vector),
#' \item \strong{annualMeanNr}: ta vector with the mean number of fish per year
#' (same unit than input catch vector),
#' \item \strong{meanBodyWeight}: a vector with the mean body weight in kg,
#' \item \strong{meanBiomassTon}: a vector with the mean biomass in tons,
#' \item \strong{YieldTon}: a vector with the yield in tons,
#' \item \strong{natLoss}: a vector with the number of fish died due
#' to natural mortality,
#' \item \strong{plot_mat}: matrix with rearranged survivors, nat losses
#' and catches for plotting;
#' }
#'
#' @importFrom graphics plot
#' @importFrom stats optimise
#'
#' @references
#' Jones, R., 1984. Assessing the effects of changes in exploitation pattern using length
#' composition data (with notes on VPA and cohort analysis). \emph{FAO Fish.Tech.Pap.},
#' (256): 118p.
#'
#' Jones, R., 1990. Length-cohort analysis: the importance of choosing the correct growth
#' parameters. \emph{Journal du Conseil: ICES Journal of Marine Science}, 46(2), 133-139
#'
#' Pope, J.G., 1972. An investigation of the accuracy of virtual population analysis using
#' cohort analysis. \emph{Res.Bull.ICNAF}, (9):65-74
#'
#' Pope, J.G., 1979. A modified cohort analysis in which constant natural mortality is
#' replaced by estimates of predation levels. \emph{ICES C.M.} 1979/H:16:7p. (mimeo)
#'
#' 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.
#'
#' References for weight-length relationship parameters (a & b):
#' Dorel, D., 1986. Poissons del'Atlantique nord-est relations taille-poids.
#' Institut Francais de Recherche pour l'Exploitation de la Mer. Nantes, France. 165 p.
#'
#' @export
VPA <- function(param,
catch_columns = NA, catch_unit = NA, catch_corFac = NA,
terminalF = NA, terminalE = NA,
analysis_type = "VPA", algorithm = "new",
plus_group = TRUE, plot = FALSE){
res <- param
if(is.na(catch_columns[1])) catch <- res$catch
if(!is.na(catch_columns[1])){
catchmat <- res$catch[,(catch_columns)]
if(length(catch_columns) > 1){
catch <- rowSums(catchmat, na.rm = TRUE)
}else catch <- catchmat
}
## PARAMETERS
if("par" %in% names(res)){
M <- res$par$M
Linf <- res$par$Linf
K <- res$par$K
t0 <- ifelse("t0" %in% names(res),res$par$t0,0)
C <- ifelse("C" %in% names(res$par), res$par$C, 0)
ts <- ifelse("ts" %in% names(res$par), res$par$ts, 0)
a <- res$par$a
b <- res$par$b
}else{
M <- res$M
Linf <- res$Linf
K <- res$K
t0 <- ifelse("t0" %in% names(res),res$t0,0)
C <- ifelse("C" %in% names(res), res$C, 0)
ts <- ifelse("ts" %in% names(res), res$ts, 0)
a <- res$a
b <- res$b
}
## maybe growth parameters in par but not others...
if(is.null(a)) a <- res$a
if(is.null(b)) b <- res$b
if(is.null(M)) M <- res$M
## CHECKS
if(is.null(a) || is.null(b)) stop("VPA requires information about the length-weight relationship. Please provide 'a' and 'b' estimates in res.")
if(is.null(M)) stop("Please provide a natural mortality estimate 'M' in res.")
#HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH#
# AGE BASED VPA AND COHORT ANALYSIS #
#HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH#
if("age" %in% names(res) == TRUE){
classes <- as.character(res$age)
# Error message if catch and age do not have same length
if(inherits(catch,'matrix') || inherits(catch,'data.frame')){
#if(length(classes) != length(catch[,1])) stop("Age/length classes and catch do not have the same length!")
if(length(classes) != length(diag(as.matrix(catch)))) warning("Age/length classes and the real cohort in the catch matrix do not have the same length. The missing age/length classes will be omitted.")
}else if(inherits(catch,'numeric')){
if(length(classes) != length(catch)) stop("Age/length classes and catch do not have the same length!")
}
if(length(M) == length(classes)){
M_vec <- M
}else if(length(M) > 1){
writeLines(noquote("The number of natural mortality values does not correspond to the number of length classes. \nOnly the first value will be used."))
M_vec <- rep(M[1],length(classes))
}else{
M_vec <- rep(M, length(classes))
}
if(!is.na(terminalF)){
terminalE <- terminalF / (terminalF + M_vec[length(M_vec)])
}else if(!is.na(terminalE)){
terminalF <- M_vec[length(M_vec)] * terminalE / (1 - terminalE)
}else{
stop("Please provide either the terminal exploitation rate (terminalE) or the terminal fishing mortality rate (terminalF)!")
}
terminalZ <- terminalF + M_vec[length(M_vec)]
# create column without plus group (sign) if present
classes.num <- do.call(rbind,strsplit(classes, split="\\+"))
classes.num <- as.numeric(classes.num[,1])
if(inherits(catch,"matrix") || inherits(catch,"data.frame")){
writeLines(noquote("A catch matrix was provided: The VPA/CA will follow the 'real' cohort, assuming yearly intervals in the catch matrix. If you want to perform the VPA/CA with a pseudo cohort please use the argument 'catch_columns' specifying which column(s) of the catch matrix to use."))
#find cohort to analyse
real.cohort <- diag(as.matrix(catch))
catch.cohort <- c(real.cohort,
rep(NA,length(classes.num) - length(real.cohort)))
if(length(classes.num) != length(real.cohort)){
catch.cohort <- real.cohort
classes.num <- classes.num[1:length(catch.cohort)]
}
if(length(M_vec) != length(classes.num)){
M_vec <- M_vec[1:length(classes.num)]
}
}
if(inherits(catch,'numeric')){
catch.cohort <- catch
}
#Correct catch if not representative for one year
if(!is.na(catch_corFac)) catch_cor <- catch.cohort * catch_corFac
if(is.na(catch_corFac)) catch_cor <- catch.cohort
# translate catch into individuals
if(catch_unit %in% c("tons", "t", "T", "Tons", "tonnes", "Tonnes")){
stop("The age-based VPA/CA is not applicable to catch in weight. Please set the argument 'catch_unit' either to NA, '000, or '000000 for catch in numbers.")
}else if(catch_unit %in% c("kg", "Kg", "KG", "kilo", "KILO", "kilogramm", "Kilogramm")){
stop("The age-based VPA/CA is not applicable to catch in weight. Please set the argument 'catch_unit' either to NA, '000, or '000000 for catch in numbers.")
}else if(catch_unit %in% c("'000","1000","1e3")){
catch_numbers <- catch_cor * 1000
}else if(catch_unit %in% c("'000000","1000000","1e6","'000.000")){
catch_numbers <- catch_cor * 1000000
}else if(!is.na(catch_unit)){
stop(paste0(catch_unit, " not known. Please use either NA, '000, or '000000 for catch in numbers."))
}else{
warning("You did not specify catch_unit. The Method assumes that catch is provided in numbers!")
catch_numbers <- catch_cor
}
#Survivors #N(L1)=(N(L2)*H(L1,L2)+C(L1,L2)) *H(L1,L2)
survivors <- rep(NA,length(classes.num))
# survivors last size class
lastLengthClass <- max(which(!is.na(catch_numbers)),na.rm=TRUE) ###
if(!plus_group) survivors[lastLengthClass] <-
catch_numbers[lastLengthClass] / (terminalE * (1 - exp(-terminalZ)))
if(plus_group) survivors[lastLengthClass] <- catch_numbers[lastLengthClass] / terminalE
# Age-based Cohort Analysis (Pope's cohort analysis)
if(analysis_type == "CA"){
# other survivors
for(x3 in (lastLengthClass-1):1){
survivors[x3] <- (survivors[x3+1] * exp((M_vec[x3]/2)) +
catch_numbers[x3] ) * exp((M_vec[x3]/2))
}
#F
FM_calc <- rep(NA,length(classes.num))
FM_calc[lastLengthClass] <- terminalF
for(x5 in 1:(lastLengthClass-1)){
FM_calc[x5] <- log(survivors[x5]/survivors[x5+1]) - M_vec[x5]
}
}
# Traditional VPA
if(analysis_type == "VPA"){
#other survivors and fishing mortality
FM_calc <- rep(NA,length(classes.num))
FM_calc[lastLengthClass] <- terminalF
for(num_class in (lastLengthClass-1):1){
sur.C <- catch_numbers[num_class]
sur.Ntplus1 <- survivors[(num_class+1)]
sur.M <- M_vec[num_class]
LHS <- sur.C / sur.Ntplus1
sur.F <- 0
if(algorithm == "old"){
LHS <- sur.C / sur.Ntplus1
sur.F <- 0
seqi <- c(1e-1,1e-2,1e-3,1e-4,1e-5,1e-6,1e-7)
#trail and error
for(y in seqi){
stepi <- y
for(x in seq(sur.F,10,stepi)){
sur.F <- x
RHS <- (sur.F/(sur.F + sur.M)) * (exp(sur.F+sur.M) - 1)
if(LHS-RHS < 0) break
}
sur.F = x-stepi
}
}
if(algorithm == "new"){
Fcalc <- function(sur.F=sur.M){
((sur.F/(sur.F+sur.M)) * (exp(sur.F+sur.M) - 1) - (sur.C / sur.Ntplus1))^2
}
tmp <- optimise(Fcalc, interval=c(0,100))
sur.F <- tmp$min
}
#fill F
FM_calc[num_class] <- sur.F
#fill survivors
survivors[num_class] <- survivors[(num_class+1)] *
exp(sur.F + sur.M)
}
}
# Z
Z <- M_vec + FM_calc
#Annual mean Nr
deads <- abs(diff(survivors))
annualMeanNr <- deads / Z[-length(survivors)]
# annualMeanNr[length(survivors)] <- NA
annualMeanNr[length(survivors)] <- survivors[length(survivors)] / Z[length(survivors)]
#FOR PLOT
#Survivors rearranged
survivors_rea <- rep(NA,length(classes.num))
for(x8 in 1:(length(survivors_rea)-1)){
survivors_rea[x8] <- survivors[x8+1]
}
survivors_rea[length(survivors_rea)] <- 0
#Calculate natural losses
natLoss <- survivors - survivors_rea - catch_numbers
#put together in dataframe
df.VPAnew <- data.frame(survivors = survivors_rea,
nat.losses = natLoss,
catch = catch_numbers,
FM_calc = FM_calc)
#transpose matrix for barplot function
df.VPAnew <- t(as.matrix(df.VPAnew))
colnames(df.VPAnew) <- classes.num
#save all in list
ret <- c(res,list(
classes.num = classes.num,
catch.cohort = catch.cohort,
FM_calc = FM_calc,
Z = Z,
survivors_L1 = survivors,
survivors_L2 = survivors_rea,
annualMeanNr = annualMeanNr,
natLoss = natLoss,
plot_mat = df.VPAnew))
class(ret) <- "VPA"
# plot results
if(plot==TRUE) try(plot(ret))
return(ret)
}
#HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH#
# LENGTH BASED VPA AND COHORT ANALYSIS #
#HHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH#
if("midLengths" %in% names(res) == TRUE &&
(inherits(catch,"matrix") || inherits(catch,"data.frame"))){
if(is.na(catch_columns[1])) stop("The length-based Cohort analysis is not applicable to length frequency data. Please provide catch as vector or use the argument 'catch_columns'.")
if(!is.na(catch_columns[1])){
catch <- rowSums(catch)
}
}
if((inherits(catch,'numeric') || inherits(catch,'integer')) &&
"midLengths" %in% names(res) == TRUE){
classes <- as.character(res$midLengths)
# Error message if catch and age do not have same length
if(inherits(catch,"matrix") || inherits(catch,"data.frame")){
if(length(classes) != length(catch[,1])) stop("Midlengths and catch do not have the same length!")
}else if(inherits(catch,'numeric')){
if(length(classes) != length(catch)) stop("Midlengths and catch do not have the same length!")
}
if(length(M) == length(classes)){
M_vec <- M
}else if(length(M) > 1){
writeLines(noquote("The number of natural mortality values does not correspond to the number of length classes. \nOnly the first value will be used."))
M_vec <- rep(M[1],length(classes))
}else{
M_vec <- rep(M, length(classes))
}
if(!is.na(terminalF)){
terminalE <- terminalF / (terminalF + M_vec[length(M_vec)])
}else if(!is.na(terminalE)){
terminalF <- M_vec[length(M_vec)] * terminalE / (1 - terminalE)
}else{
stop("Please provide either the terminal exploitation rate (terminalE) or the terminal fishing mortality rate (terminalF)!")
}
terminalZ <- terminalF + M_vec[length(M_vec)]
# correct catch with raising factor
if(!is.na(catch_corFac)) catch_cor <- catch * catch_corFac
if(is.na(catch_corFac)) catch_cor <- catch
# create column without plus group (sign) if present
classes.num <- do.call(rbind, strsplit(classes, split="\\+"))
classes.num <- as.numeric(classes.num[,1])
#calculate size class interval
interval <- classes.num[2] - classes.num[1]
# lower and upper length vectors
lowerLength <- classes.num - (interval / 2)
upperLength <- classes.num + (interval / 2)
if(plus_group) upperLength[length(upperLength)] <- Linf
#Mean body weight
# FAO manual:
meanBodyWeight <- a * ((lowerLength + upperLength)/2)^b # a * classes.num ^ b
# same as what provided in FAO manual: a * ((lowerLength + upperLength)/2)^b
#meanBodyWeight <- meanBodyWeight / 1000 # in kg
#according to Beyer (1987) (FISAT II)
# meanBodyWeight <- (1/(upperLength - lowerLength)) * (a / (b + 1)) * (upperLength^(b+1) - lowerLength^(b+1))
# translate catch in tons into numbers
if(catch_unit %in% c("tons", "t", "T", "Tons", "tonnes", "Tonnes")){
catch_numbers <- (catch_cor * 1000) / meanBodyWeight
}else if(catch_unit %in% c("kg", "Kg", "KG", "kilo", "KILO", "kilogramm", "Kilogramm")){
catch_numbers <- catch_cor / meanBodyWeight
}else if(catch_unit %in% c("'000","1000","1e3")){
catch_numbers <- catch_cor * 1000
}else if(catch_unit %in% c("'000000","1000000","1e6","'000.000")){
catch_numbers <- catch_cor * 1000000
}else if(!is.na(catch_unit)){
stop(paste0(catch_unit, " not known. Please use either 'tons' or 'kg' for catch in weight or NA, '000, or '000000 for catch in numbers."))
}else{
warning("You did not specify catch_unit. The Method assumes that catch is provided in numbers!")
catch_numbers <- catch_cor
}
# t of lower length classes
t_L1 <- t0 - (1/K) * log(1 - (lowerLength / Linf))
# t of lower upper classes
t_L2 <- t0 - (1/K) * log(1 - (upperLength / Linf))
if(upperLength[length(upperLength)] > Linf){
writeLines(noquote("Upper limit of last length class is larger than Linf, \nconsider creating lower plus group or set the argument plus_group = TRUE."))
}
# delta t
# dt <- t_L2 - t_L1
dt <- (1/K) * log((Linf - lowerLength)/(Linf - upperLength))
#Survivors
survivors <- rep(NA, length(classes.num))
# survivors last size class
lastLengthClass <- max(which(!is.na(catch_numbers)),na.rm=TRUE)
if(plus_group) survivors[length(survivors)] <- catch_numbers[length(survivors)] / terminalE
if(!plus_group){
survivors[length(survivors)] <- catch_numbers[length(survivors)] /
(terminalE * (1 - exp(-terminalZ * dt[length(survivors)])))
}
### Jones' Length-based Cohort Analysis
if(analysis_type == "CA"){
# H (L1,L2) #H(L1,L2)=((Linf-L1)/Linf-L2)^(M/2K)
H <- ((Linf - lowerLength)/(Linf - upperLength))^(M_vec/(2*K))
# other survivors
for(x3 in (length(survivors)-1):1){
survivors[x3] <- (survivors[x3+1] * H[x3] + catch_numbers[x3]) * H[x3]
}
# F/Z #F(L1,L2)/Z(L1,L2)=C(L1,L2)/(N(L1)-N(L2))
deads <- abs(diff(survivors))
F_Z <- catch_numbers[-length(survivors)] / deads
F_Z[length(survivors)] <- terminalE
#F # F = M * (F_Z / 1-F_Z)
FM_calc <- M_vec * F_Z / (1 - F_Z)
}
### Length-based VPA
if(analysis_type == "VPA"){
#other survivors and fishing mortality
FM_calc <- rep(NA,length(classes.num))
FM_calc[lastLengthClass] <- terminalF
for(num_class in (lastLengthClass-1):1){
sur.C <- catch_numbers[num_class]
sur.Ntplus1 <- survivors[(num_class+1)]
sur.M <- M_vec[num_class]
sur.dt <- dt[num_class]
LHS <- sur.C / sur.Ntplus1
sur.F <- 0
if(algorithm == "old"){
LHS <- sur.C / sur.Ntplus1
sur.F <- 0
seqi <- c(1e-1,1e-2,1e-3,1e-4,1e-5,1e-6,1e-7)
#trail and error
for(y in seqi){
stepi <- y
for(x in seq(sur.F,10,stepi)){
sur.F <- x
RHS <- (sur.F/(sur.F + sur.M)) * (exp(sur.F+sur.M) * sur.dt - 1)
if(LHS-RHS < 0) break
}
sur.F = x-stepi
}
}
if(algorithm == "new"){
Fcalc <- function(sur.F=sur.M){
((sur.F/(sur.F+sur.M)) * (exp((sur.F+sur.M) * sur.dt) - 1) - (sur.C / sur.Ntplus1))^2
}
tmp <- optimise(f = Fcalc, interval=c(0,100))
sur.F <- tmp$min
}
#fill F
FM_calc[num_class] <- sur.F
#fill survivors
survivors[num_class] <- survivors[(num_class+1)] *
exp((sur.F + sur.M) * sur.dt)
}
}
# Z
Z <- M_vec + FM_calc
#Annual mean Nr
deads <- abs(diff(survivors))
annualMeanNr <- deads / Z[-length(survivors)]
# annualMeanNr[length(survivors)] <- NA
annualMeanNr[length(survivors)] <- survivors[length(survivors)] / Z[length(survivors)]
#Mean biomass
meanBiomass <- annualMeanNr * meanBodyWeight
meanBiomassTon <- meanBiomass / 1000
#Yield
yield <- catch_numbers * meanBodyWeight
yieldTon <- yield / 1000
#FOR PLOT
#Survivors rearranged
survivors_rea <- rep(NA,length(classes.num))
for(x8 in 1:(length(survivors_rea)-1)){
survivors_rea[x8] <- survivors[x8+1]
}
survivors_rea[length(survivors_rea)] <- 0
#Calculate natural losses
natLoss <- survivors - survivors_rea - catch_numbers
#put together in dataframe
df.VPAnew <- data.frame(survivors = survivors_rea,
nat.losses = natLoss,
catch = catch_numbers,
FM_calc = FM_calc,
meanBodyWeight = meanBodyWeight,
meanBiomassTon = meanBiomassTon)
#transpose matrix for barplot function
df.VPAnew <- t(as.matrix(df.VPAnew))
colnames(df.VPAnew) <- classes.num
#save all in list
ret <- c(res,list(
classes.num = classes.num,
FM_calc = FM_calc,
Z = Z,
meanBodyWeight = meanBodyWeight,
survivors_L1 = survivors,
survivors_L2 = survivors_rea,
catch_numbers = catch_numbers,
annualMeanNr = annualMeanNr,
meanBiomassTon = meanBiomassTon,
yieldTon = yieldTon,
natLoss = natLoss,
plot_mat = df.VPAnew))
class(ret) <- "VPA"
# plot results
if(plot == TRUE & all(!is.na(survivors)) & all(!is.na(natLoss))) try(plot(ret))
return(ret)
}# stop("Please choose analysis_type = 'CA' for length composition data!")
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.