Merge pull request #38 from pfmc-assessments/retro_figs

Add retrospective figures

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

@@ -13,9 +13,9 @@ on:
   # The push build trigger runs jobs when new commits are pushed up to github.
   push:
     # specifying branches means the workflow will only run on pushes to the branches listed, in this case, only main
-    branches: 
+    branches:
       - main
-      - r4ss_1.46.1
+      - r4ss_1.50.0
   # The pull_request build trigger runs jobs when a pull request is made or commits are pushed to the pull request.
   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.
@@ -52,14 +52,13 @@ jobs:
           extra-packages: any::rcmdcheck
           needs: check
 
-      - name: Get the latest release SS3 executable and move to expected location
+      - name: Get the latest SS3 executable and move to expected location
         run: |
-          curl https://api.github.com/repos/nmfs-ost/ss3-source-code/releases/latest | grep "browser_download_url" | grep -Eo 'https://[^\"]*' | grep "ss_linux" | xargs wget
-          mv ss_linux ss
-          sudo chmod a+x ss
-          cp ss inst/extdata/simple/ss
-          rm ss
-          
+          wget -O ss3 https://github.com/nmfs-ost/ss3-source-code/releases/latest/download/ss3_linux
+          sudo chmod a+x ss3
+          cp ss3 inst/extdata/simple_small/ss3
+          rm ss3
+
       - name: run r-cmd-check using R
         env:
           _R_CHECK_CRAN_INCOMING_REMOTE_: false

DESCRIPTION

@@ -12,6 +12,7 @@ LazyData: true
 Depends:
     R (>= 3.5)
 Imports:
+    cli,
     dplyr,
     kableExtra,
     knitr,
@@ -27,6 +28,6 @@ Suggests:
     testthat
 Remotes: 
     github::r4ss/r4ss
-RoxygenNote: 7.2.3
+RoxygenNote: 7.3.2
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)

NAMESPACE

@@ -1,15 +1,25 @@
 # Generated by roxygen2: do not edit by hand
 
 export("%>%")
+export(check_profile_range)
+export(get_jitter_quants)
+export(get_param_values)
+export(get_retro_quants)
 export(get_settings)
 export(get_settings_profile)
 export(get_summary)
 export(jitter_wrapper)
+export(plot_jitter)
+export(plot_profile)
+export(plot_retro)
 export(pngfun)
 export(pretty_decimal)
-export(profile_plot)
 export(profile_wrapper)
 export(rerun_profile_vals)
 export(retro_wrapper)
+export(round_any)
 export(run_diagnostics)
+export(run_jitter)
+export(run_profile)
+export(run_retro)
 importFrom(magrittr,"%>%")

R/check_profile_range.R

@@ -0,0 +1,73 @@
+#' Check the profile range specified
+#'
+#'
+#' @template mydir
+#' @template model_settings
+#'
+#' @author Chantel Wetzel
+#' @return
+#' Nothing is explicitly returned from `check_profile_range()`
+#' @export
+#'
+#'
+check_profile_range <- function(mydir, model_settings){
+  # Read in the base model
+  rep <- r4ss::SS_output(
+    file.path(mydir, model_settings$base_name),
+    covar = FALSE,
+    printstats = FALSE,
+    verbose = FALSE
+  )
+
+  N <- nrow(model_settings$profile_details)
+  for (aa in 1:N){
+    profile_details <- model_settings[["profile_details"]][aa, ]
+    para <- profile_details[, "parameters"]
+    est <- rep$parameters[rep$parameters$Label == para, "Value"]
+
+   # Determine the parameter range
+   if (profile_details$param_space == "relative") {
+     range <- c(
+       est + profile_details$low,
+       est + profile_details$high
+     )
+   }
+   if (profile_details$param_space == "multiplier") {
+     range <- c(
+       est - est * profile_details$low,
+       est + est * profile_details$high
+     )
+   }
+   if (profile_details$param_space == "real") {
+     range <- c(
+       profile_details$low,
+       profile_details$high
+     )
+   }
+   step_size <- profile_details$step_size
+
+   # Create parameter vect from base down and the base up
+   if (est != round_any(est, step_size, f = floor)) {
+     low <- rev(seq(
+       round_any(range[1], step_size, f = ceiling),
+       round_any(est, step_size, f = floor), step_size
+     ))
+   } else {
+     low <- rev(seq(
+       round_any(range[1], step_size, f = ceiling),
+       round_any(est, step_size, f = floor) - step_size, step_size
+     ))
+   }
+
+   if (est != round_any(est, step_size, f = ceiling)) {
+     high <- c(est, seq(round_any(est, step_size, f = ceiling), range[2], step_size))
+   } else {
+     high <- c(seq(round_any(est, step_size, f = ceiling), range[2], step_size))
+   }
+
+   vec <- c(low, high)
+   cli::cli_inform(
+     "Profiling over {para} across values of {vec}."
+   )
+  }
+}

