Merge pull request #24 from pfmc-assessments/profile_control

Profile control

.github/workflows/r-cmd-check.yml

@@ -20,7 +20,7 @@ on:
   pull_request:
     # specifying branches means the workflow will only run when the pull request is to the merge into the branch listed, in this case, main.
     branches:
-      -main
+      - main
 # jobs specifies the jobs to run. you can have multiple jobs per github actions workflow, but in this case we have only one, which is
 # named run-R-CMD-check.
 jobs:

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: nwfscDiag
 Type: Package
 Title: Generate standard NWFSC assessment diagnostics
-Version: 1.1.1
+Version: 1.1.2
 Author: Chantel Wetzel
 Maintainer: Chantel Wetzel <[email protected]>
 Description: Package that can automates diagnositics for SS3 models by running jitters, retrospective, and profiles.  
@@ -24,6 +24,6 @@ Suggests:
     devtools
 Remotes: 
     github::r4ss/r4ss
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.3
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)

R/get_settings.R

@@ -17,6 +17,7 @@
 #' settings <- list()
 #'
 get_settings <- function(settings = NULL, verbose = FALSE) {
+
   if (is.vector(settings)) settings <- as.list(settings)
 
   Settings_all <- list(
@@ -42,9 +43,14 @@ get_settings <- function(settings = NULL, verbose = FALSE) {
     retro_yrs = -1:-5,
     overwrite = TRUE,
     show_in_console = FALSE,
+
+    # Plot target lines
+    btarg = NULL, 
+    minbthresh = NULL,
 
     # Profile Settings
     remove_files = TRUE,
+    oldctlfile = "control.ss_new",
     newctlfile = "control_modified.ss",
     linenum = NULL,
     string = NULL,
@@ -58,7 +64,8 @@ get_settings <- function(settings = NULL, verbose = FALSE) {
     whichruns = NULL,
     prior_check = FALSE,
     read_like = TRUE,
-    init_values_src = 0
+    init_values_src = 0,
+    subplots = c(1, 3)
   )
 
   Settings_all$profile_details <- get_settings_profile()
@@ -80,7 +87,7 @@ get_settings <- function(settings = NULL, verbose = FALSE) {
     if (!is.numeric(Settings_all$profile_details$low) &
       !is.numeric(Settings_all$profile_details$high) &
       !is.numeric(Settings_all$profile_details$step_size)) {
-      stop("There is a non-numeric value in the low, high, or step size fiedl of the get_settings_profile data frame.")
+      stop("There is a non-numeric value in the low, high, or step size field of the get_settings_profile data frame.")
     }
     if (sum(!Settings_all$profile_details$param_space %in% c("real", "relative", "multiplier")) > 0) {
       stop("The param_space column should be either real or relative in the get_settings_profile data frame.")

R/get_settings_profile.R

@@ -21,11 +21,9 @@
 #' @param param_space options: real, mulitplier, relative indicates how to interpret the low and high bound values.
 #' real indicates bounds in the parameter space, relative indicates how far to go from the base parameter, and
 #' multiplier indicates that low and high bounds are set at x\% above and below the base parameter.
-#' @param use_prior_like options: Option to include or exclude the prior likelihoods in the likelihood
-#' profiles. A value of 0 corresponds to the Stock Sythesis were this would exclude
-#' and a value of 1 would include the prior likelihood contribution. Parameters with priors used for estimation
-#' (e.g., natural mortality, steepness) are often profiled across and including or excluding the prior likelihood
-#' contribution may be wanted in specific instances.  The default setting excludes the prior likelihood contributions.
+#' @param use_prior_like Deprecated: The use_prior_like input is no longer needed since r4ss now 
+#' automatically plots the likelihood profile with and without any parameter prior likelihood contributions
+#' regardless of the setting in the user starter file. 
 #'
 #' @return A matrix of low, high, and step size values for the default parameters
 #' that should be profiled. The goal is to provide users with a template
@@ -41,11 +39,10 @@
 #' # Define each parameter in real space
 #' get_settings_profile(
 #'   parameters = c("NatM_uniform_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"),
-#'   # 						 low =  c(0.02, 0.25, 8),
-#'   # 						 high = c(0.07, 1.0,  11),
-#'   # 						 step_size = c(0.005, 0.05, 0.25),
-#'   # 						 param_space = c('real', 'real', 'real'),
-#'   use_prior_like = c(1, 1, 0)
+#'    						 low =  c(0.02, 0.25, 8),
+#'    						 high = c(0.07, 1.0,  11),
+#'    						 step_size = c(0.005, 0.05, 0.25),
+#'    						 param_space = c('real', 'real', 'real')
 #' )
 #'
 #' # Example 2: Run a profile for natural mortality one with the prior likelihood and one without
@@ -54,8 +51,7 @@
 #'   low = c(0.40, 0.40),
 #'   high = c(0.40, 0.40),
 #'   step_size = c(0.005, 0.005),
-#'   param_space = c("multiplier", "multiplier"),
-#'   use_prior_like = c(0, 1)
+#'   param_space = c("multiplier", "multiplier")
 #' )
 #' }
 #'
@@ -64,22 +60,29 @@ get_settings_profile <- function(parameters = c("NatM_uniform_Fem_GP_1", "SR_BH_
                                  high = c(0.40, 1.0, 2),
                                  step_size = c(0.01, 0.05, 0.25),
                                  param_space = c("multiplier", "real", "relative"),
-                                 use_prior_like = c(0, 0, 0)) {
+                                 use_prior_like = lifecycle::deprecated()) 
+{
+
   if (length(parameters) != length(low) |
     length(parameters) != length(high) |
     length(parameters) != length(step_size) |
-    length(parameters) != length(param_space) |
-    length(parameters) != length(use_prior_like)) {
+    length(parameters) != length(param_space)) {
     stop("Error: input vectors do match in length.")
   }
 
+  if (lifecycle::is_present(use_prior_like)) {
+    lifecycle::deprecate_warn(
+      when = "1.1.2",
+      what = "get_settings_profile(use_prior_like)"
+    )
+  }
+
   out <- data.frame(
     parameters = parameters,
     low = low,
     high = high,
     step_size = step_size,
-    param_space = param_space,
-    use_prior_like = use_prior_like
+    param_space = param_space
   )
 
   return(out)

R/get_summary.R

@@ -15,7 +15,7 @@
 
 get_summary <- function(mydir, para, vec, name, profilemodels, profilesummary) {
 
-  # Need to identify a way to determine if a model estmates male growth parameters as offsets from females
+  # Need to identify a way to determine if a model estimates male growth parameters as offsets from females
 
   # get output
   outputs <- profilemodels
@@ -30,8 +30,8 @@ get_summary <- function(mydir, para, vec, name, profilemodels, profilesummary) {
     "likelihood" = sapply(sapply(outputs, "[[", "likelihoods_used", simplify = FALSE), "[", 1, 1),
     "gradient" = sapply(outputs, "[[", "maximum_gradient_component"),
     "SB0" = sapply(quants, "[[", "SSB_Virgin", "Value"),
-    "SBfinal" = sapply(quants, "[[", paste0("SSB_", outputs[[1]]$endyr), "Value"),
-    "Deplfinal" = sapply(quants, "[[", paste0("Bratio_", outputs[[1]]$endyr), "Value"),
+    "SBfinal" = sapply(quants, "[[", paste0("SSB_", outputs[[1]]$endyr + 1), "Value"),
+    "Deplfinal" = sapply(quants, "[[", paste0("Bratio_", outputs[[1]]$endyr + 1), "Value"),
     # "Fmsy" = sapply(quants, "[[", "annF_MSY", "Value"),
     "Nparsonbounds" = apply(status, 2, function(x) sum(x %in% c("LO", "HI"))),
     stringsAsFactors = FALSE
@@ -57,17 +57,17 @@ get_summary <- function(mydir, para, vec, name, profilemodels, profilesummary) {
     parmlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Parm_devs", 1:n]),
     R0 = as.numeric(x$pars[x$pars$Label == "SR_LN(R0)", 1:n]),
     SB0 = as.numeric(x$SpawnBio[x$SpawnBio$Label == "SSB_Virgin", 1:n]),
-    SBfinal = as.numeric(x$SpawnBio[x$SpawnBio$Label == paste0("SSB_", endyr), 1:n]),
-    deplfinal = as.numeric(x$Bratio[x$Bratio$Label == paste0("Bratio_", endyr), 1:n]),
+    SBfinal = as.numeric(x$SpawnBio[x$SpawnBio$Label == paste0("SSB_", endyr + 1), 1:n]),
+    deplfinal = as.numeric(x$Bratio[x$Bratio$Label == paste0("Bratio_", endyr + 1), 1:n]),
     yieldspr = as.numeric(x$quants[x$quants$Label == "Dead_Catch_SPR", 1:n]),
     steep = as.numeric(x$pars[x$pars$Label == "SR_BH_steep", 1:n]),
-    mfem = as.numeric(x$pars[x$pars$Label == "NatM_p_1_Fem_GP_1", 1:n]),
+    mfem = as.numeric(x$pars[x$pars$Label == "NatM_uniform_Fem_GP_1", 1:n]),
     lminfem = as.numeric(x$pars[x$pars$Label == "L_at_Amin_Fem_GP_1", 1:n]),
     lmaxfem = as.numeric(x$pars[x$pars$Label == "L_at_Amax_Fem_GP_1", 1:n]),
     kfem = as.numeric(x$pars[x$pars$Label == "VonBert_K_Fem_GP_1", 1:n]),
     cv1fem = as.numeric(x$pars[grep("young_Fem_GP_1", x$pars$Label), 1:n]),
     cv2fem = as.numeric(x$pars[grep("old_Fem_GP_1", x$pars$Label), 1:n]),
-    mmale = as.numeric(x$pars[x$pars$Label == "NatM_p_1_Mal_GP_1", 1:n]),
+    mmale = as.numeric(x$pars[x$pars$Label == "NatM_uniform_Mal_GP_1", 1:n]),
     lminmale = as.numeric(x$pars[x$pars$Label == "L_at_Amin_Mal_GP_1", 1:n]),
     lmaxmale = as.numeric(x$pars[x$pars$Label == "L_at_Amax_Mal_GP_1", 1:n]),
     kmale = as.numeric(x$pars[x$pars$Label == "VonBert_K_Mal_GP_1", 1:n]),
@@ -81,10 +81,10 @@ get_summary <- function(mydir, para, vec, name, profilemodels, profilesummary) {
   if (para == "SR_LN(R0)") {
     colnames(new_out) <- paste0("R0 ", vec)
   }
-  if (para == "NatM_p_1_Fem_GP_1") {
+  if (para == "NatM_uniform_Fem_GP_1") {
     colnames(new_out) <- paste0("M_f ", vec)
   }
-  if (para == "NatM_p_1_Mal_GP_1") {
+  if (para == "NatM_uniform_Mal_GP_1") {
     colnames(new_out) <- paste0("M_m ", vec)
   }
   if (para == "SR_BH_steep") {

R/profile_plot.R

@@ -18,6 +18,7 @@
 #' @seealso [profile_wrapper] and [rerun_profile_vals] call `profile_plot`.
 
 profile_plot <- function(mydir, rep, para, profilesummary) {
+
   label <- ifelse(para == "SR_LN(R0)", expression(log(italic(R)[0])),
     ifelse(para %in% c("NatM_p_1_Fem_GP_1", "NatM_uniform_Fem_GP_1"), "Natural Mortality (female)",
       ifelse(para %in% c("NatM_p_1_Mal_GP_1", "NatM_uniform_Mal_GP_1"), "Natural Mortality (male)",
@@ -130,24 +131,30 @@ profile_plot <- function(mydir, rep, para, profilesummary) {
 
   dev.off()
 
-  maxyr <- min(profilesummary$endyrs)
+  maxyr <- min(profilesummary$endyrs + 1)
   minyr <- max(profilesummary$startyrs)
-
   est <- rep$parameters[rep$parameters$Label == para, "Value", 2]
   sb0_est <- rep$derived_quants[rep$derived_quants$Label == "SSB_Virgin", "Value"]
   sbf_est <- rep$derived_quants[rep$derived_quants$Label == paste0("SSB_", maxyr), "Value"]
   depl_est <- rep$derived_quants[rep$derived_quants$Label == paste0("Bratio_", maxyr), "Value"]
 
   x <- as.numeric(profilesummary$pars[profilesummary$pars$Label == para, n])
-  like <- as.numeric(profilesummary$likelihoods[profilesummary$likelihoods$Label == "TOTAL", n] - rep$likelihoods_used[1, 1])
+  # determine whether to include the prior likelihood component in the likelihood profile
+  starter <- r4ss::SS_readstarter(file = file.path(mydir, "starter.ss"))
+  like <- as.numeric(profilesummary$likelihoods[profilesummary$likelihoods$Label == "TOTAL", n] -
+    ifelse(starter$prior_like == 0, 
+      profilesummary$likelihoods[profilesummary$likelihoods$Label == "Parm_priors", n],
+      0) - 
+    rep$likelihoods_used[1, 1])
+
   ylike <- c(min(like) + ifelse(min(like) != 0, -0.5, 0), max(like))
   sb0 <- as.numeric(profilesummary$SpawnBio[na.omit(profilesummary$SpawnBio$Label) == "SSB_Virgin", n])
   sbf <- as.numeric(profilesummary$SpawnBio[na.omit(profilesummary$SpawnBio$Yr) == maxyr, n])
   depl <- as.numeric(profilesummary$Bratio[na.omit(profilesummary$Bratio$Yr) == maxyr, n])
 
   # Get the relative management targets - only grab the first element since the targets should be the same
-  btarg <- as.numeric(profilesummary$btargs[1])
-  thresh <- ifelse(btarg == 0.40, 0.25, 0.125)
+  btarg <- as.numeric(profilesummary$btarg[1])
+  thresh <- as.numeric(profilesummary$minbthresh[1]) # ifelse(btarg == 0.40, 0.25, ifelse(btarg == 0.25, 0.125, -1))    
 
   pngfun(wd = mydir, file = paste0("parameter_panel_", para, ".png"), h = 7, w = 7)
   par(mfrow = c(2, 2), mar = c(4, 4, 2, 2), oma = c(1, 1, 1, 1))
@@ -163,11 +170,13 @@ profile_plot <- function(mydir, rep, para, profilesummary) {
   # parameter vs. final depletion
   plot(x, depl, type = "l", lwd = 2, xlab = label, ylab = "Fraction of unfished", ylim = c(0, 1.2))
   points(est, depl_est, pch = 21, col = "black", bg = "blue", cex = 1.5)
-  abline(h = c(thresh, btarg), lty = c(2, 2), col = c("red", "darkgreen"))
-  legend("bottomright",
-    legend = c("Management target", "Minimum stock size threshold"),
-    lty = 2, col = c("red", "darkgreen"), bty = "n"
-  )
+  abline(h = c(btarg, thresh), lty = c(2, 2), col = c("darkgreen", "red"))
+  if(btarg > 0){
+    legend("bottomright",
+           legend = c("Management target", "Minimum stock size threshold"),
+           lty = 2, col = c("darkgreen", "red"), bty = "n"
+    )    
+  }
 
   # parameter vs. SB0
   plot(x, sb0, type = "l", lwd = 2, xlab = label, ylab = expression(SB[0]), ylim = c(0, max(sb0)))
@@ -202,7 +211,10 @@ profile_plot <- function(mydir, rep, para, profilesummary) {
       ifelse(est == x, " (base)", "")
     ),
     ylimAdj = 1.15,
-    plotdir = mydir, subplots = c(1, 3),
+    btarg = btarg,
+    minbthresh = thresh,
+    plotdir = mydir, 
+    subplots = profilesummary$subplots,
     pdf = FALSE, print = TRUE, plot = FALSE,
     filenameprefix = paste0(para, "_trajectories_")
   )