merging old and new VPA-CA methodlogy

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))

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)