R/get_jitter_quants.R

@@ -0,0 +1,71 @@
+#' Create output files with jitter quantities
+#'
+#' Tables generated to visualize results
+#'
+#'
+#' @template mydir
+#' @param output List of model output created by [run_jitter()].
+#' @template model_settings
+#'
+#' @author Chantel Wetzel
+#' @return
+#' Nothing is explicitly returned from [get_jitter_quants()].
+#'
+#'
+#' @export
+
+get_jitter_quants <- function(mydir, model_settings, output) {
+  jitter_dir <- output[["plotdir"]]
+  like <- output[["like"]]
+  est <- output[["est"]]
+  profilesummary <- output[["profilesummary"]]
+
+  outputs <- output$profilemodels
+  quants <- lapply(outputs, "[[", "derived_quants")
+  status <- sapply(sapply(outputs, "[[", "parameters", simplify = FALSE), "[[", "Status")
+  bounds <- apply(status, 2, function(x) rownames(outputs[[1]]$parameters)[x %in% c("LO", "HI")])
+  out <- data.frame(
+    "run" = gsub("replist", "", names(outputs)),
+    "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_", profilesummary$endyrs[1]), "Value"),
+    "Nparsonbounds" = apply(status, 2, function(x) sum(x %in% c("LO", "HI"))),
+    "Lowest NLL" = ifelse(min(like) == like, "Best Fit", 0),
+    stringsAsFactors = FALSE
+  )
+
+  # Write a md file to be included in a stock assessment document
+  # Text was pirated from @chantelwetzel-noaa's 2021 dover assessment
+  file_md <- file.path(jitter_dir, "model-results-jitter.md")
+  sink(file_md)
+  on.exit(sink(), add = TRUE)
+  cat(
+    sep = "",
+    "Model convergence was in part based on starting the minimization process ",
+    "from dispersed values of the maximum likelihood estimates to determine if the ",
+    "estimation routine results in a smaller likelihood.\n",
+    "Starting parameters were jittered using the built-in functionality of ",
+    "Stock Synthesis, where you specify a jitter fraction.\n",
+    "Here we used a jitter fraction of ",
+    round(model_settings$jitter_fraction, 2), " and the jittering was repeated ",
+    xfun::numbers_to_words(model_settings$Njitter), " times.\n",
+    "A better, i.e., lower negative log-likelihood, fit was ",
+    ifelse(
+      sum(like - est < 0) == 0,
+      "not found",
+      paste0("found for ", xfun::numbers_to_words(sum(like - est < 0)), " fits")
+    ), ".\n",
+    "Several models resulted in similar log-likelihood values ",
+    "with little difference in the overall model estimates, ",
+    "indicating a relatively flat likelihood surface around the maximum likelihood estimate.\n",
+    "Through the jittering analysis performed here and ",
+    "the estimation of likelihood profiles, ",
+    "we are confident that the base model as presented represents the ",
+    "best fit to the data given the assumptions made.\n"
+  )
+
+  # write tables
+  utils::write.csv(x = table(unlist(bounds)), file = file.path(jitter_dir, "jitter_parsonbounds.csv"), row.names = FALSE)
+  utils::write.csv(x = out, file = file.path(jitter_dir, "jitter_results.csv"), row.names = FALSE)
+}

R/get_param_values.R

@@ -0,0 +1,72 @@
+#' Generate likelihood profiles
+#' To be called from the run_diagnostics function after creating
+#' the model settings using the get_settings function.
+#'
+#'
+#' @template mydir
+#' @param para A character string specifying the SS3 parameter name that the
+#'   profile pertains to. The parameter name should match the name in the
+#'   control.ss_new file from SS3.
+#' @param vec Vector of parameter values or retrospective runs the summary object contains
+#' @param summary List created by the [r4ss::SSsummarize()] function
+#'
+#' @author Chantel Wetzel & Kelli Johnson
+#' @export
+
+get_param_values <- function(mydir, para = NULL, vec, summary) {
+
+  x <- summary
+  n <- x$n
+  endyr <- x$endyrs[1] + 1
+  out <- data.frame(
+    totlikelihood = as.numeric(x$likelihoods[x$likelihoods$Label == "TOTAL", 1:n]),
+    surveylike = as.numeric(x$likelihoods[x$likelihoods$Label == "Survey", 1:n]),
+    discardlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Discard", 1:n]),
+    lengthlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Length_comp", 1:n]),
+    agelike = as.numeric(x$likelihoods[x$likelihoods$Label == "Age_comp", 1:n]),
+    recrlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Recruitment", 1:n]),
+    forerecrlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Forecast_Recruitment", 1:n]),
+    priorlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Parm_priors", 1:n]),
+    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]),
+    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_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_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]),
+    cv1male = as.numeric(x$pars[grep("young_Mal_GP_1", x$pars$Label), 1:n]),
+    cv2male = as.numeric(x$pars[grep("old_Mal_GP_1", x$pars$Label), 1:n]),
+    stringsAsFactors = FALSE
+  )
+
+  out <- t(out)
+  colnames(out) <- vec
+
+  if(!is.null(para)) {
+  	name <- para
+  	if (para == "SR_LN(R0)") {
+  	  colnames(out) <- paste0("R0 ", vec)
+  	}
+  	if (para == "NatM_uniform_Fem_GP_1") {
+  	  colnames(out) <- paste0("M_f ", vec)
+  	}
+  	if (para == "NatM_uniform_Mal_GP_1") {
+  	  colnames(out) <- paste0("M_m ", vec)
+  	}
+  	if (para == "SR_BH_steep") {
+  	  colnames(out) <- paste0("h ", vec)
+  	}
+  }
+
+  utils::write.csv(x = out, file = file.path(mydir, paste0(name, "_quant_table.csv")), row.names = TRUE)
+}

