From 657b8aa1cf873bbba84e07dc370556633bf53767 Mon Sep 17 00:00:00 2001 From: tokami Date: Thu, 30 Mar 2017 11:12:07 -0500 Subject: [PATCH] merging old and new VPA-CA methodlogy --- R/VPA.R | 21 +++++++++++++-------- R/plot.VPA.R | 30 +++++++++++++++++++++++++++--- man/VPA.Rd | 8 +++++--- man/plot.VPA.Rd | 6 ++++-- 4 files changed, 49 insertions(+), 16 deletions(-) diff --git a/R/VPA.R b/R/VPA.R index f99d65f..c810529 100644 --- a/R/VPA.R +++ b/R/VPA.R @@ -68,12 +68,13 @@ #' #_______________________________________________ #' # Virtual Popuation Analysis with age-composition data #' data(whiting) -#' output <- VPA(param = whiting, terminalE = 0.5, analysis_type = "VPA") +#' 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, terminalF = 0.5, analysis_type = "CA", plot= TRUE) +#' VPA(whiting, terminalE = 0.5, catch_columns = 3, analysis_type = "CA", +#' plot= TRUE, plus_group = TRUE) #' #'#_______________________________________________ #' # Virtual population analysis with length-composition data @@ -237,7 +238,7 @@ VPA <- function(param, # other survivors for(x3 in (lastLengthClass-1):1){ survivors[x3] <- (survivors[x3+1] * exp((M_vec[x3]/2)) + - catch_numbers[x3] ) * exp((M[x3]/2)) + catch_numbers[x3] ) * exp((M_vec[x3]/2)) } #F @@ -410,11 +411,13 @@ VPA <- function(param, lowerLength <- classes.num - (interval / 2) upperLength <- classes.num + (interval / 2) - #Mean body weight according to Beyer (1987) - meanBodyWeight <- (1/(upperLength - lowerLength)) * (a / (b + 1)) * (upperLength^(b+1) - lowerLength^(b+1)) + #Mean body weight + # FAO manual: + meanBodyWeight <- a * classes.num ^ b + # same as what provided in FAO manual: a * ((lowerLength + upperLength)/2)^b meanBodyWeight <- meanBodyWeight / 1000 # in kg - # old - # meanBodyWeight <- a * classes.num ^ b # same as what provided in FAO manual: a * ((lowerLength + upperLength)/2)^b + #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")){ @@ -556,7 +559,9 @@ VPA <- function(param, df.VPAnew <- data.frame(survivors = survivors_rea, nat.losses = natLoss, catch = catch_numbers, - FM_calc = FM_calc) + FM_calc = FM_calc, + meanBodyWeight = meanBodyWeight, + meanBiomassTon = meanBiomassTon) #transpose matrix for barplot function df.VPAnew <- t(as.matrix(df.VPAnew)) diff --git a/R/plot.VPA.R b/R/plot.VPA.R index 8320620..59fd303 100644 --- a/R/plot.VPA.R +++ b/R/plot.VPA.R @@ -4,6 +4,7 @@ #' mortality resulting from the \link{VPA} model. #' #' @param x list of the class \code{"VPA"} containing the results of the VPA model. +#' @param yaxis indicating which variable should be displayed on the y axis, either "numbers" or "biomass". #' @param display_last_class logical; should last age/length class be displayed in graph? #' @param xlabel Label of the x axis #' @param ylabel1 Label of the first y axis @@ -30,7 +31,9 @@ #' #' @export -plot.VPA <- function(x, display_last_class = TRUE, +plot.VPA <- function(x, + yaxis = "numbers", + display_last_class = TRUE, xlabel = NA, ylabel1 = "Population", ylabel2 = "Fishing mortality", @@ -41,8 +44,29 @@ plot.VPA <- function(x, display_last_class = TRUE, natLoss <- pes$plot_mat[2,] catch <- pes$plot_mat[3,] FM_calc <- pes$plot_mat[4,] + if(dim(pes$plot_mat)[1] == 5){ + meanBodyWeight <- pes$plot_mat[5,] + } + if(dim(pes$plot_mat)[1] == 6){ + meanBiomassTon <- pes$plot_mat[6,] + } classes.num <- as.numeric(colnames(pes$plot_mat)) - df.VPAnew <- pes$plot_mat[-4,] + + #put together in dataframe + if(yaxis == "numbers"){ + df.VPAnew <- data.frame(survivors = survivors, + nat.losses = natLoss, + catch = catch) + } + if(yaxis == "biomass"){ + df.VPAnew <- data.frame(survivors = c(meanBiomassTon[-1],0), + nat.losses = natLoss * meanBodyWeight, + catch = catch * meanBodyWeight) + } + + #transpose matrix for barplot function + df.VPAnew <- t(as.matrix(df.VPAnew)) + colnames(df.VPAnew) <- classes.num if(is.na(xlabel)){ if("age" %in% names(pes)){ @@ -58,7 +82,7 @@ plot.VPA <- function(x, display_last_class = TRUE, #save x axis positions - max_sur <- round(max(survivors+natLoss+catch,na.rm=TRUE),digits=0) + max_sur <- round(max(colSums(df.VPAnew),na.rm=TRUE),digits=0) dim_sur <- 10 ^ (nchar(max_sur)-1) max_FM <- ceiling(max(FM_calc,na.rm=TRUE)) max_clas <- max(classes.num,na.rm=TRUE) diff --git a/man/VPA.Rd b/man/VPA.Rd index 2d447d9..8ac4753 100644 --- a/man/VPA.Rd +++ b/man/VPA.Rd @@ -16,7 +16,8 @@ VPA(param, catch_columns = NA, catch_unit = NA, catch_corFac = NA, \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], + \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), \item \code{b}: length-weight relationship coefficent (W = a * L^b), \item \code{catch}: catch as vector for pseudo cohort analysis, @@ -106,12 +107,13 @@ The main difference between virtual population analysis (VPA) and cohort #_______________________________________________ # Virtual Popuation Analysis with age-composition data data(whiting) -output <- VPA(param = whiting, terminalE = 0.5, analysis_type = "VPA") +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, terminalF = 0.5, analysis_type = "CA", plot= TRUE) +VPA(whiting, terminalE = 0.5, catch_columns = 3, analysis_type = "CA", + plot= TRUE, plus_group = TRUE) #_______________________________________________ # Virtual population analysis with length-composition data diff --git a/man/plot.VPA.Rd b/man/plot.VPA.Rd index 94c89bd..efaa1b8 100644 --- a/man/plot.VPA.Rd +++ b/man/plot.VPA.Rd @@ -4,12 +4,14 @@ \alias{plot.VPA} \title{VPA plot} \usage{ -\method{plot}{VPA}(x, display_last_class = TRUE, xlabel = NA, - ylabel1 = "Population", ylabel2 = "Fishing mortality", ...) +\method{plot}{VPA}(x, yaxis = "numbers", display_last_class = TRUE, + xlabel = NA, ylabel1 = "Population", ylabel2 = "Fishing mortality", ...) } \arguments{ \item{x}{list of the class \code{"VPA"} containing the results of the VPA model.} +\item{yaxis}{indicating which variable should be displayed on the y axis, either "numbers" or "biomass".} + \item{display_last_class}{logical; should last age/length class be displayed in graph?} \item{xlabel}{Label of the x axis}