R/define_S4_class_profile.R

In SAMtool: Stock Assessment Methods Toolkit

#' Class-`prof`
#'
#' An S4 class that contains output from [profile].
#'
#' @name prof-class
#' @docType class
#'
#' @slot Model Name of the assessment model.
#' @slot Name Name of Data object.
#' @slot Par Character vector of parameters that were profiled.
#' @slot MLE Numeric vector of the estimated values of the parameters (corresponding to `Par`) from the assessment.
#' @slot grid A data.frame of the change in negative log-likelihood (`nll`) based on the profile of the parameters.
#' @seealso [plot.prof] [profile]
#' @author Q. Huynh
#' @export prof
#' @exportClass prof
prof <- setClass("prof", slots = c(Model = "character", Name = "character", Par = "character", MLE = "vector", grid = "data.frame"))


#' @name plot.prof
#' @aliases plot,prof,missing-method
#' @title Plot profile object
#' @description Generates a profile plot generated by [profile]. If a two-parameter profile is performed, then a contour
#' plot of the likelihood surface is returned.
#' @param x An object of class [prof-class] returned by [profile].
#' @param contour_levels Integer, passed to `nlevels` argument of [contour][graphics::contour].
#' @param ... Miscellaneous. Not used.
#' @return A likelihood profile plot, either a one-dimensional line plot or a two-dimensional contour plot.
#' @author Q. Huynh
#' @exportMethod plot
setMethod("plot", signature(x = "prof", y = "missing"),
          function(x, contour_levels = 20, ...) {
            joint_profile <- length(x@MLE) == 2
            if (joint_profile && !requireNamespace("reshape2", quietly = TRUE)) {
              stop("Please install the reshape2 package.", call. = FALSE)
            }
            
            if (x@Model == "RCM" && length(x@Par) == 1) {
              if (!requireNamespace("ggplot2", quietly = TRUE)) {
                stop("Please install the ggplot2 package.", call. = FALSE)
              }
              g <- parse(text = paste0('reshape2::melt(x@grid, id.vars = x@Par) %>% 
                mutate(Type = ifelse(grepl("Fleet", variable), "Fishery", 
                                     ifelse(grepl("Index", variable), "Index", as.character(variable)))) %>%
                group_by(', x@Par, ', Type) %>% summarise(value = sum(value, na.rm = TRUE), .groups = "drop") %>%
                group_by(Type) %>%
                mutate(value = value - min(value, na.rm = TRUE)) %>%
              ggplot(aes(', x@Par, ', value, colour = Type)) + geom_point() + geom_line() + theme_bw() + 
                labs(y = "Change in neg. log-likelihood")')) %>% eval()
              return(g)
                
            } else if (joint_profile) {
              z.mat <- reshape2::acast(x@grid, as.list(x@Par), value.var = "nll")
              x.mat <- as.numeric(dimnames(z.mat)[[1]])
              y.mat <- as.numeric(dimnames(z.mat)[[2]])
              contour(x = x.mat, y = y.mat, z = z.mat - min(z.mat, na.rm = TRUE), 
                      xlab = x@Par[1], ylab = x@Par[2], nlevels = contour_levels)
              points(x@MLE[1], x@MLE[2], col = "red", cex = 1.5, pch = 16)

              if (x@MLE[1] >= min(x.mat) && x@MLE[1] <= max(x.mat) && x@MLE[2] >= min(y.mat) && x@MLE[2] <= max(y.mat)) {
                sub <- "Red point indicates model estimate."
              } else sub <- NULL
              title("Contour plot of joint likelihood profile", sub = sub)

            } else {
              profile_par <- x@Par
              x.plot <- getElement(x@grid, profile_par)
              y.plot <- getElement(x@grid, "nll") - min(getElement(x@grid, "nll"), na.rm = TRUE)
              plot(x.plot, y.plot, xlab = profile_par, ylab = "Change in neg. log-likelihood", typ = "n")
              abline(h = 0, col = "grey")
              lines(x.plot, y.plot)
              points(x.plot, y.plot, pch = 16)
              abline(v = x@MLE, lty = 2)

              if (x@MLE >= min(x.plot) && x@MLE <= max(x.plot)) {
                sub <- "Dotted line indicates model estimate."
              } else sub <- NULL
              title(paste("Likelihood profile of", profile_par), sub = sub)
            }
            invisible()
          })


#' @rdname profile
#' @title Profile likelihood of assessment models
#' @aliases profile profile,Assessment-method
#' @description Profile the likelihood for parameters of assessment models.
#'
#' @param fitted,Assessment An object of class [Assessment-class].
#' @param figure Logical, indicates whether a figure will be plotted.
#' @param ... A sequence of values of the parameter(s) for the profile. See details and example below.
#' See details for name of arguments to be passed on.
#' @details
#'
#' For the following assessment models, possible sequence of values for profiling are:
#' \itemize{
#' \item DD_TMB and DD_SS: `R0` and `h`
#' \item SP and SP_SS: `FMSY` and `MSY`
#' \item DD and cDD_SS: `R0` and `h`
#' \item SCA and SCA_Pope: `R0` and `h`
#' \item SCA2: `meanR`
#' \item VPA: `F_term`
#' \item SSS: `R0`
#' }
#' 
#' For RCM: `D` (spawning biomass depletion), `R0`, and `h` are used. If the Mesnil-Rochet stock-recruit function
#' is used, can also profile `MRRmax` and `MRgamma`.
#' @author Q. Huynh
#' @return An object of class [prof-class] that contains a data frame of negative log-likelihood values from the profile and, optionally,
#' a figure of the likelihood surface.
#' @examples
#' \donttest{
#' output <- SCA(Data = MSEtool::SimulatedData)
#' 
#' # Profile R0 only
#' pro <- profile(output, R0 = seq(1000, 2000, 50))
#' 
#' # Profile both R0 and steepness
#' pro <- profile(output, R0 = seq(1000, 2000, 100), h = seq(0.8, 0.95, 0.025))
#'
#' # Ensure your grid is of proper resolution. A grid that is too coarse
#' # will likely distort the shape of the likelihood surface.
#' }
#' @exportMethod profile
setMethod("profile", signature(fitted = "Assessment"),
          function(fitted, figure = TRUE, ...) {
            dots <- list(...)
            if (!length(dots)) stop("No parameters for profile was found. See help.")
            dots[["Assessment"]] <- fitted

            func <- get(paste0('profile_likelihood_', fitted@Model))
            res <- do.call2(func, dots)
            if (figure) plot(res)
            return(res)
          })

#' @rdname profile
#' @exportMethod profile
setMethod("profile", signature(fitted = "RCModel"),
          function(fitted, figure = TRUE, ...) {
            dots <- list(...)
            if (!length(dots)) stop("No parameters for profile was found. See help.")
            dots[["Assessment"]] <- fitted
            
            res <- profile_likelihood_RCM(fitted, ...)
            if (figure) plot(res)
            return(res)
          })