R/PMobjects.R

In Blue-Matter/OMtool: Management Strategy Evaluation Toolkit

# All performance metric functions share a single help file: PerformanceMetric.Rd

#' Performance Metrics Methods
#'
#' Performance metric (PMs) methods for your management strategy evaluation.
#'
#' @name PerformanceMetric
#' @param MSEobj An object of class MSE
#' @param Ref Reference point for calculating the performance metric. See details.
#' @param Yrs Numeric vector of length 2 with year indices to summarize performance.
#' If NULL, the performance is summarized over all projection years.
#'
#' @details Performance Metric definitions:
#'
#' \tabular{ll}{
#' \code{P10} \tab Probability B > 0.1 BMSY \cr
#' \code{P50} \tab Probability B > 0.5 BMSY \cr
#' \code{P100} \tab Probability B > BMSY \cr
#' \code{PNOF} \tab Probability F < FMSY \cr
#' \code{LTY} \tab Probability Long-Term Yield > 0.5 Relative Yield \cr
#' \code{STY} \tab Probability Short-Term Yield > 0.5 Relative Yield \cr
#' \code{AAVY} \tab Probability AAVY < 0.2 (Average Annual Variability in Yield) \cr
#' \code{AAVE} \tab Probability AAVE < 0.2 (Average Annual Variability in Effort) \cr
#' \code{Yield} \tab Average Yield (relative to Reference Yield) \cr
#' }
#'
#' Argument `Ref` provides the ratio relative to the reference point for calculating
#' the performance metric. For biomass-based PMs (`P10`, `P50`, `P100`), this is the fraction of
#' BMSY. For `PNOF`, the fraction of FMSY. For `Yield` (and `LTY`/`STY`), the fraction of the
#' Reference Yield. For `AAVY` is it the maximum acceptable variability in yield (i.e, default
#' for `AAVY` is `Ref=0.2`)
#'
#' The `Yrs` argument defines the number of years to calculate the performance statistic over.
#' A value of `NULL`, the default for `AAVY`, `AAVE`, `P10`, `P50`, `P100`, and `PNOF`, means that the
#' performance metric is calculated over all projection years. A numeric vector of length two is used
#' to specify the first and last year, e.g, if `Yrs=c(1,10)` the performance statistic is calculated
#' over the first 10 projection years. A numeric vector of length one with positive or negative value
#' respectively can be used to specify the first *x* or last *x* years, e.g, `Yrs=10` is first 10 years,
#' and `Yrs=-10` is the last 10 years. See \code{\link{ChkYrs}} for more details.
#'
#' By default Long-Term Yield (`LTY`) is the Yield in the last ten years of the projection period in the MSE,
#' and Short-Term Yield (`STY`) is that in the first 10 years of the projection period.
#'
# #' @templateVar url performance-metrics
# #' @templateVar ref NULL
# #' @template userguide_link
#'
#' @return An object of class `PMobj`
#' @examples
#' \dontrun{
#' myMSE <- runMSE()
#' P10(myMSE)
#' P50(myMSE)
#' P100(myMSE)
#' PNOF(myMSE)
#' LTY(myMSE)
#' STY(myMSE)
#' AAVY(myMSE)
#' AAVE(myMSE)
#' Yield(myMSE)
#' }
NULL

#' Check the years to summarize performance
#'
#' @param Yrs Numeric vector of length 2 with year indices to summarize performance.
#' If NULL, the performance is summarized over all projection years.
#' `Yrs` can also be length one, in which case if it is positive it is the first `Yrs` and
#' if negative the last `Yrs` of the projection period.
#' @param MSEobj An object of class `MSE`
#'
#' @return A numeric vector of length 2 with year indices to summarize performance
#' @examples
#' \dontrun{
#' MSE <- runMSE()
#' ChkYrs(NULL, MSE) # returns c(1, MSE@proyears)
#' ChkYrs(c(2,5), MSE) # returns c(2,5)
#' ChkYrs(c(70,80), MSE) # returns c(MSE@proyears-10,MSE@proyears)
#' ChkYrs(5, MSE) # returns c(1,5)
#' ChkYrs(-5, MSE) # returns c(46,50)
#' }
#'
#' @keywords internal
#' @export
ChkYrs <- function(Yrs, MSEobj) {
   if (!methods::is(MSEobj, "MSE") & !methods::is(MSEobj, "MMSE"))
      stop("Require object of class `MSE` or class `MMSE`", call. = FALSE)

 if (is.null(Yrs)) {
   y.st <- 1
   y.end <- MSEobj@proyears
 } else {
   if (length(Yrs) == 1) {
     if (Yrs == 0) stop("Yrs must be postive or negative", call.=FALSE)
     if (Yrs < 0) {
       y.st <- MSEobj@proyears + Yrs[1] + 1
       y.end <- MSEobj@proyears
     } else {
       y.st <- 1
       y.end <- y.st + Yrs[1] - 1
     }
   } else {
     if (length(Yrs)>2) stop("Yrs must be numeric vector of length 1 or 2", call.=FALSE)

     y.st <- Yrs[1]
     y.end <- Yrs[2]

     if (Yrs[1] > Yrs[2]) stop("Yrs[1] is > Yrs[2]", call.=FALSE)
     if (any(Yrs < 1)) stop("Yrs must be positive", call.=FALSE)

     if (Yrs[2] > MSEobj@proyears) {
       message('Yrs[2] is greater than MSEobj@proyears. Setting Yrs[2] = MSEobj@proyears')
       y.end <- MSEobj@proyears
     }
     if (Yrs[1] > MSEobj@proyears) {
       message('Yrs[1] is greater than MSEobj@proyears. Setting Yrs[1] = Yrs[2] - Yrs[1]')
       y.st <- max(1,y.end - (Yrs[2] - Yrs[1]))
     }
   }
 }

 return(c(y.st, y.end))
}