R/get_retro_quants.R

@@ -0,0 +1,68 @@
+#' Get model values across retrospective runs
+#'
+#' @template mydir
+#' @template model_settings
+#' @param output List of model output created by [run_retro()].
+#'
+#' @author Chantel Wetzel
+#' @return
+#' Nothing is explicitly returned from [get_retro_quants()].
+#' The following objects are saved to the disk.
+#'
+#' * `mohnsrho.csv` with the following columns:
+#'   * type: the type of Mohn's rho
+#'     * [Mohn (1999)](https://academic.oup.com/icesjms/article/56/4/473/624639),
+#'     * Woods Hole Mohn's rho [(Legault 2009)](https://archive.nefmc.org/tech/council_mtg_docs/Sept%202009/Herring/Doc%209_Retro%20Working%20Group%20Report.pdf) used by the [Northeast Fisheries Science Center (NEFSC)](https://www.fisheries.noaa.gov/about/northeast-fisheries-science-center), and
+#'     * [Hurtado-Ferro et al. (2015)](https://doi.org/10.1093/icesjms/fsu198) used by the [Alaska Fisheries Science Center (AFSC)](https://www.fisheries.noaa.gov/about/alaska-fisheries-science-center)
+#'   * Quantity: the stock assessment quantity of interest
+#'   * values: the Mohn's rho values
+#'
+#'   `apply(utils::read.csv(file.path("..", paste0(mod_loc, "_retro"), "retrofigures4doc.csv")), 1, function(x) do.call(sa4ss::add_figure, as.list(x)))`
+#'
+#' * `mohnsrho.tex` for use with `sa4ss::read_child()`
+#'   inside of an environment with `results = "asis"`
+#'   to include a table of Mohn's rho values in a document.
+#'
+#'   `sa4ss::read_child(file.path(paste0(params$model, "_retro"), "mohnsrho.tex"))`
+#'
+#'
+#' @export
+
+get_retro_quants <- function(mydir,  model_settings, output) {
+  retro_dir <- output[["plotdir"]]
+  endyrvec <- output[["endyrvec"]]
+  retroSummary <- output[["retroSummary"]]
+  rhosall <- output[["rhosall"]]
+  rhos <- output[["rhos"]]
+
+  get_param_values(
+    mydir = retro_dir,
+    para = "retro",
+    vec = c("Base Model", paste0("Retro -", 1:(length(endyrvec)-1))),
+    summary = retroSummary
+  )
+
+  utils::write.csv(
+    x = data.frame(
+      caption = paste(
+        "Retrospective patterns for",
+        c("spawning stock biomass (\\emph{SSB})", "fraction unfished"),
+        "when up to", xfun::numbers_to_words(max(abs(model_settings$retro_yr))),
+        "years of data were removed from the base model.",
+        "Mohn's rho (Mohn, 1999) values were",
+        "recalculated for each peel given the removal of another year of data.",
+        "See Table \\ref{tab:RetroMohnsrho} for other derivations of Mohn's rho."
+      ),
+      alt_caption = sprintf("Each successive peel of data led to a Mohn's rho of %s for %s.",
+                            lapply(c("SSB", "Bratio"), function(x) {
+                              knitr::combine_words(sprintf("%.2f", (rhosall[rownames(rhosall) == x, ])))
+                            }),
+                            c("SSB", "fraction unfished")
+      ),
+      label = c("RetroSsb", "RetroFractionunfished"),
+      filein = file.path("..", retro_dir, c("compare2_spawnbio_uncertainty.png", "compare4_Bratio_uncertainty.png"))
+    ),
+    file = file.path(retro_dir, "retrofigures4doc.csv"),
+    row.names = FALSE
+  )
+}