diff --git a/.github/workflows/r-cmd-check.yml b/.github/workflows/r-cmd-check.yml index 678ea96..944cd77 100644 --- a/.github/workflows/r-cmd-check.yml +++ b/.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: diff --git a/DESCRIPTION b/DESCRIPTION index 932accf..1f8bc95 100644 --- a/DESCRIPTION +++ b/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 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) diff --git a/R/get_settings.R b/R/get_settings.R index 15eb553..9cec654 100644 --- a/R/get_settings.R +++ b/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.") diff --git a/R/get_settings_profile.R b/R/get_settings_profile.R index 5737c34..d0a670f 100644 --- a/R/get_settings_profile.R +++ b/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) diff --git a/R/get_summary.R b/R/get_summary.R index 44ab6e9..aa1ca7d 100644 --- a/R/get_summary.R +++ b/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") { diff --git a/R/profile_plot.R b/R/profile_plot.R index f14858e..53df251 100644 --- a/R/profile_plot.R +++ b/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_") ) diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index 15594a8..fe021b9 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -1,22 +1,34 @@ -#' Run [r4ss::profile] based on `model_settings` +#' Run [r4ss::profile()] based on `model_settings` #' -#' Generate likelihood profiles -#' To be called from the run_diagnostics function after creating -#' the model settings using the get_settings function. -#' -#' -#' @seealso The following functions interact with `profile_wrapper`: -#' * [run_diagnostics]: calls `profile_wrapper` -#' * [r4ss::profile]: the workhorse of `profile_wrapper` that does the parameter profiles +#' Setting up the specifications, running the profile using [r4ss::profile()], +#' and generating figures and tables can be tedious, error prone, and time +#' consuming. Thus, `profile_wrapper()` aims to further decrease the work +#' needed to generate a profile that is easily included in a assessment report +#' for the Pacific Fisheries Management Council. See the See Also section for +#' information on a workflow to use this function. #' +#' @seealso +#' The following functions interact with `profile_wrapper()`: +#' * [get_settings()]: populates a list of settings for `model_settings` +#' * [run_diagnostics()]: calls `profile_wrapper()` +#' * [r4ss::profile()]: the workhorse of `profile_wrapper()` that does the +#' parameter profiles #' #' @template mydir #' @template model_settings #' -#' @author Chantel Wetzel. +#' @author Chantel Wetzel and Ian Taylor. #' @return -#' Nothing is explicitly returned from `profile_wrapper`. -#' The following objects are saved to the disk. +#' Nothing is explicitly returned from `profile_wrapper()`. +#' The following objects are saved to the disk: +#' * `*_profile_output.Rdata` +#' * `*_trajectories_*` from [r4ss::SSplotComparisons()] +#' * `piner_panel_*.png` +#' * `parameter_panel_*.png` +#' * `run_diag_warning.txt` +#' * a copy of the control file saved to `model_settings$newctlfile` +#' * `backup_oldctlfile.ss` +#' * `backup_ss.par` #' @export profile_wrapper <- function(mydir, model_settings) { @@ -37,10 +49,9 @@ profile_wrapper <- function(mydir, model_settings) { for (aa in 1:N) { para <- model_settings$profile_details$parameters[aa] - prior_used <- model_settings$profile_details$use_prior_like[aa] # Create a profile folder with the same naming structure as the base model # Add a label to show if prior was used or not - profile_dir <- file.path(mydir, paste0(model_settings$base_name, "_profile_", para, "_prior_like_", prior_used)) + profile_dir <- file.path(mydir, paste0(model_settings$base_name, "_profile_", para)) dir.create(profile_dir, showWarnings = FALSE) # Check for existing files and delete @@ -56,21 +67,27 @@ profile_wrapper <- function(mydir, model_settings) { ), file = "run_diag_warning.txt") message(paste0("Running profile for ", para, ".")) - # Copy the control file to run from the copy - if (!file.exists(file.path(profile_dir, "control.ss_new"))) { - orig_dir <- getwd() - setwd(profile_dir) - r4ss::run( - dir = profile_dir, - exe = model_settings$exe, - extras = model_settings$extras - ) - setwd(orig_dir) + # check for whether oldctlfile exists + if (!file.exists(file.path(profile_dir, model_settings$oldctlfile))) { + # if the oldctlfile is control.ss_new, and doesn't exist, + # run the model to create it + if (model_settings$oldctlfile == "control.ss_new") { + if (model_settings$verbose) { + message("running model to get control.ss_new file") + } + r4ss::run( + dir = profile_dir, + exe = model_settings$exe, + extras = model_settings$extras + ) + } else { + stop("Can not find ", model_settings$oldctlfile) + } } - # Use the SS_parlines funtion to ensure that the input parameter can be found + # Use the SS_parlines function to ensure that the input parameter can be found check_para <- r4ss::SS_parlines( - ctlfile = "control.ss_new", + ctlfile = model_settings$oldctlfile, dir = profile_dir, verbose = FALSE, version = model_settings$version, @@ -79,16 +96,18 @@ profile_wrapper <- function(mydir, model_settings) { if (sum(check_para) == 0) { print(para) - stop(paste0("The input profile_custom does not match a parameter in the control.ss_new file.")) + stop("The input profile_custom does not match a parameter in the file ", + model_settings$oldctlfile) } - file.copy(file.path(profile_dir, "control.ss_new"), file.path(profile_dir, model_settings$newctlfile)) + # Copy oldctlfile to newctlfile before modifying it + file.copy(file.path(profile_dir, model_settings$oldctlfile), + file.path(profile_dir, model_settings$newctlfile)) + # Change the control file name in the starter file starter <- r4ss::SS_readstarter(file = file.path(profile_dir, "starter.ss")) - starter$ctlfile <- "control_modified.ss" + starter$ctlfile <- model_settings$newctlfile starter$init_values_src <- model_settings$init_values_src - # make sure the prior likelihood is calculated for non-estimated quantities - starter$prior_like <- prior_used r4ss::SS_writestarter(mylist = starter, dir = profile_dir, overwrite = TRUE) # Read in the base model @@ -141,26 +160,60 @@ profile_wrapper <- function(mydir, model_settings) { vec <- c(low, high) num <- sort(vec, index.return = TRUE)$ix - profile <- r4ss::profile( - dir = profile_dir, - oldctlfile = "control.ss_new", - newctlfile = model_settings$newctlfile, - linenum = model_settings$linenum, - string = para, - profilevec = vec, - usepar = model_settings$usepar, - globalpar = model_settings$globalpar, - parlinenum = model_settings$parlinenum, - parstring = model_settings$parstring, - saveoutput = model_settings$saveoutput, - overwrite = model_settings$overwrite, - whichruns = model_settings$whichruns, - prior_check = model_settings$prior_check, - read_like = model_settings$read_like, - exe = model_settings$exe, - verbose = model_settings$verbose, - extras = model_settings$extras - ) + # backup original control.ss_new file for use in second half of profile + file.copy(file.path(profile_dir, model_settings$oldctlfile), + file.path(profile_dir, "backup_oldctlfile.ss"), + overwrite = model_settings$overwrite) + # backup original par file for use in second half of profile + # if usepar = TRUE + file.copy(file.path(profile_dir, "ss.par"), + file.path(profile_dir, "backup_ss.par"), + overwrite = model_settings$overwrite) + + # loop over down, then up + for (iprofile in 1:2) { + whichruns <- which(vec %in% if(iprofile == 1){low} else {high}) + #if (iprofile == 1) { + # whichruns <- which(vec %in% low) + # if (!is.null(model_settings$whichruns)) { + # whichruns <- intersect(model_settings$whichruns, whichruns) + # } + #} + #if (iprofile == 2) { + # whichruns <- which(vec %in% high) + if (!is.null(model_settings$whichruns)) { + whichruns <- intersect(model_settings$whichruns, whichruns) + } + if (iprofile == 2) { + # copy backup back to use in second half of profile + file.copy(file.path(profile_dir, "backup_oldctlfile.ss"), + file.path(profile_dir, model_settings$oldctlfile)) + # copy backup back to use in second half of profile + file.copy(file.path(profile_dir, "backup_ss.par"), + file.path(profile_dir, "ss.par"), + overwrite = model_settings$overwrite) + } + profile <- r4ss::profile( + dir = profile_dir, + oldctlfile = model_settings$oldctlfile, + newctlfile = model_settings$newctlfile, + linenum = model_settings$linenum, + string = para, + profilevec = vec, + usepar = model_settings$usepar, + globalpar = model_settings$globalpar, + parlinenum = model_settings$parlinenum, + parstring = model_settings$parstring, + saveoutput = model_settings$saveoutput, + overwrite = model_settings$overwrite, + whichruns = whichruns, # values set above + prior_check = model_settings$prior_check, + read_like = model_settings$read_like, + exe = model_settings$exe, + verbose = model_settings$verbose, + extras = model_settings$extras + ) + } # Save the output and the summary name <- paste0("profile_", para) @@ -169,7 +222,12 @@ profile_wrapper <- function(mydir, model_settings) { profilemodels <- r4ss::SSgetoutput(dirvec = profile_dir, keyvec = num) profilesummary <- r4ss::SSsummarize(biglist = profilemodels) - + if(!is.null(model_settings$btarg)){ + profilesummary$btarg <- model_settings$btarg + profilesummary$minbthresh <- model_settings$minbthresh + } + profilesummary$subplots <- model_settings$subplots + profile_output <- list() profile_output$mydir <- profile_dir profile_output$para <- para diff --git a/R/rerun_profile_vals.R b/R/rerun_profile_vals.R index a9efc94..c5e9a06 100644 --- a/R/rerun_profile_vals.R +++ b/R/rerun_profile_vals.R @@ -15,7 +15,7 @@ #' @author Chantel Wetzel. #' @export #' -#' @example +#' @examples #' rerun_profile_vals(mydir = file.path(directory, "base_model"), #' model_settings = model_settings, #' para_name = "NatM_uniform_Fem_GP_1", @@ -46,23 +46,25 @@ rerun_profile_vals <- function(mydir, dir.create(temp_dir, showWarnings = FALSE) file.copy(file.path(profile_dir, "ss.exe"), temp_dir, overwrite = TRUE) - file.copy(file.path(profile_dir, "control_modified.ss"), temp_dir) - file.copy(file.path(profile_dir, "control.ss_new"), temp_dir) + file.copy(file.path(profile_dir, model_settings$newctlfile), temp_dir) + file.copy(file.path(profile_dir, model_settings$oldctlfile), temp_dir) file.copy(file.path(profile_dir, data_file_nm), temp_dir) file.copy(file.path(profile_dir, "starter.ss_new"), temp_dir) file.copy(file.path(profile_dir, "forecast.ss_new"), temp_dir) file.rename(file.path(temp_dir, "starter.ss_new"), file.path(temp_dir, "starter.ss")) file.rename(file.path(temp_dir, "forecast.ss_new"), file.path(temp_dir, "forecast.ss")) - # Use the SS_parlines funtion to ensure that the input parameter can be found + # Use the SS_parlines function to ensure that the input parameter can be found check_para <- r4ss::SS_parlines( + ctlfile = model_settings$oldctlfile, dir = temp_dir, verbose = FALSE, active = FALSE )$Label == para if (sum(check_para) == 0) { - stop(paste0("The input profile_custom does not match a parameter in the control.ss_new file.")) + stop("The input profile_custom does not match a parameter in the ", + model_settings$oldctlfile) } load(file.path(profile_dir, paste0(para_name, "_profile_output.Rdata"))) @@ -80,8 +82,8 @@ rerun_profile_vals <- function(mydir, for (i in run_num) { setwd(temp_dir) r4ss::SS_changepars( - ctlfile = "control_modified.ss", - newctlfile = "control_modified.ss", + ctlfile = model_settings$newctlfile, + newctlfile = model_settings$newctlfile, strings = para, newvals = vec[i], estimate = FALSE diff --git a/R/retro_wrapper.R b/R/retro_wrapper.R index 9cf7268..1069804 100644 --- a/R/retro_wrapper.R +++ b/R/retro_wrapper.R @@ -145,6 +145,8 @@ retro_wrapper <- function(mydir, model_settings) { ifelse(abs(model_settings$retro_yrs) == 1, "", "s") ) ), + btarg = model_settings$btarg, + minbthresh = model_settings$minbthresh, plotdir = retro_dir, legendloc = "topright", print = TRUE, @@ -158,6 +160,8 @@ retro_wrapper <- function(mydir, model_settings) { endyrvec = endyrvec, legendloc = "topleft", plotdir = retro_dir, + btarg = model_settings$btarg, + minbthresh = model_settings$minbthresh, print = TRUE, plot = FALSE, pdf = FALSE ), subplot = c(2, 4), diff --git a/README.md b/README.md index 9e27347..372105a 100644 --- a/README.md +++ b/README.md @@ -2,12 +2,14 @@ # nwfscDiag: Diagnostic Package for West Coast Groundfish Assessments -The package provides the functionality to conduct model diagnostics for SS models. The standard diagnostic included in this package are standard required analysis for West Coast Groundfish stock assessments. The package was designed to perform model diagnostics and create plots and tables in a standardized format. The standardized approach will facilitate the use of these outputs in the assessment template approach [sa4ss](https://github.com/pfmc-assessments/sa4ss). +The package provides the functionality to conduct model diagnostics for Stock Synthesis (SS3) models. The standard diagnostic included in this package are standard required analysis for U.S. West Coast Groundfish stock assessments managed by the Pacific Fisheries Management Council. The package was designed to perform model diagnostics and create plots and tables in a standardized format. The standardized approach will facilitate the use of these outputs in the assessment template approach [sa4ss](https://github.com/pfmc-assessments/sa4ss). The diagnostics created by the package are: -- jitter runs to ensure model convergence at the MLE -- retrospective runs to examine model sensitivity to recent data -- likelihood profiles across parameters +- jitter runs to ensure model convergence at the MLE, +- retrospective runs to examine model sensitivity to recent data, and +- likelihood profiles across parameters. + +This package does not maintain backward compatibility with previous versions of Stock Synthesis. However, if needed user can download older package versions that may work with older versions (3.30.+) of Stock Synthesis. ## Installation @@ -17,12 +19,12 @@ install.packages("remotes") remotes::install_github("pfmc-assessments/nwfscDiag") ``` ## Running the code -The package depends upon a few other packages and they should be installed upon installation of the package. The dependant packages are: +The package depends upon a few other packages and they should be installed upon installation of the package. The dependent packages are: ``` install.packages('dplyr') remotes::install_github('r4ss/r4ss') ``` -A new version of r4ss package was released on July 29, 2022. This release included some significant changes that included changes to function names and function inputs. The current version of the nwfscDiag 1.1.1 package is designed to work with the latest release of r4ss. Please see release version 1.0.1 to use earlier versions of r4ss. +A new version of r4ss package was released on July 29, 2022. This release included some significant changes that included changes to function names and function inputs. The current version of the nwfscDiag 1.1.2 package is designed to work with the latest release of r4ss. Please see release version 1.0.1 to use earlier versions of r4ss. ## Reporting problems diff --git a/man/get_settings_profile.Rd b/man/get_settings_profile.Rd index 98c17af..8808007 100644 --- a/man/get_settings_profile.Rd +++ b/man/get_settings_profile.Rd @@ -5,12 +5,12 @@ \title{Get Default Settings For Profiles} \usage{ get_settings_profile( - parameters = c("NatM_p_1_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), + parameters = c("NatM_uniform_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), low = c(0.4, 0.25, -2), high = c(0.4, 1, 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() ) } \arguments{ @@ -26,11 +26,9 @@ get_settings_profile( 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.} -\item{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.} +\item{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.} } \value{ A matrix of low, high, and step size values for the default parameters @@ -59,22 +57,20 @@ parameter. A user can select any of the options for specifying a parameter range # Define each parameter in real space get_settings_profile( - parameters = c("NatM_p_1_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) + 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') ) # Example 2: Run a profile for natural mortality one with the prior likelihood and one without get_settings_profile( - parameters = c("NatM_p_1_Fem_GP_1", "NatM_p_1_Fem_GP_1"), + parameters = c("NatM_uniform_Fem_GP_1", "NatM_uniform_Fem_GP_1"), 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") ) } diff --git a/man/profile_wrapper.Rd b/man/profile_wrapper.Rd index 2a208f2..37ff168 100644 --- a/man/profile_wrapper.Rd +++ b/man/profile_wrapper.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/profile_wrapper.R \name{profile_wrapper} \alias{profile_wrapper} -\title{Run \link[r4ss:profile]{r4ss::profile} based on \code{model_settings}} +\title{Run \code{\link[r4ss:profile]{r4ss::profile()}} based on \code{model_settings}} \usage{ profile_wrapper(mydir, model_settings) } @@ -17,21 +17,36 @@ Typically, created using \link{get_settings} but users can create their own list as long as it has all of the necessary components for the function at hand.} } \value{ -Nothing is explicitly returned from \code{profile_wrapper}. -The following objects are saved to the disk. +Nothing is explicitly returned from \code{profile_wrapper()}. +The following objects are saved to the disk: +\itemize{ +\item \verb{*_profile_output.Rdata} +\item \verb{*_trajectories_*} from \code{\link[r4ss:SSplotComparisons]{r4ss::SSplotComparisons()}} +\item \code{piner_panel_*.png} +\item \code{parameter_panel_*.png} +\item \code{run_diag_warning.txt} +\item a copy of the control file saved to \code{model_settings$newctlfile} +\item \code{backup_oldctlfile.ss} +\item \code{backup_ss.par} +} } \description{ -Generate likelihood profiles -To be called from the run_diagnostics function after creating -the model settings using the get_settings function. +Setting up the specifications, running the profile using \code{\link[r4ss:profile]{r4ss::profile()}}, +and generating figures and tables can be tedious, error prone, and time +consuming. Thus, \code{profile_wrapper()} aims to further decrease the work +needed to generate a profile that is easily included in a assessment report +for the Pacific Fisheries Management Council. See the See Also section for +information on a workflow to use this function. } \seealso{ -The following functions interact with \code{profile_wrapper}: +The following functions interact with \code{profile_wrapper()}: \itemize{ -\item \link{run_diagnostics}: calls \code{profile_wrapper} -\item \link[r4ss:profile]{r4ss::profile}: the workhorse of \code{profile_wrapper} that does the parameter profiles +\item \code{\link[=get_settings]{get_settings()}}: populates a list of settings for \code{model_settings} +\item \code{\link[=run_diagnostics]{run_diagnostics()}}: calls \code{profile_wrapper()} +\item \code{\link[r4ss:profile]{r4ss::profile()}}: the workhorse of \code{profile_wrapper()} that does the +parameter profiles } } \author{ -Chantel Wetzel. +Chantel Wetzel and Ian Taylor. } diff --git a/man/rerun_profile_vals.Rd b/man/rerun_profile_vals.Rd index 3d2c944..d80e4b3 100644 --- a/man/rerun_profile_vals.Rd +++ b/man/rerun_profile_vals.Rd @@ -31,6 +31,14 @@ you would like to rerun.} Generate likelihood profiles To be called from the run_diagnostics function after creating the model settings using the get_settings function. +} +\examples{ +rerun_profile_vals(mydir = file.path(directory, "base_model"), + model_settings = model_settings, + para_name = "NatM_uniform_Fem_GP_1", + run_num = c(3, 4), + data_file_nm = "data.ss") + } \author{ Chantel Wetzel. diff --git a/tests/testthat/test-runs.R b/tests/testthat/test-runs.R index 1857c2e..9ce6d53 100644 --- a/tests/testthat/test-runs.R +++ b/tests/testthat/test-runs.R @@ -52,7 +52,7 @@ test_that("Do profile using the simple model", { file.exists( file.path(path, paste0(model_settings$base_name, "_profile_", - profile_para[a], "_prior_like_", get[a, "use_prior_like"]), + profile_para[a]), paste0("profile_", profile_para[a], "_quant_table.csv"))) ) }