R/define_S4_class_profile.R

In tcarruth/MSEtool: Management Strategy Evaluation Toolkit

#' Class-\code{prof}
#'
#' An S4 class that contains output from \link{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 \code{Par}) from the assessment.
#' @slot grid A data.frame of the change in negative log-likelihood (\code{nll}) based on the profile of the parameters.
#' @seealso \link{plot.prof} \link{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 \link{profile}. If a two-parameter profile is performed, then a contour
#' plot of the likelihood surface is returned.
#' @param x An object of class \linkS4class{prof} returned by \link{profile}.
#' @param contour_levels Integer, passed to \code{nlevels} argument of \link[graphics]{contour}.
#' @param ... Miscellaneous. Not used.
#' @author Q. Huynh
#' @importFrom reshape2 acast
#' @importFrom graphics contour
#' @exportMethod plot
setMethod("plot", signature(x = "prof", y = "missing"),
          function(x, contour_levels = 20, ...) {
            joint_profile <- length(x@MLE) == 2

            if(joint_profile) {
              z.mat <- 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, 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")
              plot(x = x.plot, y = y.plot, xlab = profile_par, ylab = "Change in neg. log-likelihood", typ = "l")
              abline(h = 0, col = "grey")
              points(getElement(x@grid, profile_par), getElement(x@grid, "nll"), 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
#' @export
setGeneric("profile", function(fitted, ...) standardGeneric("profile"))

#' @name profile
#' @title Profile likelihood of assessment models
#' @aliases profile,Assessment-method
#' @description Profile the likelihood for parameters of assessment models.
#'
#' @param fitted,Assessment An object of class \linkS4class{Assessment}.
#' @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: \code{R0} and \code{h}
#' \item SP and SP_SS: \code{FMSY} and \code{MSY}
#' \item DD and cDD_SS: \code{R0} and \code{h}
#' \item SCA and SCA_Pope: \code{R0} and \code{h}
#' \item SCA2: \code{meanR}
#' \item VPA: \code{F_term}
#' \item SSS: \code{R0}
#' }
#' @author Q. Huynh
#' @return An object of class \linkS4class{prof} that contains a data frame of negative log-likelihood values from the profile and, optionally,
#' a figure of the likelihood surface.
#' @examples
#' \donttest{
#' output <- DD_TMB(Data = DLMtool::Red_snapper)
#' pro <- profile(output, R0 = seq(0.75, 1.25, 0.025), h = seq(0.9, 0.99, 0.01))
#' pro <- profile(output, R0 = seq(0.75, 1.25, 0.025)) # Profile R0 only
#'
#' # 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) == 0) stop("No parameters for profile was found. See help.")

            f <- get(paste0('profile_likelihood_', fitted@Model))
            res <- f(fitted, ...)
            if(figure) plot(res)
            return(res)
          })