#' @rdname PerformanceMetric
#' @export
P10 <- function(MSEobj=NULL, Ref=0.1, Yrs=NULL) {
 Yrs <- ChkYrs(Yrs, MSEobj)

 PMobj <- new("PMobj")
 PMobj@Name <- "Spawning Biomass relative to SBMSY"

 if (Ref !=1) {
   PMobj@Caption <- paste0('Prob. SB > ', Ref, ' SBMSY (Years ', Yrs[1], ' - ', Yrs[2], ')')
 } else {
   PMobj@Caption <- paste0('Prob. SB > SBMSY (Years ', Yrs[1], ' - ', Yrs[2], ')')
 }

 PMobj@Ref <- Ref
 PMobj@Stat <- MSEobj@SB_SBMSY[,,Yrs[1]:Yrs[2]] # Performance Metric statistic of interest - here SB/SBMSY
 PMobj@Prob <- calcProb(PMobj@Stat > PMobj@Ref, MSEobj) # calculate probability Stat > 0.1 nsim by nMP

 PMobj@Mean <- calcMean(PMobj@Prob) # calculate mean probability by MP
 PMobj@MPs <- MSEobj@MPs
 PMobj

}
class(P10) <- "PM"


#' @rdname PerformanceMetric
#' @export
P50 <- P10
formals(P50)$Ref <- 0.5
class(P50) <- "PM"

#' @rdname PerformanceMetric
#' @export
P100 <- P10
formals(P100)$Ref <- 1
class(P100) <- "PM"


#' @rdname PerformanceMetric
#' @export
PNOF <- function(MSEobj=NULL, Ref=1, Yrs=NULL) {
 Yrs <- ChkYrs(Yrs, MSEobj)
 PMobj <- new("PMobj")
 PMobj@Name <- "Probability of not overfishing (F<FMSY)"
 if (Ref !=1) {
   PMobj@Caption <- paste0('Prob. F < ', Ref, ' FMSY (Years ', Yrs[1], ' - ', Yrs[2], ')')
 } else {
   PMobj@Caption <- paste0('Prob. F < FMSY (Years ', Yrs[1], ' - ', Yrs[2], ')')
 }

 PMobj@Stat <- MSEobj@F_FMSY[,,Yrs[1]:Yrs[2]] # Performance Metric statistic of interest - here F/FMSY
 PMobj@Ref <- Ref
 PMobj@Prob <- calcProb(PMobj@Stat < PMobj@Ref, MSEobj) # calculate probability Stat < 1 nsim by nMP


 PMobj@Mean <- calcMean(PMobj@Prob) # calculate mean probability by MP
 PMobj@MPs <- MSEobj@MPs
 PMobj

}
class(PNOF) <- "PM"

#' @rdname PerformanceMetric
#' @export
LTY <- function(MSEobj=NULL, Ref=0.5, Yrs=-10) {
 Yrs <- ChkYrs(Yrs, MSEobj)
 PMobj <- new("PMobj")
 PMobj@Name <- paste0("Average Yield relative to Reference Yield (Years ", Yrs[1], "-", Yrs[2], ")")
 if (Ref != 1) {
   PMobj@Caption <- paste0('Prob. Yield > ', Ref, ' Ref. Yield (Years ', Yrs[1], "-", Yrs[2], ")")
 } else {
   PMobj@Caption <- paste0('Prob. Yield > Ref. Yield (Years ', Yrs[1], "-", Yrs[2], ")")
 }

 RefYd <- array(MSEobj@OM$RefY, dim=dim(MSEobj@Catch[,,Yrs[1]:Yrs[2]]))

 PMobj@Stat <- MSEobj@Catch[,,Yrs[1]:Yrs[2]]/RefYd
 PMobj@Ref <- Ref
 PMobj@Prob <- calcProb(PMobj@Stat > PMobj@Ref, MSEobj)

 PMobj@Mean <- calcMean(PMobj@Prob) # calculate mean probability by MP
 PMobj@MPs <- MSEobj@MPs
 PMobj

}
class(LTY) <- "PM"

#' @rdname PerformanceMetric
#' @export
STY <- LTY
formals(STY)$Yrs <- 10
class(STY) <- "PM"



#' @rdname PerformanceMetric
#' @export
Yield <- function(MSEobj=NULL, Ref=1, Yrs=NULL) {
 Yrs <- ChkYrs(Yrs, MSEobj)
 PMobj <- new("PMobj")
 PMobj@Name <- paste0("Yield relative to Reference Yield (Years ", Yrs[1], "-", Yrs[2], ")")
 PMobj@Caption <- paste0("Mean Relative Yield (Years ", Yrs[1], "-", Yrs[2], ")")

 RefYd <- array(MSEobj@OM$RefY, dim=dim(MSEobj@Catch[,,Yrs[1]:Yrs[2]]))

 PMobj@Stat <- MSEobj@Catch[,,Yrs[1]:Yrs[2]]/RefYd
 PMobj@Ref <- Ref
 PMobj@Prob <- calcProb(PMobj@Stat, MSEobj) # no probability to calculate

 PMobj@Mean <- calcMean(PMobj@Prob) # calculate mean probability by MP
 PMobj@MPs <- MSEobj@MPs
 PMobj

}
class(Yield) <- "PM"

#' @rdname PerformanceMetric
#' @export
AAVY <- function(MSEobj=NULL, Ref=0.2, Yrs=NULL) {
 Yrs <- ChkYrs(Yrs, MSEobj)
 PMobj <- new("PMobj")
 PMobj@Name <- paste0("Average Annual Variability in Yield (Years ", Yrs[1], "-", Yrs[2], ")")
 PMobj@Caption <- paste0('Prob. AAVY < ', Ref*100, "% (Years ", Yrs[1], "-", Yrs[2], ")")

 y1<- Yrs[1]:(Yrs[2]-1) # year index
 y2<-(Yrs[1]+1):Yrs[2]

 if (MSEobj@nMPs > 1) {
   AAVY <- apply(((((MSEobj@Catch[, , y1] - MSEobj@Catch[, , y2])/MSEobj@Catch[, , y2])^2)^0.5), c(1, 2), mean)
   # AAVY <- apply(((MSEobj@Catch[,,y1]-MSEobj@Catch[,,y2])^2)^0.5,c(1,2),mean)/apply(MSEobj@Catch[,,y2],c(1,2),mean)
 } else {
   AAVY <- array(apply(((((MSEobj@Catch[,1,y1]-MSEobj@Catch[,1,y2])/MSEobj@Catch[,1,y2])^2)^0.5),1,mean))
 }

 PMobj@Stat <- AAVY
 PMobj@Ref <- Ref
 PMobj@Prob <- calcProb(PMobj@Stat < Ref, MSEobj)  # probability AAVY < 0.2

 PMobj@Mean <- calcMean(PMobj@Prob) # calculate mean probability by MP
 PMobj@MPs <- MSEobj@MPs
 PMobj

}
class(AAVY) <- "PM"

#' @rdname PerformanceMetric
#' @export
AAVE <- function(MSEobj=NULL, Ref=0.2, Yrs=NULL) {
 Yrs <- ChkYrs(Yrs, MSEobj)
 PMobj <- new("PMobj")
 PMobj@Name <- paste0("Average Annual Variability in Effort (Years ", Yrs[1], "-", Yrs[2], ")")
 PMobj@Caption <- paste0('Prob. AAVE < ', Ref*100, "% (Years ", Yrs[1], "-", Yrs[2], ")")

 y1<- Yrs[1]:(Yrs[2]-1) # year index
 y2<-(Yrs[1]+1):Yrs[2]

 if (MSEobj@nMPs > 1) {
   AAVE <- apply(((((MSEobj@Effort[, , y1] - MSEobj@Effort[, , y2])/MSEobj@Effort[, , y2])^2)^0.5), c(1, 2), mean)
 } else {
   AAVE <- array(apply(((((MSEobj@Effort[,1,y1]-MSEobj@Effort[,1,y2])/MSEobj@Effort[,1,y2])^2)^0.5),1,mean))
 }

 PMobj@Stat <- AAVE
 PMobj@Ref <- Ref
 PMobj@Prob <- calcProb(PMobj@Stat < Ref, MSEobj)  # probability AAVE < 0.2

 PMobj@Mean <- calcMean(PMobj@Prob) # calculate mean probability by MP
 PMobj@MPs <- MSEobj@MPs
 PMobj

}
class(AAVE) <- "PM"




#' Calculate Probability
#'
#' @param PM A PM method
#'
#' @export
#' @keywords internal
#'
calcProb <- function(PM, MSEobj) {
 if (MSEobj@nMPs > 1) {
   mar <- 2
 } else mar <- 1
 mar <- 1:mar
 apply(PM, mar, mean, na.rm=TRUE)
}


#' Calculate Mean Probability
#'
#' @param Prob Prob slot from an object of class PMobj
#'
#' @export
#' @keywords internal
#'
calcMean <- function(Prob) {
 if ('matrix' %in% class(Prob)) return(apply(Prob , 2, mean, na.rm=TRUE))
 if ('numeric' %in% class(Prob)) return(mean(Prob, na.rm=TRUE))
}